Numerical simulations of the planar contraction flow for a polyethylene melt using the XPP model

J. Non-Newtonian Fluid Mech. 117 (2004) 73–84

Numerical simulations of the planar contraction flow for a polyethylene melt using the XPP model Wilco M.H. Verbeeten, Gerrit W.M. Peters∗ , Frank P.T. Baaijens Materials Technology, Faculty of Mechanical Engineering, Eindhoven University of Technology, P.O. Box 513, 5600 MB Eindhoven, The Netherlands Received 10 February 2003; received in revised form 19 December 2003

Abstract The Discrete Elastic Viscous Stress Splitting Technique in combination with the Discontinuous Galerkin (DEVSS/DG) method is used to simulate a low density polyethylene melt flowing in a transient contraction flow problem. Numerical results using the original and a modified form of the eXtended Pom–Pom (XPP) model are compared to numerical results obtained with the exponential form of the Phan–Thien Tanner (PTT-a) and the Giesekus model and to experimental data of velocities and stresses. These models are known to be well capable of predicting all characteristic features encountered experimentally. Curiously, the eXtended Pom–Pom mode, formulated with the same non-affine or irreversible stretch dynamics as the original Pom–Pom model, encounters convergence problems using the DEVSS/DG method. A slight modification of the stretch dynamics such that it becomes consistent with other viscoelastic models and in agreement with a modification of the stretch dynamics based on non-equilibrium thermodynamics by Van Meerveld [J. Non-Newtonian Fluid Mech. 108 (2002) 291] gives a more numerical stable behaviour and steady state could be reached. From this it is clear that physical and numerical issues still play a mixed role in numerical viscoelastic flow problems. © 2004 Elsevier B.V. All rights reserved. Keywords: DEVSS/DG method; Constitutive models; Differential models; eXtended Pom–Pom model; PTT-a model; Giesekus model; Polymer melt; Contraction flow; Numerical–experimental comparison; Convergence problems

1. Introduction In recent years, the performance of numerical methods for viscoelastic flow analysis has improved significantly. A review on mixed finite element methods for viscoelastic flow analysis is given in Baaijens [3]. Furthermore, constitutive models have become available that are able to quantitatively describe the behaviour of polyethylene melts over a wide range of rheological data, e.g. the molecular stress function (MSF) model by Wagner et al. [29] and the eXtended Pom–Pom model (XPP) [7,14,18,27]. Both a robust numerical method and a reliable constitutive model are crucial in predictive modelling of viscoelastic materials. The robustness and efficiency of the Discrete Elastic Viscous Stress Splitting Technique in combination with the Discontinuous Galerkin method (DEVSS/DG) has been shown in several studies [2,4,5,11]. A drawback of the method is that it is not unconditionally stable in transient calculations ∗

Corresponding author. E-mail address: [email protected] (G.W.M. Peters). URL: http://www.mate.tue.nl/. 0377-0257/$ – see front matter © 2004 Elsevier B.V. All rights reserved. doi:10.1016/j.jnnfm.2003.12.003

at high Weissenberg numbers, as shown by Bogaerds et al. [10]. They found a temporal instability in a simple couette flow for a single-mode upper convected Maxwell (UCM) model depending on the amount of elasticity. However, the DEVSS/DG method has proven to be very robust and accurate in non-smooth problems [2]. This contrary to numerical techniques that apply the SUPG method for handling the convective terms in the constitutive equation. As a constitutive model, the multi-mode eXtended Pom–Pom model is capable of quantitatively describing the rheological behaviour of polyethylene melts over a full range of well-defined rheometric experiments [27,28]. Good agreement is found for transient and steady state shear, including reversed flow, and different elongational flows with this model. More conventional models, such as the Giesekus and exponential Phan–Thien Tanner (PTT-a) models yield less satisfactory results in describing the non-linear behaviour in both shear and elongation simultaneously [4,28]. In a recent study, Bishko et al. [6] calculated the flow of a highly branched polymer in a 4:1 contraction geometry using a single-mode version of the original differential Pom–Pom model. Their results were promising and they predicted in

74

W.M.H. Verbeeten et al. / J. Non-Newtonian Fluid Mech. 117 (2004) 73–84

a good qualitative sense all the specific features observed in experiments for low density polyethylene (LDPE) melts. Verbeeten et al. [28] performed calculations for a commercial grade LDPE melt flowing through two benchmark complex flow geometries, the confined flow around a cylinder and the cross-slot flow, using the XPP model and the DEVSS/DG method. Comparison between numerical and experimental results showed a good qualitative, and to a large extent also quantitative, agreement with the velocities and stresses in the two complex flow geometries. Motivated by the results of Verbeeten et al. [27], Bishko et al. [6] and Verbeeten et al. [28], a contraction flow problem will be investigated. Furthermore, this prototype industrial flow problem is also of interest due to the presence of a corner singularity, which makes it numerically more interesting, and because it is the most widely investigated geometry, closely related to real industrial applications. The objective of this study is to investigate the performance of the multi-mode eXtended Pom–Pom model representing a commercial grade LDPE melt using the DEVSS/DG method. Results will be compared against experimental data taken from Schoonen [22] and the performance of the PTT-a and Giesekus models. Unfortunately, for transient calculations using the combination of the XPP model as given in Verbeeten et al. [27] and the chosen numerical method convergence problems occur for the contraction flow problem at the Weissenberg number calculated in this study. Unphysical negative backbone stretch values are encountered near the re-entrant corner of the contraction, i.e. the geometrical singularity, and along the top wall after the re-entrant corner leading to divergence of the numerical solution. Similar negative backbone stretch values were also found by Verbeeten et al. [28] at the front and back stagnation points in the flow around a cylinder. The contraction flow problem is known to have a thin shear stress boundary layer, situated along the top wall after the re-entrant corner, like the stress boundary layer in the cylinder flow problem. Additionally, the contraction flow problem has a geometrical singularity. Although the DEVSS/DG method is known for not being unconditionally stable, it is, however, difficult to demonstrate if a temporal or physical instability is the cause of divergence, since no analytical solution is at hand. Standard solutions, such as mesh refinement, time step reduction, addition of solvent viscosity, and rounding off the re-entrant corner to eliminate the geometrical singularity did not solve the divergence problem. Van Meerveld [25] looked at the integral Pom–Pom model from a GENERIC point of view. He proposed the small change in the evolution equation for the backbone stretch. By incorporating this modification, a more stable behaviour was achieved and the steady state solution could be reached. In Section 2, the problem definition and the constitutive equations are outlined. The computational method is shortly described in Section 3 and can be found in more detail in Bogaerds et al. [11]. The material characterisation and the

performance of the three constitutive models in simple rheological flows is given in Section 4. Section 5 describes the performance of both the XPP and PTT-a model in the transient contraction flow problem. The conclusions are given in Section 6. Depending on the used model, its structure, and its strain hardening behaviour, calculations with the DEVSS/DG method can reach the steady state solution.

2. Problem definition Isothermal and incompressible fluid flows, neglecting inertia, are described by the equations for conservation of momentum (1) and mass (2):   · T = 0, ∇

(1)

 ·u ∇  = 0,

(2)

 the gradient operator, T the Cauchy stress tensor and with ∇ u  the velocity field. The Cauchy stress tensor T is defined for polymer melts by: T = −pI + 2ηs D +

M


τi.

(3)

i=1

Here, p is the pressure, I the unit tensor, ηs denotes the viscosity of the purely viscous (or solvent) contribution, D = 1/2(L + LT ) the rate of deformation tensor, in which u L = (∇  )T is the velocity gradient tensor and (·)T denotes the transpose of a tensor. The visco-elastic contribution of the ith relaxation mode is denoted by τ i and M is the total number of different modes. The multi-mode approach is necessary for most polymeric fluids to give a realistic description of stresses over a broad range of deformation rates. The visco-elastic contribution τ i still has to be defined by a constitutive model. 2.1. Constitutive models Within the scope of this work, a sufficiently general expression for the extra stress tensor of the ith relaxation mode is given by: ∇

τ i + Bi · τ i + τ i · BTi + Gi (Bi + BTi ) = 2Gi (1 − ξi )D, (4) with Bi the averaged, macroscopic slip tensor, Gi the plateau modulus, and ξi a slip parameter that, if non-zero, accounts for the Gordon–Schowalter derivative. The upper convected time derivative of the extra stress τ i ∇ is defined as: ∇

τ i = τ˙ i − L · τ i − τ i · LT ∂τ i  i − L · τ i − τ i · LT , = +u  · ∇τ ∂t

(5)

with t denoting time. The slip tensor Bi defines the chosen constitutive equation.

W.M.H. Verbeeten et al. / J. Non-Newtonian Fluid Mech. 117 (2004) 73–84

The eXtended Pom–Pom model [27] is a modification of the differential version of the original Pom–Pom model by McLeish and Larson [18] incorporating branch-point withdrawal as introduced by Blackwell et al. [7]. The basic idea of this model is that the rheological properties of entangled polymer melts depend on the topological structure of the polymer molecules. The simplified topology consists of a backbone tube with a number of dangling arms at both ends and is called a Pom–Pom molecule. A key feature is the separation of relaxation times for the stretch and orientation. For the Pom–Pom model, in which the slip parameter is not used and equals ξi = 0, the extra stress tensor is expressed as: τ i = σ i − Gi I = 3Gi Λ2i S i − Gi I,

(6)

with Λi the backbone tube stretch, defined as the length of the backbone tube divided by the length at equilibrium, and S i denoting the orientation tensor. The evolution equations for the orientation tensor S i and backbone tube stretch Λi are given by: ∇

S i + 2[D : S i ]S i + Bi · S i + S i · BTi − 2[Bi : S i ]S i = 0, ˙ i = Λi [D : S i ] − Λi [Bi : S i ]. Λ

75

result was found by Van Meerveld [25] who analyzed the integral Pom–Pom model using the thermodynamical framework GENERIC. He pointed out that it may be of importance for contraction/expansion flows, in which Wapperom and Keunings [30] detected a negative affine motion. Physically, this change in the stretch dynamics means a faster relaxation towards its equilibrium state when the backbone is stretched or compressed. This could be favourable in circumstances that the stretch increases very fast, as is the case during uniaxial elongation and close to geometrical singularities. In fact we have detected that this seemingly small change in the constitutive equation marks the difference between reaching converged results and divergence of the numerical scheme. The slip tensor Bi for the modified XPP model, i.e. including the modified stretch dynamics, is then given by:  1 − αi − 3αi Λ4i Tr(S i · S i ) αi Bi = σi + 2Gi λb,i 2λb,i Λ2i   1 Gi (1 − αi ) −1 ri eνi (Λi −1) 1− 2 I− + σi . λb,i 2λb,i Λi

(7)

(10)

(8)

Due to its different structure from most differential constitutive equations and non-zero initial conditions, it may be cumbersome to implement the above double-equation formulation consisting of Eqs. (7), (8) and (10) in FE packages. However, the XPP model can also be written in a fully equivalent single-equation fashion by substituting the slip tensor of Eq. (10) into Eq. (4). A numerically consistent implementation of the backbone stretch values are then calculated from Eq. (6):

In the original Pom–Pom model [18] and in the eXtended Pom–Pom model introduced by Verbeeten et al. [27], the slip tensor Bi is given by:  1 − αi − 3αi Λ4i Tr(S i · S i ) αi σi + Bi = 2Gi λb,i 2λb,i Λ2i    1 Gi (1 − αi ) −1 ri eνi (Λi −1) 1− I− + σ i . (9) λb,i Λi 2λb,i Here, αi is an adjustable parameter, λb,i the relaxation time of the backbone tube orientation equal to the linear relaxation time λi , and σ i is defined according to Eq. (6). Both the plateau modulus Gi and the relaxation time λi are obtained from the linear relaxation spectrum determined by dynamic measurements, while αi controls the anisotropic drag, similar to the Giesekus model. ri = λb,i /λs,i is the ratio between the orientation relaxation time and the stretch relaxation time λs,i , and νi = 2/qi is a parameter denoting the influence of the surrounding polymer chains on the backbone tube stretch with qi the amount of arms at the end of a backbone. The expression (9) for the slip tensor Bi was chosen such that the non-affine or irreversible stretch dynamics to be proportional to [Λi − 1]. This is determined only by the term with λs in Eq. (9) and this term does not contribute to the orientation evolution, Eq. (7). However, this choice is not consistent with the irreversible stretch dynamics as found for the other viscoelastic models such as the PTT-a, Giesekus and UCM model (see the end of this section for the formal expressions). To make it consistent the term [Λi − 1] should be replaced by [Λi − (1/Λi )]. This introduces a singularity for Λi = 0 and guarantees that Λi > 0. The same

Λ2i = 1 + and Λi =

Tr(τ i ) . 3Gi

 1+

|Tr(τ i )| ≥ 1, 3Gi

(11)

(12)

which is only present in the exponential term and secures its minimum at Λi = 1. Notice, that for this definition of the backbone stretch, we avoid numerical problems if Tr(τ i )/3Gi < −1. Theoretically, Tr(τ i ) is always 0 or positive. However, numerically we have encountered these unrealistic values near the re-entrant corner in the contraction flow and at the front and back stagnation points in the flow around a cylinder for high enough flow rates, both for the XPP model as well as the Giesekus and PTT-a models. For the Giesekus model, ξi = 0 and the tensor Bi is defined as: αi (1 − 2αi ) Gi (1 − αi ) −1 Bi = σi + I− σi , (13) 2Gi λi 2λi 2λi with αi an adjustable parameter, which controls the anisotropic drag. Notice that for αi = 0, the upper convected Maxwell model is retrieved.

76

W.M.H. Verbeeten et al. / J. Non-Newtonian Fluid Mech. 117 (2004) 73–84



For the exponential PTT model, Bi is defined as: Bi = −ξi D +

e(εi Tr(τ i )/Gi ) Gi e(εi Tr(τ i )/Gi ) −1 I− σi . 2λi 2λi

¯ + Dv , 2ηs D + 2¯η(D − D)

˙ i = Λi [(1 − ξi )D : S i ] − Λ

eεi Tr(τ i )Gi 2λi

 Λi −

1 Λi

τ i − (∇ · v, p) = 0, (18)

 .



i=1

(14)

with εi an adjustable parameter, controlling the amount of strain softening. Similar to the XPP model, both the Giesekus and PTT-a models can be split up in an orientation and stretch equation. Substitution of Eq. (13) (for the Giesekus model) or Eq. (14) (for the PTT-a model) into the Eqs. (7) and (8) result in the evolution equations for the backbone orientation and backbone stretch. To illustrate the differences and similarities between the three models, the backbone stretch equations are shown for the Giesekus:  1 ˙ Λi = Λi [D : S i ] − 3Λ2i αi Tr(S i · S i ) 2λi  (1 − αi ) , (15) +(1 − 2αi )Λi − Λi and for the PTT-a model:

M


(16)

If Eq. (11) is used, the expression for the PTT-a model can also be written as:   ∗ 2 eνi (Λi −1) 1 ˙ i = Λi [(1 − ξi )D : S i ] − Λ Λi − , 2λi Λi νi∗ = 3εi . (17) Notice that the stretch relaxation for the Giesekus and PTT-a models are also proportional to [Λi − (1/Λi )]. 3. Computational method The DEVSS/DG method, Baaijens et al. [4], is applied. This method efficiently handles multiple relaxation modes and is robust near geometrical singularities. It is a combination of the Discrete Elastic Viscous Stress Splitting Technique of Guénette and Fortin [12] and the Discontinuous Galerkin method, developed by Lesaint and Raviart [16]. Application to the single-equation XPP model is discussed in this section. The weak formulation and solution strategy for the Giesekus and PTT-a models is similar and can be found in more detail in Bogaerds et al. [11] and Baaijens et al. [4]. 3.1. DEVSS/DG method Consider the problem described by the Eqs. (1), (2), (4) and (10). Following common FE techniques, these equations are converted into a mixed weak formulation: 3.1.1. Problem DEVSS/DG ¯ τ i ) for any time t such that for all admisFind ( u, p, D, sible test functions (v, q, E, si ),

(q, ∇ · u  ) = 0,

(19)

¯ − D) = 0, (E, D

(20)



   αi 1 1 si , τ i + τi · τi + 2ri eνi (Λi −1) 1 − 2 λb,i Gi λb,i Λi  αi Tr(τ i · τ i ) 1 τi + 2− Λi 3G2i Λ2i    1 1 1 νi (Λi −1) 2ri e + 1− 2 + 2 λb,i Λi Λi   αi Tr(τ i · τ i ) − − 1 Gi I − 2Gi D 3G2i Λ2i K
− si : u  · n (τ i − τ ext i )dΓ = 0 ∇

e e=1 Γin

∀ i ∈ {1, 2, . . . , M}.

(21)

Here, (·, ·) denotes the L2 -inner product on the domain Ω,  v+(∇  v)T ), η¯ an auxiliary viscosity chosen to be Dv = 1/2(∇ M η¯ = i=1 Gi λi , Γine is the inflow boundary of element Ωe , n the unit vector pointing outward normal on the boundary of the element (Ωe ) and τ ext i denote the extra stress tensor of the neighbouring, upstream element. Time discretisation of the constitutive equation is attained using an semi-implicit Euler scheme. In this scheme most of the variables are updated implicitly with the exception of the ext external component τ ext i which is taken explicitly (i.e. τ i = ext τ i (tn ). Furthermore, the equation for the discrete rate of ¯ is decoupled from the ‘Stokes’ problem and deformation D thus updated explicitly. 3.2. Solution strategy A two-dimensional flow domain is divided into quadrilateral elements as a spatial discretisation. The interpolation of the velocity field u  is bi-quadratic, the pressure p ¯ are interpolated and discrete rate of deformation fields D bi-linearly, while the extra stress tensor τ i field has a discontinuous bi-linear interpolation. The Eqs. (18)–(21) are integrated over elements using a 3 × 3-Gauss quadrature rule. The non-linear equations are solved using a one step Newton–Raphson iteration. Within the linearized set of equations, the extra stress related variable τ i forms a block diagonal structure due to the fact that the components in the upstream elements are taken explicitly. This enables elimination of the extra stress variable at element level. A further reduction of the system matrix is obtained by a decoupling

W.M.H. Verbeeten et al. / J. Non-Newtonian Fluid Mech. 117 (2004) 73–84

77

with N being the total number of nodes in the mesh. For calculations presented in this paper, the convergence criterium is chosen to be = 10−7 .

approach. First, the ‘Stokes’ problem ( u, p) is solved, after which the updated solution is used to find a new approx¯ Notice, that this introduces a small extra imation for D. ¯ in the inconsistency for the stabilization term 2¯η(D − D) ¯ momentum Eq. (18), since D is taken at time tn+1 , while D is taken at time tn . Finally, the nodal increments of the orientation tensor and backbone stretch are retrieved element by element. To solve the non-symmetrical system of the ‘Stokes’ problem ( u, p) an iterative solver is used based on the Bi-CGSTAB method of Van der Vorst [24]. The symmet¯ is solved rical set of algebraic equations of problem (D) using a conjugate-gradient solver. Incomplete LU preconditioning is applied to both solvers. In order to enhance the computational efficiency of the Bi-CGSTAB solver, static condensation of the centre-node velocity variables results in filling of the zero block diagonal matrix in the ‘Stokes’ problem and, in addition, achieves a further reduction of global degrees of freedom. The above decoupling procedure enables the use of iterative solvers. The use of the Bi-CGSTAB solver in combination with the incomplete LU ¯ was unsuccessful. preconditioning on the problem ( u, p, D) For more details, we refer to [11,28]. To solve the above sets of algebraic equations, both essential and natural boundary conditions must be imposed on the boundaries of the flow channels. At the entrance of the flow channel the velocity profile is prescribed. At the exit, either the velocity profile is prescribed together with a pressure point that is set to zero, or the velocities perpendicular to the outflow directions are suppressed. At the entrance, the values of the extra stress tensor of a few elements downstream the inflow channel are prescribed along the inflow boundary, imposing a periodic boundary condition for the stress variable. Calculations were stopped when the variables had reached steady state values. The steady state convergence criterium is established by calculating the least square difference between the solution of the actual and previous iteration step, divided by the least square actual solution, and checking if it is smaller than a certain convergence criterium : 

N n+1 − uni )2 i=1 (ui 

≤ , (22) N n+1 2 ) i=1 (ui

4. Material characterisation Since a multi-mode approach is necessary for a realistic description of polymer melts, we will investigate multi-mode calculations. Parameters are identified for a commercial grade polymer melt, a low density polyethylene (DSM, Stamylan LD 2008 XC43), subsequently referred to as LDPE. This long-chain branched material has been extensively characterised by Schoonen [22] and is also given in Peters et al. [21]. The parameters for a four mode Maxwell model are obtained from dynamic measurements at a temperature of T = 170 ◦ C. The non-linearity parameters qi and the ratio ri = λb,i /λs,i for the XPP model and αi for the Giesekus model are determined using the transient uniaxial elongational data [31] only. For the PTT-a model εi and ξi are identified by using transient uniaxial elongational and shear data. Since second planar elongational or second normal stress difference data is not available for this material, the anisotropy parameter αi in the XPP model is chosen as 0.3/qi , similar to Verbeeten et al. [28]. All parameters are given in Table 1. Fig. 1 shows the rheological behaviour of the three constitutive models and measured data in uniaxial and simple shear behaviour. The ‘steady state’ data points in the uniaxial elongational viscosity plot are the end points of the transient curves. Since the samples break at the end of the transient measurements, it is hard to say if steady state has truly been reached. However, sometimes it is claimed that steady state coincides with break-up of the samples [17]. All three models can satisfactorily predict the uniaxial elongational data in a transient sense. However, notice that the PTT-a model is showing an upswing towards the quasi-steady state values that is a little too soon in time. Furthermore, for the high elongation rate ε˙ = 30 s−1 , for which no experimental data is available, the Giesekus model shows a much higher upswing compared to the other two models. For the quasi-steady state data, the XPP model predicts the most smooth and realistic curve. With this set of relaxation times and these parameters, the PTT-a model demonstrates a more non-smooth behaviour, reflecting the

Table 1 Linear and non-linear parameters for fitting of the DSM Stamylan LD 2008 XC43 LDPE melt i

1 2 3 4

Maxwell parameters

mXPP

Giesekus

Exp. PTT

Gi (Pa)

λi (s)

qi

λb,i /λs,i

αi

αi

εi

ξi

7.2006 × 104 1.5770 × 104 3.3340 × 103 3.0080 × 102

3.8946 × 10−3 5.1390 × 10−2 5.0349 × 10−1 4.5911 × 100

1 1 2 10

1.0 1.0 1.0 1.0

0.30 0.30 0.15 0.03

0.30 0.30 0.25 0.04

0.30 0.20 0.02 0.02

0.08 0.08 0.08 0.08

T = 170 ◦ C. νi = 2/qi . Activation energy: E0 = 48.2 kJ/mol. WLF-shift parameters: C1 = 14.3 K, C2 = 480.8 K.

78

W.M.H. Verbeeten et al. / J. Non-Newtonian Fluid Mech. 117 (2004) 73–84

Fig. 1. Transient (a) and quasi-steady state (b) uniaxial viscosity ηu , transient (c) and steady state (d) shear viscosity ηs , and transient (e) and steady state (f) first normal stress coefficient Ψ1 at T = 170 ◦ C of the mXPP, Giesekus and PTT-a models for DSM Stamylan LD2008 XC43 LDPE melt.

individual relaxation times due to the excessive elongational thinning behaviour of the model at high strain rates. This non-smooth behaviour could be reduced by incorporating more modes, however at the cost of a significant increase of the computational time. As expected and also indicated by the transient model predictions at the high elongation rate, the Giesekus model predicts the wrong shape of the steady state curve, as it ends in a plateau for high strain rates. Low density polyethylene melts are in general elongational thinning in a steady state sense at high strain rates (see for

example data from Münstedt and Laun [19] and Hachmann [13]). For simple shear, all models show a rather good agreement with the measured data (Fig. 1(c) and (d)). For intermediate shear rates, the Giesekus and PTT-a models somewhat over predict the measurements. At shear rates over 102 s−1 , the curves for all models drop rapidly (see Fig. 1d), as here only the last mode of the relaxation spectrum is active, showing a too steep slope of −1. For the contraction flow problem at the chosen Weissenberg number, using the definition for the

W.M.H. Verbeeten et al. / J. Non-Newtonian Fluid Mech. 117 (2004) 73–84

y/h

0

H h

flow direction 0

x/h

Fig. 2. Schematic of half the planar contraction flow. The origin of the coordinate system is at the centre of the contraction plane. T = 170 ◦ C, λ¯ = 1.7415 s, u¯ 2D = 7.47 mm/s, h = 0.775 mm, H = 2.55 mm, d = 40 mm.

√ effective strain rate ε˙ eff = 1/2 2D : D, maximum shear and strain rates of γ˙ max ≈ 50 s−1 and ε˙ max ≈ 40 s−1 were observed around the re-entrant corner. Thus, calculations all fall within the range of the chosen relaxation spectrum capable of predicting the experimental data. In order to obtain a good fit of the shear data for the PTT-a model, the non-linearity parameter ξ is larger than ε for the two longest relaxation times. Although difficult to observe in Fig. 1(c), this results in oscillations in transient start-up of simple shear (see also Appendix A in Verbeeten et al. [28]). For the first normal stress difference (Fig. 1(e) and (f)), the XPP model shows the best agreement for intermediate shear rates, while the PTT-a model performs better at low and high shear rates. Note that the experimental data keeps on decreasing at larger times. This is most probably due to instabilities at the free surface in the cone-plate experiments, see for example Pahl et al. [20]. Although similar discrepancies between the experimental data and the calculated results exist for all models, in our opinion, the XPP model shows the best overall agreement with experimental rheological data. This opinion is formed mostly due to the most smooth and realistic curve for the steady state elongational viscosity and its capability for simultaneously describing the shear and elongational experimental data of other low density polyethylene melts in quantitative sense.

5. Planar contraction flow problem The planar contraction benchmark flow is known to have regions with pure simple shear flow, pure planar elongation and combinations thereof. It also has a corner singularity which makes it numerically challenging. A schematic of half the flow geometry and its characteristics are given in Fig. 2. Experiments were performed at a temperature of T = 170 ◦ C, resulting in a viscosity relaxation time for the material of averaged



M M 2G / λ¯ = λ i=1 i i i=1 λi Gi = 1.7415 s. The downstream and upstream half heights of the geometry are h = 0.775 mm and H = 2.55 mm, respectively, giving an actual contraction ratio of 3.29. The average two-dimensional mean upstream velocity is equal to U¯ 2D = 2.27 mm/s, and

79

consequently the average downstream velocity is u¯ 2D = 7.47 mm/s. The depth of the flow cell equals d = 40 mm to create a nominally two-dimensional flow, resulting in a depth-to-height ratio of 7.84 upstream and 25.81 downstream. This depth-to-height ratio reduces the influence of the front and back walls to about 6% maximum on the isochromatic stress patterns with the flow rate investigated in this work (see [11,23]). This is a sufficiently small deviation to compare experimental results to two-dimensional calculations. By taking half the downstream channel height and the downstream two-dimensional mean velocity as characteristic values, the flow averaged dimensional Weissenberg number, indicating the shear strength of the flow, results in: Wi =

λ¯ u¯ 2D = 16.8. h

(23)

The maximum averaged Weissenberg number of the flow, using the longest relaxation time, is: Wimax =

λmax u¯ 2D = 44.3. h

(24)

If related to the local shear rate γ˙ = ∂u/∂x, higher Weissenberg numbers can even be found in the flow. Particle tracking velocimetry (PTV) was used to obtain the experimental velocities. Stresses were measured using the field-wise flow induced birefringence (FIB) method recording isochromatic fringe patterns. To link these fringe patterns to calcuted stresses, we have used the semi-emperical stress optical rule. For two-dimensional flows, the stress optical rule simplifies to:  kϑ 2 + N2 = τFRG = 4τxy , (25) 1 dC with τFRG defined as the isochromatic fringe stress, τxy the shear stress in the xy-plane, N1 = τxx − τyy the first normal stress difference, ϑ the wave length of the light used in the measurements, d the light path length in the birefringent medium, which in this case is the depth of the flow cell, C = 1.47 × 10−9 Pa−1 the stress optical coefficient (identified by Peters et al. [21]), and k the fringe order of the observed dark fringe bands where extinction of the light occurs. For more details on the experimental aspects and more experimental results, we refer to Schoonen [22]. 5.1. Original XPP model For the original eXtended Pom–Pom mode, as given by Verbeeten et al. [27,28], convergence problems occurred and no steady state results could be reached. It was difficult to detect the exact reason for the onset of divergence, as various possibilities interact at the same time. Divergence problems may be accounted for by a temporal instability of the DEVSS/DG method as found in simple couette flows for a single-mode UCM model by Bogaerds et al. [10]. Moreover, it could be due to the highly strain

80

W.M.H. Verbeeten et al. / J. Non-Newtonian Fluid Mech. 117 (2004) 73–84 Table 2 Mesh characteristics of the planar contraction flow

hardening behaviour or anisotropy with the chosen set of parameters. The onset of divergence may have been provoked by Tr(τ i )/3Gi < −1 resulting in negative backbone stretch values near the re-entrant corner singularity. This occurrence was detected for both the original XPP model as well as for the Giesekus and PTT-a models. However, for the PTT-a model the negative backbone stretch values disappeared after sufficient iteration steps and steady state results could be reached. This contrary to the Giesekus and the original XPP model, where negative backbone stretch values were detected well before the onset of divergence. These negative backbone stretch values were only encountered for the mode with the largest relaxation time. Various attempts were made to overcome the convergence problems for the contraction flow problem. Mesh refinement, decreasing the time step, changing the boundary conditions, or adding a little solvent viscosity at best delayed the onset of divergence. Reaching full convergence within each time step using multiple Newton iterations instead of only one, did not resolve the convergence problems either. By choosing α = 0, which means no anisotropy, calculations diverged even sooner, which rejects the idea that anisotropy has a negative effect on the numerical stability. The influence of a rounded corner was investigated to see if the corner singularity was the cause of divergence. Surprisingly, by rounding of the re-entrant corner divergence occurred even sooner in time. Higher stresses and stress gradients appeared in the simulations and that may influence convergence in a negative sense. Although fully equivalent to the single-equation formulation, the double-equation formulation of the original eXtended Pom–Pom model shows different numerical response. Unfortunately, for the double-equation formulation, negative backbone stretch values were also encountered, although they appeared only just before onset of divergence. Possibly the convergence problem is caused by an insufficiently fine mesh or too large time step. Hence, the numerical scheme may not be capable of coping with the high stress gradients due to the high strain hardening behaviour of

#Elements #Nodes #DOF( u, p) ¯ #DOF(D) #DOF(S,.)=4 × #elements× 20

M2

2.5

2

2 y/h

y/h

2.5

1.5

1.5

1

1

0.5

0.5

0 x/h

(a)

1246 5189 9235 4047 99680

A comparison is made between the experimental data and numerical results for the modified eXtended Pom–Pom (mXPP) and PTT-a models. The Giesekus model was not able to reach steady state. We assume this to be caused by predicting a steep upswing and high uniaxial elongational viscosity for high strain rates, which is significantly higher than for the other two models, as is depicted in Fig. 1a. First, mesh refinement has been investigated for the mXPP model using two different meshes. Only half of the geometry is meshed due to symmetry. Details of the FE meshes are depicted in Fig. 3 having mesh characteristics given in Table 2. The contraction is placed at x/ h = 0, while the inlet is at x/ h = −20 and the outlet at x/ h = 30. The periodic boundary for the inlet stresses is situated at x/ h = −10. At the entrance and exit, a fully developed velocity profile is prescribed. No slip boundary conditions are imposed on the

3

-1

1000 4161 7403 3243 80000

5.2. Modified XPP model

3

-2

Mesh M2

the constitutive models. However, finer meshes and smaller time steps result in a significant and undesirable increase of computational time. The influence of the entrance and exit lengths and the auxiliary viscosity η¯ also remains to be investigated. A clearly more favourable solution is the incorporation of the modified stretch dynamics in the XPP model. It does not affect the computational effort and seems to influence numerical stability in a positive way. By including the modified stretch dynamics, steady state results could be reached, which will be presented here.

M1

0 -3

Mesh M1

1

2

0 -3

-2

-1

0 x/h

(b)

Fig. 3. Detail of FE meshes of half the planar contraction flow: (a) Mesh M1; (b) Mesh M2.

1

2

W.M.H. Verbeeten et al. / J. Non-Newtonian Fluid Mech. 117 (2004) 73–84 4

4

x 10

Characteristic stress values along centerline Wi = 16.8 mesh M1 mesh M2

[Pa]

3.5

N21 + 4τ  √ 2xy

top confining wall, and the velocity in y-direction is suppressed for the symmetry line y/ h = 0. Mesh M1 has the elements more or less aligned in the flow direction, using a three-quarter division circle to mesh around the corner singularity. The elements are gradually refined towards the re-entrant corner, resulting in a very fine mesh near the singularity. Due to the corner singularity and a thin stress layer at the upper wall after the re-entrant corner, this is a numerically difficult region. Therefore, mesh M2 has refinement at the downstream region after the re-entrant corner. Fig. 4 shows that mesh convergence has been reached with these two meshes. For the mXPP and PTT-a models using mesh M2 and time step 0t = 0.001 s, the velocities and stresses are compared to experimental data in Figs. 5–7. As can be seen in Fig. 5, in the upstream region before the contraction, both models are in excellent agreement with the experimental data. In the contraction region, agreement between calculated and experimental results is very good, although the PTT-a model is somewhat overpredicting the experimental data. In the downstream region after the contraction, both models are overpredicting the experimental

81

3 2.5 2 1.5 1 0.5 0 -5

0

5

10

x/h

Fig. 4. Calculated isochromatic fringe stresses in the contraction flow over the centerline for the meshes M1 and M2 using the mXPP model at Wi = 16.8. T = 170 ◦ C.

velocity results, especially over the centerline, however agreement is still good. The highest overshoot is predicted by the PTT-a model. Both models predict the same constant maximum velocity in the downstream channel far away Velocities over centerline 12

Dimensionless velocities

Wi = 16.8 Experiment XPP exp. PTT

10 8 u [mm/s]

y/h

2

0

6 4

-2 -6

-4

-2

0 x/h

2

(a) Cross-section velocities

4

6

Wi = 16.8 XPP exp. PTT

2 0 -6

-4

-2

0 x/h

2

4

6

(b) Centerline velocity

Fig. 5. Calculated and measured velocity profiles in the planar contraction flow over four cross-sections (a) (x/ h = −6, −4, −1, 1), and over the centerline (b) (y/ h = 0) for the mXPP and PTT-a models at Wi = 16.8. T = 170 ◦ C.

Fig. 6. Calculated and measured isochromatic fringe patterns of the contraction flow for the mXPP model (a) and PTT-a model (b) at Wi = 16.8. T = 170 ◦ C.

82

W.M.H. Verbeeten et al. / J. Non-Newtonian Fluid Mech. 117 (2004) 73–84

x 10

4 Characteristic

stress values along centerline

5

3 Wi = 16.8 Experiment XPP exp. PTT

[Pa]

4

x 10

Characteristic stress values at y/h=1 Wi = 16.8 Experiment XPP exp. PTT

2.5 [Pa]

5

2

N21 + 4τ √  2xy

N21 + 4τ √  2xy

3

2

1

0

1.5 1 0.5

-5

0

5

0 -5

10

-4

-3

(a) Centerline (y / h = 0) 5

Characteristic stress values at cross section x=0 Wi = 16.8 Experiment XPP exp. PTT

[Pa]

2.5

N21 + 4τ √  2xy

4

15

2 1.5 1

[Pa]

x 10

-1

0

(b) Cross-section (y / h = 1)

10

N21 + 4τ √  2xy

3

-2 x/h

x/h

5

x 10

Characteristic stress values at x/h=4.8 Wi = 16.8 Experiment XPP exp. PTT

0.5 0 -1

-0.5

0 y/h

0.5

1

(d) Cross-section (x / h = 0)

0 -1

-0.5

0 y/h

0.5

1

(d) Cross-section (x / h = 4.8)

Fig. 7. Calculated and measured stresses in the planar contraction flow over the cross-sections (y/ h = 0) (a), (y/ h = 1) (b), (x/ h = 0) (c), and (x/ h = 4.8) (d) for the mXPP and PTT-a models at Wi = 16.8. T = 170 ◦ C.

from the contraction region, which is slightly higher than the experimental data. This overprediction can be ascribed to the higher aspect-ratio of the downstream channel. Given a certain flow-rate in a duct, the mean velocity in the plane of symmetry will relatively become lower for increasing aspect-ratio [22]. In the full field isochromatic fringe patterns (Fig. 6), the top parts are the experimental results, while the bottom halfs are the computational solutions. Agreement between numerical and experimental results is remarkably good: both models predict the ‘butterfly’ shaped fringes and the recirculation zones upstream of the contraction as experimentally observed. Note that some other contraction flow studies, see e.g. Ahmed et al. [1], Kiriakidis et al. [15], predict similar ‘butterlfy’ patterns but do not observe these in the experiments; this is due to the poor aspect-ratio in the upstream region in these studies (≈1), where the flow can not be considered 2D anymore. This has already been mentioned by Schoonen [22]. The PTT-a model is excellently predicting the fringe patterns both in number and position, while the mXPP model is slightly underpredicting the stresses in that region. In the downstream region, the mXPP

model shows a better agreement with the experimental data than the PTT-a model. A more quantitative comparison of the stresses is given in Fig. 7. For the stresses along the centerline (Fig. 7(a)), both models are underpredicting the measured stresses. Here, the PTT-a model shows the best agreement. Although the mXPP model shows the best description for the uniaxial elongational data, it shows the least accurate prediction over the centerline. The discrepancy between experiments and simulations can be accounted for by a number of reasons. First, the parameters of the model are determined on a logarithmic scale mostly using the uniaxial viscosity data. Given the difficulty of performing these experiments, the experimental error is rather large, which also introduces an error in fitting the model parameters. By plotting the stresses over the centerline on a linear scale, possible errors will become even more evident. Second, in the simulations, the transient start-up elongational behaviour possibly lacks behind in time. The results shown here are steady state results, however, a material element flowing over the centerline experiences deformation in a transient sense rather than a steady state deformation. Hence, the transient behaviour

W.M.H. Verbeeten et al. / J. Non-Newtonian Fluid Mech. 117 (2004) 73–84

of the consitutive models is of importance in this complex flow. And third, the three-dimensional effects are most pronounced in this region as velocities increase. To draw more confident conclusions regarding these last effects, full three-dimensional calculations would provide more insight. Over the line y/ h = 1 towards the contraction singularity (Fig. 7(b)), the mXPP model shows the best agreement. It can be seen that both models encounter oscillations near the geometrical singularity, which seems to be more pronouned for the PTT-a model. Possibly more mesh refinement in that region could avoid these oscillations. However, that would also mean an increase in computational time, which is already quite considerable. For the cross-section x/ h = 0 and x/ h = 4.8 (Fig. 7(c) and (d)), the PTT-a model shows a better agreement near the centerline, while the mXPP model predicts the experimental data better away from the centerline near the upper walls. The PTT-a model shows a ‘long-wave’ sort of oscillation towards the upper confining walls. Since this is a shear dominated region, these oscillations may be ascribed to the oscillations in the start-up of transient shear flows for this model (see Fig. 1(c) and (e)). However, further investigation is necessary regarding this phenomenon to underline this statement with more confidence. In general, both the mXPP and PTT-a models show a good to excellent agreement between the numerical and experimental results, both for the velocity and stress data.

6. Discussions and conclusions The transient flow in a contraction problem is investigated using the DEVSS/DG method in combination with the eXtended Pom–Pom model with the modified stretch dynamics. Steady state numerical results using a sufficiently fine mesh and time step 0t = 0.001 s are compared to the velocity and stress experimental data and the well-known conventional exponential form of the Phan–Thien Tanner model. Results of both constitutive models are in good agreement with the experimental data and all the typical experimental features are also predicted by the numerical models. Although the mXPP model shows a good agreement with uniaxial viscosity data, experimental stresses over the centerline are surprisingly underpredicted, which is also the case for the PTT-a model. This discrepancy is expected to arise from three-dimensional effects. More investigation is needed to confirm that assumption. Numerical simulations using the Giesekus model showed convergence problems and no steady state results could be reached. This was also encountered for the for the original XPP model. Negative backbone stretch values were encountered near the re-entrant corner singularity. It is difficult to point out the exact cause of the onset of divergence. It may be due to the locally very high stress gradients which the mesh is unable to resolve, given a certain time step. On the other hand, the temporal instability of the DG method as

83

was found by Bogaerds et al. [10] may influence divergence of the problem. Recently, Bogaerds et al. [8] showed that instabilities are also dependent on the shape of the elements in the mesh; long shaped elements in the flow direction increase stability. No numerical remedy was encountered to overcome divergence of the calculations. Mesh refinement, decreasing the time step, or changing the boundary conditions at best delayed divergence. By increasing 0t and having more iterations per time step, or adding solvent viscosity divergence was encountered even sooner. Rounding off the re-entrant corner to remove the geometrical singularity also caused more difficulties in resolving the problem. Loss of convergence was probably encountered sooner because of an enhanced convection into the downstream flow channel of high stresses imposed by the re-entrant corner. Simulations with both the single-equation as well as the double-equation version of the eXtended Pom–Pom model were performed. Although fully equivalent, they show different numerical responses. For the double-equation version, negative backbone stretch values are only encountered just before onset of divergence. Since the backbone stretch values in the single-equation version are derived from Tr(τ i ) instead of calculated directly, negative values are encountered already earlier in the calculations. This leads to complex backbone stretch values and numerical inconsistencies. The Giesekus and PTT-a models also suffer from this numerical artefact, which for the Giesekus model has led to divergence of the numerical solution. For transient calculations, the DEVSS/DG method becomes more unstable at higher elasticities (see [10]), i.e. higher flow rates in the flow geometry or higher stresses as a consequence of higher strain hardening behaviour. This is consistent with observations of Verbeeten [26] that calculations with the single-equation formulation converge if a parameter set is chosen that describes a lower strain hardening behaviour of the material, and a steady state solution is reached. This holds both for a geometry with a corner singularity as well as a rounded re-entrant corner. However, to our surprise, steady state solutions for the double-equation formulation were never reached, no matter if a low or high strain hardening parameter set was chosen. Most probably, the extra stress variables Λi and S i stay about the same order in magnitude, hardly influenced by a change of material parameters. In the presented solution method, the ‘Stokes’ problem ¯ are ( u, p) and the discrete rate of deformation problem (D) ¯ in decoupled. As a result, the stabilization term 2¯η(D − D) the momentum equation((18)) is inconsistent: D is taken at ¯ is taken at time tn . This decoupling may time tn+1 , while D influence stability of the numerical scheme. By choosing a ¯ can different preconditioner, the coupled problem ( u, p, D) be solved while still using iterative solvers. The DEVSS/DG method shows divergence problems, depending on the constitutive model, the mesh, the time step, and the material behaviour. More stable methods may

84

W.M.H. Verbeeten et al. / J. Non-Newtonian Fluid Mech. 117 (2004) 73–84

improve convergence when solving complex flows. Recently, Bogaerds et al. [9] proposed a new time marching scheme for the time dependent simulation of viscoelastic flows governed by constitutive equations of the differential form. This method does not show the temporal instability problems in simple shear flows which are encountered for the DEVSS/DG method. Moreover, it uses a different preconditioner which is able to solve the coupled problem ¯ However, since it is a SUPG-based method, it is ( u, p, D). less robust for non-smooth flow geometries.

References [1] R. Ahmed, R.F. Liang, M.R. Mackley, The experimental observation and numerical prediction of planar entry flow and die swell for molten polyethylenes, J. Non-Newtonian Fluid Mech. 59 (1995) 129–153. [2] F.P.T. Baaijens, An iterative solver for the DEVSS/DG method with application to smooth and non-smooth flows of the upper convected Maxwell fluid, J. Non-Newtonian Fluid Mech. 75 (1998) 119–138. [3] F.P.T. Baaijens, Mixed finite element methods for viscoelastic flow analysis: a review, J. Non-Newtonian Fluid Mech. 79 (1998) 361– 386. [4] F.P.T. Baaijens, S.H.A. Selen, H.P.W. Baaijens, G.W.M. Peters, H.E.H. Meijer, Viscoelastic flow past a confined cylinder of a LDPE melt, J. Non-Newtonian Fluid Mech. 68 (1997) 173–203. [5] C. Béraudo, A. Fortin, T. Coupez, Y. Demay, B. Vergnes, J.F. Agassant, A finite element method for computing the flow of multi-mode viscoelastic fluids: comparison with experiments, J. Non-Newtonian Fluid Mech. 75 (1998) 1–23. [6] G.B. Bishko, O.G. Harlen, T.C.B. McLeish, T.M. Nicholson, Numerical simulation of the transient flow of branched polymer melts through a planar contraction using the ‘pom-pom’ model, J. Non-Newtonian Fluid Mech. 82 (1999) 255–273. [7] R.J. Blackwell, T.C.B. McLeish, O.G. Harlen, Molecular drag-strain coupling in branched polymer melts, J. Rheol. 44 (1) (2000) 121– 136. [8] A.C.B. Bogaerds, A.M. Grillet, G.W.M. Peters, F.P.T. Baaijens, Stability analysis of polymer shear flows using the eXtended Pom-Pom constitutive equations, J. Non-Newtonian Fluid Mech. 108 (1–3) (2002). [9] A.C.B. Bogaerds, M.A. Hulsen, G.W.M. Peters, F.P.T. Baaijens, Time dependent finite element analysis of the linear stability of viscoelastic flows with interfaces, J. Non-Newtonian Fluid Mech. 116 (1) (2003) 33–54. [10] A.C.B. Bogaerds, W.M.H. Verbeeten, F.P.T. Baaijens, Successes and failures of discontinuous Galerkin methods in viscoelastic fluid analysis, in: B. Cockburn, G.E. Karniadakis, C.-W. Shu (Eds.), Discontinuous Galerkin Methods, Springer, 1999, pp. 264–270. [11] A.C.B. Bogaerds, W.M.H. Verbeeten, G.W.M. Peters, F.P.T. Baaijens, 3D viscoelastic analysis of a polymer solution in a complex flow, Comp. Meth. Appl. Mech. Eng. 180 (3–4) (1999) 413–430. [12] R. Guénette, M. Fortin, A new mixed finite element method for computing viscoelastic flows, J. Non-Newtonian Fluid Mech. 60 (1995) 27–52.

[13] P. Hachmann, Multiaxiale Dehnung von Polymerschmelzen, Ph.D. thesis, Dissertation ETH Zürich No. 11890, 1996. [14] N.J. Inkson, T.C.B. McLeish, O.G. Harlen, D.J. Groves, Predicting low density polyethylene melt rheology in elongational and shear flows with “pom-pom” constitutive equations, J. Rheol. 43 (4) (1999) 873–896. [15] D.G. Kiriakidis, H.J. Park, E. Misoulis, B. Vergnes, J.F. Agassant, A study of stress distribution in contraction flows of an LLDPE melt, J. Non-Newtonian Fluid Mech. 47 (1993) 339–356. [16] P. Lesaint, P.A. Raviart, On a Finite Element Method for Solving the Neutron Transport Equation, Academic Press, New York, 1974. [17] G.H. McKinley, O. Hassager, The Considere condition and rapid stretching of linear and branched polymer melts, J. Rheol. 43 (5) (1999) 1195–1212. [18] T.C.B. McLeish, R.G. Larson, Molecular constitutive equations for a class of branched polymers: the pom-pom polymer, J. Rheol. 42 (1) (1998) 81–110. [19] H. Münstedt, H.M. Laun, Elongational behaviour of a low density polyethylene melt II, Rheol. Acta 18 (1979) 492–504. [20] M. Pahl, W. Gleissle, H.M. Laun, Praktische Rheologie der Kuststoffe und Elastomere, VDI-Gesellschaft Kunststofftechnik, fourth ed., 1995, Chapter 3.8.3, pp. 86–88. [21] G.W.M. Peters, J.F.M. Schoonen, F.P.T. Baaijens, H.E.H. Meijer, On the performance of enhanced constitutive models for polymer melts in a cross-slot flow, J. Non-Newtonian Fluid Mech. 82 (1999) 387– 427. [22] J.F.M. Schoonen, Determination of Rheological Constitutive Equations using Complex Flows, Ph.D. thesis, Eindhoven University of Technology, 1998. [23] J.F.M. Schoonen, F.H.M. Swartjes, G.W.M. Peters, F.P.T. Baaijens, H.E.H. Meijer, A 3D numerical/experimental study on a stagnation flow of a polyisobuthylene solution, J. Non-Newtonian Fluid Mech. 79 (1998) 529–562. [24] H.A. van der Vorst, Bi-CGSTAB: a fast and smoothly converging variant of Bi-CG for the solution of nonsymmetrical linear systems, SIAM J. Scientific Stat. Comput. 13 (1992) 631–644. [25] J. Van Meerveld, Note on thermodynamic consistency of the integral pom-pom model, J. Non-Newtonian Fluid Mech. 1080 (2002) 291– 299. [26] W.M.H. Verbeeten, Computational Polymer Melt Rheology, Ph.D. thesis, Eindhoven University of Technology, 2001. [27] W.M.H. Verbeeten, G.W.M. Peters, F.P.T. Baaijens, Differential constitutive equations for polymer melts: The extended Pom-Pom model, J. Rheol. 45 (4) (2001) 823–843. [28] W.M.H. Verbeeten, G.W.M. Peters, F.P.T. Baaijens, Viscoelastic analysis of complex polymer melt flows using the eXtended Pom-Pom model, J. Non-Newtonian Fluid Mech. 108 (2002) 301– 326. [29] M.H. Wagner, P. Rubio, H. Bastian, The molecular stress function model for polydisperse polymer melts with dissipative convective constraint release, J. Rheol. 45 (6) (2001) 1387–1412. [30] P. Wapperom, R. Keunings, Numerical simulation of branched polymer melts in transient complex flow using pom-pom models, J. Non-Newtonian Fluid Mech. 97 (2001) 267–281. [31] W.F. Zoetelief, J.H.M. Palmen, Experimental results/personal communication, 2000 (DSM Research).