Numerical simulation of flow past a cylinder using models of XPP type

J. Non-Newtonian Fluid Mech. 156 (2009) 7–20

Contents lists available at ScienceDirect

Journal of Non-Newtonian Fluid Mechanics journal homepage: www.elsevier.com/locate/jnnfm

Numerical simulation of flow past a cylinder using models of XPP type N.J. Inkson a , T.N. Phillips a,∗ , R.G.M. van Os b a b

School of Mathematics, Cardiff University, Cardiff CF24 4AG, United Kingdom Etudes et Productions Schlumberger, 1 rue Henri Becquerel, 92140 Clamart, France

a r t i c l e

i n f o

Article history: Received 11 April 2007 Received in revised form 6 June 2008 Accepted 11 June 2008 Keywords: Spectral element method XPP model Flow around a cylinder Multi-mode model

a b s t r a c t We present stable and accurate spectral element methods for predicting the steady-state flow of branched polymer melts past a confined cylinder. The fluid is modelled using a modification of the pom-pom model known as the single eXtended Pom-Pom (XPP) model, where we have included a multi-mode model of a commercial low-density polyethylene. We have analyzed the XPP model and found interesting multiple solutions for certain choices of the parameters which indicate possible problems with the model. The operator-integration-factor-splitting technique is used to discretize the governing equations in time, while the spectral element method is used in space. An iterative solution algorithm that decouples the computation of velocity and pressure from that of stress is used to solve the discrete equations. Appropriate preconditioners are developed for the efficient solution of these problems. Local upwinding factors are used to stabilize the computations. Numerical results are presented demonstrating the performance of the algorithm and the predictions of the model. The influence of the model parameters on the solution is described and, in particular, the dependence of the drag on the cylinder as function of the Weissenberg number. © 2008 Elsevier B.V. All rights reserved.

1. Introduction The mathematical modelling of polymer melts has improved over the last decade due to progress in understanding of the molecular origin of the polymer physics for branched polymers. This new theory leads to the development of the pom-pom constitutive model for branched polymers by McLeish and Larson [14]. This model has been successful in modelling the nonlinear behaviour in both shear and planar extension simultaneously in contrast to traditional models such as the Giesekus [7,8] and K-BKZ [1] models. The success of the pom-pom model is due to the incorporation of molecular dynamics into the modelling process based on considerations of the polymer physics. While the Giesekus model can give qualitatively both shear-thinning and extension-hardening behaviour, it lacks the flexibility to quantitatively fit polymer melt data. The K-BKZ model can fit uniaxial extensional and shear viscosities but will simultaneously fail to predict the correct planar extensional viscosity. Modifications to the original pom-pom model [2,26] have allowed quantitative agreement between experiments and model predictions obtained using numerical simulation techniques. The inclusion of multiple modes [4,10,27,28] allows the fitting of

∗ Corresponding author. Tel.: +44 2920 874194; fax: +44 2920 874199. E-mail addresses: Nathanael [email protected] (N.J. Inkson), [email protected] (T.N. Phillips), [email protected] (R.G.M. van Os). 0377-0257/$ – see front matter © 2008 Elsevier B.V. All rights reserved. doi:10.1016/j.jnnfm.2008.06.004

low-density polyethylene (LDPE) rheology data. The rheology of such a complex branched material spans frequencies over several orders of magnitude. This corresponds to the wide separation of relaxation times associated with the randomly branched topology of the molecules. In this way different parts of the molecule can be assigned different relaxation times, which relax on different time-scales. This is modelled mathematically using multiple Maxwell modes each with a characteristic relaxation time. The importance of reptation in the mathematical description of linear polymers was identified by de Gennes [5] and later quantified by Doi and Edwards [6]. This concept assumes that individual polymer molecules are constrained to move in a tube formed by neighbouring polymer molecules. Effectively, the perpendicular motion of a polymer chain is constrained within a tube of given radius generated along its backbone, while parallel to the tube a linear polymer chain is free to diffuse. The key features in the mathematical modelling of branched polymers are the topology of the polymer molecules and the related relaxation processes which enable the molecule to escape its tube. For branched polymers more complicated relaxation processes are involved than for linear polymers. For star polymers it was shown that reptation does not occur. Rather, a retraction of the star polymer arm into the tube and then a withdrawal occurs which allows the polymer to relax its stress [15]. The pom-pom constitutive model of McLeish and Larson [14] identified that H-shaped molecules and pom-poms would relax by a combination of the star fluctuation and reptation theory once the

8

N.J. Inkson et al. / J. Non-Newtonian Fluid Mech. 156 (2009) 7–20

arms had relaxed. This came from the realisation that polymer chain segments that are trapped between branch points, contribute the most to the extra-stress in a strongly deforming flow. The model assumes that the polymer chains are represented by a backbone segment connecting two identical pom-poms each with q arms at the branch points. The drag that the melt exerts on these arms, as a result of friction, causes the backbone to stretch. The branched arms inhibit the reptation motion of the backbone section by pinning the molecule in place at the tube junctions. The free ends of the arms are still able to move however, and the arms gradually work their way out of the tubes, in a similar fashion to star polymers, towards the branch points by a diffusion process called arm retraction [15]. The arm retraction process is accelerated by the self diluting action of the arm ends. The arms are able to explore configuration space rapidly and act to remove the constraints for neighbouring chains which results in dilated tubes, in a process known as dynamic dilution. Once the arms have relaxed, the backbone can subsequently relax by moving the branch points and reptating. The backbone can then be treated as a linear polymer with two diffuse ‘blobs’ comprising the diffusing arms, and the chain can be modelled as a dumbbell inside a tube. A key feature of the model is therefore distinct relaxation times for these two processes, i.e. the orientation and stretch of the backbone section. This separation of relaxation times gives the flexibility of the pom-pom model to quantitatively describe polymer melt rheology. Due to the balance of tension in the backbone of the molecule, and in the arms, the stretch is not allowed to be greater than q. At this point, the arms are retracted into the backbone tube and this ensures that the stretch is bounded, although retraction lengths are found to be very small and can be neglected [10]. Some problems have been identified with the original differential form of the pom-pom model, including: the prediction of a zero second normal-stress difference and an unphysical discontinuity in the gradient of the extensional viscosity when the stretch is at a maximum (equal to q) at steady state. A modification of the stretch equation by Blackwell et al. [2] allows for the branch point to move and withdraw into the tube. This has the effect of smoothing the transition to the maximum extensional viscosity plateau and removes the discontinuity in the gradient of the extensional viscosity. Verbeeten et al. [26] introduced a modified model, the eXtended Pom-Pom (XPP) model, to circumvent these problems. The constraint that the stretch cannot be larger than the number of arms is removed. A Giesekus-like term has been added to the orientation equation. This term adds a non-zero second normal-stress difference in order to describe the behaviour of polymer melts more completely. The orientation equation is bounded for high strain rates (it was not in the original model). The XPP model has been used along with a finite element method to model LDPE melts (see Verbeeten et al. [27]). We performed a numerical study to investigate the use of a high-order spectral element method in predicting the flow of an LDPE past a confined cylinder in a planar channel. A multi-mode XPP model is used to realistically describe LDPE rheology. A single mode analysis of the flow around a cylinder was presented in [24]. This paper studied the velocity profile and stresses in planar channel Poiseuille flow at different Weissenberg numbers, and the effect of changing the ratio of the stretch relaxation time to the orientational relaxation time in the extended pom-pom model. The flow around a cylinder in a planar channel was also considered and the variation on the drag on the cylinder was found to decrease with increasing Weissenberg number. Material parameters were again varied to see the effect of pom-pom arm number, Newtonian solvent effects and Reynolds number. The effect of changing individual model parameters has been established for flow around a cylinder.

In this paper we consider a multi-mode version of the model with parameters selected corresponding to a commercial low-density polyethylene. Instead of individually varying the Weissenberg and Reynolds numbers we vary the characteristic velocity which we define as the average inflow velocity. In this way both non-dimensional numbers change for each simulation. The chosen LDPE was originally studied by Verbeeten et al. [28]. In their study, the finite element method was used to discretize the weak formulation of the problem. Enhanced stability was achieved using the DEVSS technique in combination with the discontinuous Galerkin method (DEVSS/DG). As in Van Os and Phillips [24] we use a numerical scheme based on a decoupled approach in which the field equations are treated separately from the constitutive equation within each time step. This scheme provides a considerable improvement in terms of stability over a previous scheme proposed by the authors [22]. The improvements are attributed to the particular decoupled scheme that is implemented and the use of a suitable upwinding strategy. The combined effect of these improvements was to increase the maximum attainable Weissenberg number by more than an order of magnitude. In the coupled scheme of Van Os and Phillips [22] most of the terms in the constitutive equation were treated explicitly in order to obtain a symmetric positive definite system for velocity, pressure and extra-stress, which could be solved efficiently using the conjugate gradient method. In this paper the velocity and pressure are decoupled from the computation of the extra-stress. This results in a symmetric positive definite system for velocity and pressure which is retained through the standard spectral element discretization of the mass and momentum equations. The use of a decoupled scheme allows treatment of the constitutive equation in a more implicit fashion. Preconditioned iterative methods for are used for solving the discrete equations at each time step. A Schur complement method reduces the size of the Helmholtz operator. This method eliminates the interior degrees of freedom associated with each spectral element to leave a system that serves to determine the unknowns at element interfaces. The Schur complement is stored as an LU decomposition and used as preconditioner for the conjugate gradient method. The most efficient preconditioner for the Uzawa operator was found to be based on an element-wise LU decomposition of the spectral element pressure matrix. The non-symmetric operator arising from the spectral element discretization of the constitutive equation is solved using the Bi-CGStab iterative method. The stress mass matrix is a simple yet efficient preconditioner for this operator. Upwinding techniques were implemented in order to improve further the stability of the spectral element approximation for high values of the Weissenberg number. The choice of upwinding method was the ‘locally upwinded spectral technique’ (LUST) developed by Owens et al. [17], which employs local upwinding factors. This was found to enhance stability. This paper is organised as follows: Section 2 introduces the governing equations for the XPP model, and describes the nondimensionalization of the equations used. Section 3 describes some of the important rheological properties of the original and extended pom-pom models. The temporal and spatial discretization of the governing equations is described in Section 4. The decoupled scheme used to solve the discrete equations and the locally upwinded spectral technique is also given in this section. Details of the solution algorithm and the performance of the preconditioners are given in Section 5. In Section 6 numerical results are presented for the flow past a cylinder. The conclusions drawn from this study are presented in Section 7.

N.J. Inkson et al. / J. Non-Newtonian Fluid Mech. 156 (2009) 7–20

2. Governing equations for the XPP fluid The flow of polymer melts is governed by equations consisting of two field equations and a constitutive equation. The constitutive equation used in this paper is given by the single equation version of the extended pom-pom model in multi-mode form. In the single equation version of the model the stress is directly related to the strain, as with all classical closed form continuum constitutive models. In the double equation version of the model, an equation of state for the orientation of the polymer molecule is provided together with an evolution equation for the stretch of the molecule, akin to the original pom-pom model. The flow of an isothermal, incompressible fluid is described by the equations for conservation of mass and momentum:

∇ · u = 0,





∂u + u · ∇u ∂t

(1)

 = −∇ p + ∇ · T,

(2)

where u is the velocity field,  is the density of the fluid, p is the pressure, and T is the extra-stress tensor. The extra-stress tensor T is decomposed into polymeric ( p ) and viscous (2s d) contributions T =  p + 2s d,

(3)

where s is the solvent viscosity and d is the rate of deformation tensor, given by d=

1 (∇ u + (∇ u)T ). 2

(4)

Furthermore,  p itself may be decomposed into a number of discrete modes according to p =

n 

i,

(5)

i=1

parameter i in the exponential term was incorporated into the stretch relaxation time to remove the discontinuity from the gradient of the extensional viscosity. Its value was found to be inversely proportional to the number of arms qi , i =

T=

(10)

 i = Gi (32i si − I),

 i + 2s d.

(6)

It remains to define a constitutive equation, which will serve to determine the stress  i , i = 1, . . . , n, associated with each mode. The evolution equations for the orientation tensor si , and the backbone stretch i , for the double extended pom-pom model of Verbeeten et al. [26] are given by si + fi (i , si )si +

˛i − 1 32i b,i

I+

3˛i 2i b,i

si · si = 0,

(7)

and ( − 1)  ( −1) ∂i + u · ∇ i = i (d : si ) − i e i i , s,i ∂t

∇

fi (i ,  i ) i + b,i  i + Gi (fi (i ,  i ) − 1)I +

1 2i b,i

(1 − ˛i − 3˛i 4i Isi ·si ).

˛i  ·  = 2b,i Gi d, Gi i i

in which the function fi (i ,  i ) is given by fi (i ,  i ) = 2



b,i  ( −1) 1 e i i 1− s,i i



+

1



2i

1−

˛i I i · i

(12)



3Gi2

.

(13)

In the single equation version of the XPP model, the backbone stretch i is expressed explicitly in terms of the extra-stress by taking the trace of (11). This yields the relation i =

1+

I i 3Gi

.

(14)

The constitutive equations for the XPP model are made dimensionless. The non-dimensional variables are defined as follows: x∗ =

x , L

u∗ =

u , U

t∗ =

Ut , L

p∗ =

Lp , U

∗ =

L , U

where L is a characteristic length scale,  is the total viscosity and U is a characteristic velocity. As we are considering a multi-mode description of the fluid, a Weissenberg number, Wei , is associated with each mode through the corresponding relaxation times Wei =

(8)

s,i U . L

(15)

The Reynolds number, Re and the viscosity ratio parameters ˇ and ˜ i are defined by ˇ

respectively, where the function fi (i , si ) is given by fi (i , si ) = 2d : si +

(11)

where Gi is the shear modulus from a linear relaxation spectrum. Note that the orientation tensor exhibits a non-zero third normalstress component, even in 2D problems. The equation for this component may be eliminated from the system of equations since it is known that the trace of the orientation tensor, Is , is unity. Eqs. (7)–(9) form the double equation version of the extended pom-pom model. Eq. (7) may be reformulated in terms of the extrastress tensor associated with the ith mode rather than in terms of the orientation tensor by multiplying (7) by 2i . This results in the equivalent single equation model, which has the same structure as the Giesekus or Phan-Thien–Tanner models [18,19]. It relates the polymeric stress to the rate of deformation tensor through



i=1

∇

2 . qi

The extra-stress tensor can be evaluated from the orientation tensor using the following relation

so that n 

9

(9)

Here b,i and s,i are the relaxation times of the backbone tube orientation and the backbone stretch for each mode, respectively. The incorporation of an additional nonlinear term into the evolution equation for the orientation tensor provides the model with the features of a modified Giesekus equation, which introduces a similar ‘anisotropy parameter’, ˛i , for each mode. When ˛i = / 0a non-zero second normal-stress difference is predicted, therefore giving a more complete mathematical description of the polymer melt. The parameter i provides a measure of the branch point withdrawal process which determines the backbone tube stretch. The

Re =

UL , 

ˇ=

s , s + p

˜i = ˇ

p,i p

,

(16)

where the polymeric viscosity p is defined in terms of moduli Gi , from the linear relaxation spectrum. This modulus can be defined as a viscosity divided by a timescale. This timescale is chosen to be b,i , so the total polymeric viscosity, p , is given by p =

n 

Gi b,i

and

p,i = Gi b,i .

(17)

i=0

A further non-dimensional material parameter is needed for the pom-pom model. This is the parameter i , which is defined

10

N.J. Inkson et al. / J. Non-Newtonian Fluid Mech. 156 (2009) 7–20

as the ratio of the two relaxation times processes of the pom-pom molecule which are orientation (reptation) and stretch relaxation. εi =

s,i . b,i

(18)



Re

(19)



∂u + u · ∇u ∂t

= −∇ p + ∇ ·  + ˇ∇ 2 u,

(20)

∇

fi (i ,  i ) i + Wei  i +

˜ i (1 − ˇ) ˇ [fi (i ,  i ) − 1] I Wei

˛i Wei ˜ i (1 − ˇ)d,  i ·  i = 2ˇ ˜ ˇi (1 − ˇ)

(21)

where fi (i ,  i ), and i are given by fi (i ,  i ) =

2 i +

and i =



1−

1



1 i



1−

2i

ei (i −1)

2

Wei ˜ i (1 − ˇ) ˇ

 ˛i I · 3 i i

,

(22)

1+

Wei 1 I . ˜ i (1 − ˇ) 3 i ˇ

(23)

2 fi (i ,  i ) = i +



1−

1 2i

1 2i





1−

ei (i −1)

2

Wei ˜ i (1 − ˇ) ˇ

 ˛i I · 3 i i

.

(24)

This equation differs from (22) in that the term 1/i in the first bracket is replaced by 1/2i . This is the functional form for fi we shall use hereafter. Finally, we compare the XPP model to the double pom-pom (DPP) model of Clemeur et al. [3]. This model also has separate equations for the stretch and molecular orientation, and was derived from the double XPP model by removing the Giesekus-like term, which is thought to be problematic in some situations [3], and replacing d by  in the evolution equation for backbone stretch (8). This results in the DPP model: ∇

Si + 2i b S i + 2 2i b ( : Si )Si −

I 3

= 0,

D i ∗ − s (i : Si ) i + ( i − 1) e ( i −1) = 0, Dt

(25) (26)

where  = (∇ u)T . The extra-stress is computed using the relation  i = Gi (3 2i Si − I).

(27)

This model has a zero second normal-stress difference like the original pom-pom model. 3. Rheological properties In this section we present some of the rheological properties of the XPP model for a single mode. We begin by investigating the prediction of the model for steady simple shear flow for which u =

˛ 2 2 ( + xy )=0 G0 xx

˛ 2 2 ( + xy )=0 G0 yy

˛ (xx + yy )xy − b G0 ˙ = 0 G0 ˛ 2  =0 G0 zz

(28) (29) (30) (31)

where f (, ) and  are given by Eqs. (13) and (14), respectively, with the modal subscript omitted. Since f (, ) and  are nonlinear functions of the components of the extra-stress tensor and the constitutive equation contains quadratic terms involving these components, Eqs. (28)–(31) collectively represent a nonlinear system of algebraic equations for the extra-stress components for each given value of . ˙ In the case of uniaxial extensional flow for which u = ˙ −y/2, ˙ ˙ (x, −z/2) where ˙ is the extensional rate the constitutive Eq. (12) takes the following component form ˙ xx +G0 )+G0 (f (, ) − 1)+ f (, )xx −2b ( ˙ yy +G0 (f (, ) − 1)+ f (, )yy +b 

Recently, it has been argued [21] that, due to thermodynamical arguments, a more appropriate functional form for fi is

s

f (, )xy − b yy ˙ +

f (, )zz + G0 (f (, ) − 1) +

while the non-dimensional XPP model is given by

+

˙ b + G0 (f (, ) − 1) + f (, )xx − 2xy  f (, )yy + G0 (f (, ) − 1) +

The non-dimensional conservation equations are

∇ · u = 0,

( y, ˙ 0, 0). In this case the constitutive Eq. (12) takes the following component form

˛ 2 ( + 2 ) = 0 G0 xx xy

˛ 2 ( + 2 )+b G0 ˙ = 0 G0 yy xy

(32) (33)

with xy = 0 and yy = zz . The nonlinear systems of Eqs. (28)–(31) and (32)–(33) are ˙ respectively, using the solved for a given range of ˙ and , Newton–Raphson method in which the Jacobian is constructed analytically. For illustrative purposes we fixed the following parameter values: G0 = 1, b = 3, s = 1 and varied the values of q and ˛. Evidence of multiple solutions to the above systems of equations was found when q was sufficiently large and for all values of ˛ studied. This problem has been previously noted in the literature by Clemeur et al. [3], but has not been studied extensively. In an early paper by Verbeeten et al. [27], it was observed that, when fitting viscosity and first normal-stress difference data for LDPE, a suitable rule of thumb for choosing the value of ˛ was ˛ ≈ 0.3/q. However, in extensional flow, evidence of multiple solutions was found even when the value of ˛ was small. The occurrence of multiple solutions seems to be more apparent in uniaxial extension. There are clearly some numerical issues to be addressed concerning the use of the XPP model. It is interesting that numerical integration of the transient problem does not seem to present a problem, in that the correct value of the viscosity is obtained in all the cases we studied. Multiple solutions of the steady-state equations were found by choosing 10,000 starting values for the Newton iterations for each value of . ˙ Some of the starting values were chosen to be the solutions corresponding to the previous shear rate, while others were generated randomly. We choose values of q and ˛ that correspond to the two modes with the longest relaxation times of the LDPE studied with the XPP model in Verbeeten et al. [27]. Fig. 1(a) shows the numerical solutions for the steady shear viscosity for a discrete set of shear rates for the XPP model for q = 10, ˛ = 0.03. At the lowest shear rates only one solution is present at the correct plateau. At higher shear rates a second solution appears at a higher plateau, and at even higher rates a third solution is found. Note that if we choose a starting value of zero for the stress we will find the correct lower solution [11]. For steady-state uniaxial viscosity (see Fig. 1(b)) we find that there is a second solution, below the correct curve describing the material rheology.

N.J. Inkson et al. / J. Non-Newtonian Fluid Mech. 156 (2009) 7–20

11

Fig. 1. Predictions of the XPP model: (a) shear viscosity and (b) uniaxial extensional viscosity for q = 10, ˛ = 0.03.

Fig. 2. Predictions of the XPP model: (a) shear viscosity and (b) uniaxial extensional viscosity for q = 2, ˛ = 0.15.

For values of XPP parameters corresponding to less branching, q = 2 and ˛ = 0.15, we find that again multiple solutions are clearly seen in shear (see Fig. 2(a)), and in extension (see Fig. 2(b)). In Fig. 2(b) there appears to be two plateau values, the correct viscosity in this case is the higher one with the plateau at, E,0 = 9.0 Pa s. The lower branch of solutions then crosses the upper curve. We also found multiple solutions with the modification of f (, ) equation due to the recommendation of van Meerfeld [21](Eq. 24). We found many such cases of multiple solutions of the steady-state viscosity equations, these results are presented in a separate paper [11]. We found that if the initial value of stress is taken to be zero and if stress is computed in a transient fashion then the correct solution is usually found for the XPP model. We therefore assume in the following sections that the correct stresses are calculated for the flow past a cylinder, as we start from zero stress. For the flow simulations, we chose to model the well characterized LDPE DSM Stamylan 2008 XC43 [27]. We altered the pom-pom parameters for the characterization of the material slightly as the Van Meerfeld correction acts to reduce the maximum extensional viscosity. In order to compensate for this the value of q associated with the largest mode was increased from 10.0 to 10.8, and the ratio of the relaxation times (1/i ) was decreased from 1.1 to 1.05. The original parameter set was from the original fit of this LDPE [27]. The parameters are summarized in Table 1. We compare the XPP

model to the DPP model of Clemeur et al. [3], the parameters of which are provided in Table 2. Since we are calculating the transient uniaxial extensional viscosity, we do not expect to run into the problem of multiple solutions since we are approaching the numerical solution from below. The extensional viscosity was calculated using a fourth order Runge–Kutta numerical integration using a 4-mode XPP fluid for strain rates of 0.38, 1.12, 3.82 s−1 and compared to experimental data in Fig. 3(a) (data scanned from [27]). Note that the extensional viscosity for the strain rate of 3.82 s−1 does not reach a plateau in contrast to the experimental data. In fact the uniaxial extensional viscosity was unbounded for extensional rates above 2 s−1 , which was found by trial and error. We are unsure as to the cause of this behaviour. This could be due to another solution being reached, or it could be a numerical artifact of the integration technique due to the differential equations being stiff. We discovered that the blow-up in stress was due to the Giesekus-like term in the stress equation. However, in the simulations we reach strain rates higher than this critical strain rate before divergence in the solution occurs. A recently published variation of the double extended pom-pom model, the ‘double pom-pom model’ [3], removes the unbounded stress by setting ˛i = 0, where i = 1, . . . .,n. A plot of predicted extensional viscosity is shown in Fig. 3(b).

Table 1 Parameters for the XPP model used to describe the LDPE 2008 XC43

Table 2 Parameters for the DPP model used to describe the LDPE 2008 XC43

b,i (s)

Gi (Pa)

qi

1/i

˛i

b,i (s)

Gi (Pa)

qi

1/i

0.0038946 0.05139 0.50349 4.5911

72006 15770 3334 300.8

1 1 2 10.8

7 5 3 1.05

0.3 0.3 0.15 0.03

0.0038946 0.05139 0.50349 4.5911

72006 15770 3334 300.8

1 1 4.4 10.0

2 2 2 1.1

12

N.J. Inkson et al. / J. Non-Newtonian Fluid Mech. 156 (2009) 7–20

Fig. 3. Comparison of the transient uniaxial extensional viscosity for the LDPE 2008 XC43 obtained from experiments and model predictions (a) XPP and (b) DPP, for strain rates 0.38, 1.12, 3.82 s−1 .

For the multi-mode spectral element simulations the first mode is taken to be a Newtonian solvent. Since this mode has q = 1 it does not contribute to the extra-stress in extension. The inclusion of a solvent contribution aids the calculation and allows a larger timestep to be used for the numerical simulations. The value of this viscosity is s = G1 b,1 = 0.06756 Pa s. 4. Temporal and spatial discretization

have been calculated, n+1 respectively. Once pn+1 , un+1 and  n+1 i i is found from Eq. (23). 4.2. Spectral element discretization The spectral element method is applied to the weak formulation of the semi-discrete Eqs. (34), (35) and (37). Suitable function spaces are chosen for the dependent variables. The velocity is cho2

4.1. Temporal discretization For the temporal discretization of the equations, a scheme that decouples the computation of velocity and pressure using the field equations from the computation of the extra-stress using the constitutive equation was used. This scheme is described in detail in the paper of Van Os and Phillips [24]. The discretization of the material derivative in the momentum equations is performed using a first-order operator-integrationfactor-splitting (OIFS) scheme (see Maday et al. [13]). This yields the following semi-discrete approximation of the field equations

∇ · un+1 = 0,

(34)

Re n+1 ˜ n+1 )) = −∇ pn+1 + ∇ ·  n + ˇ∇ 2 un+1 , (u − u(t t

(35)

˜ n+1 ) is the solution at time t = t n+1 of where the approximation u(t the pure convection problem ˜ ∂u ˜ = −un · ∇ u, ∂t

t ∈ [t n , t n+1 ],

˜ with u(x, t n ) = un (x).

(36)

A fourth order explicit Runge–Kutta (RK4) method is used to solve the initial value problem (36) using a time step s satisfying t = 8 s. In the constitutive equation, the upper-convected derivative of the extra-stress together with the rate of deformation tensor are discretized implicitly: C(ni ,  ni ) n+1 + Wei (u · ∇  i − (∇ u) i −  i (∇ u)T ) i +

4.3. The discrete equations

n+1

Applying the spectral approximations to the weak formulation of the conservation Eqs. (34) and (35) yields the following set of discrete equations

˛i Wei ˜ i (1 − ˇ)dn+1  i n ·  i n+1 = 2ˇ ˜ ˇi (1 − ˇ)

+ hi (ni ,  i n , un ),

(37)

where C(i ,  i ) and hi (i ,  i , u) are defined by Wei C(i ,  i ) = fi (i ,  i ) + t

(38)

and hi (i ,  i , u) =

˜ (1 − ˇ) Wei ˇ  − (fi (i ,  i ) − 1) i I, Wei t i

sen to be in a subspace, V, of [H 1 (˝)] whose elements satisfy the prescribed velocity boundary conditions. The weak formulation of the generalized Stokes problem (34) and (35) is well-posed if the pressure space is chosen to be P = L02 (˝). In the spectral element method the weak formulation of the semi-discrete Eqs. (34), (35) and (37) is discretized using highorder polynomial approximations in space. The physical domain ˝ is divided into K non-overlapping spectral elements ˝k , 1 ≤ k ≤ K, such that ∪Kk=1 ˝k = ˝. We denote by PN (˝k ) the space of all polynomials on ˝k of degree less than or equal to N, and further define PN (˝) = { : |˝k ∈ PN (˝k )}. Each of the spectral elements is mapped onto a parent element D = [−1, 1] × [−1, 1], where each point (, ) ∈ D is associated with a point (x(, ), y(, )) ∈ ˝k . The dependent variables are approximated on D using Lagrangian interpolants of degree N in both spatial directions, based on the Gauss–Lobatto–Legendre (GLL) points. This creates a Gauss–Lobatto–Legendre grid within each spectral element. The discrete approximation spaces must satisfy a compatibility condition to ensure that the problem is well-posed. For spectral elements, Maday and Patera [12] have shown that this LBB condition is satisfied when the velocity approximation space is the polynomial space PN (˝), and the pressure approximation space is PN−2 (˝). The extra-stress is also approximated by polynomials in the space PN (˝) as well. A continuous approximation is used for the extra-stress.

(39)

=0 DN un+1 N

(40)

T n+1 − DN pN = gN − BN  nN (Ml CN + ˇEN )un+1 N

(41)

where DN and BN are the discrete divergence operators acting on T and BT are gradient operators velocity and stress, respectively, DN N acting on pressure and velocity, EN is the discrete Laplace operator, CN is the discrete velocity mass matrix, and gN includes the boundary conditions and the remaining source terms.

N.J. Inkson et al. / J. Non-Newtonian Fluid Mech. 156 (2009) 7–20

The discrete Helmholtz operator, HN , is defined by HN = Ml CN + ˇEN .

of the quadratic equation (42)

Eliminating the velocity from (41) using the discrete continuity Eq. (40) yields the pressure equation



T n+1 pN = −DN HN−1 gN − BN  nN , DN HN−1 DN

(43)

T on the left-hand side of (43) is The operator UN = DN HN−1 DN known as the Uzawa operator. The pressure equation can be expressed in terms of this operator, i.e.

UN pn+1 = bN , N

(44)

where bN is the right-hand side of (43). Once the pressure has been found by solving (44), the velocity is determined from (41). The discrete weak formulation of the constitutive Eq. (37) is ˜ i (1 − ˇ)BT un+1 = hN , SN  n+1 −ˇ N N N

(45)

in which SN is the non-symmetric discrete operator that contains all the terms in the left-hand side of the constitutive Eq. (37). The righthand side hN is the discrete form of (39) in the weak formulation. 4.4. The locally upwinded spectral technique Upwinding techniques such as SUPG and the locally upwinded spectral technique developed by Owens et al. [17] are used to improve the stability of the spectral element calculations. In both of these techniques a modified test space is used for the constitutive equation. In particular, if  is a test function in the original weak formulation of the constitutive equation then  + hu · ∇  is the modified test function within an SUPG or LUST formulation where h is the scalar upwinding factor. The difference between SUPG and LUST lies in the choice of the upwinding factor h. In SUPG this is a global quantity and for spectral element methods it is common to use h = 1/N 2 , where N is the order of the spectral approximation. The upwinding factors in LUST are chosen to achieve the property of superconsistency in the quadrature rules used to approximate the weak formulation of the constitutive equation. More precisely, each of the upwinding factors are determined using the quadrilateral mini-elements formed by four GLL nodes inside a spectral element. Let xij be one of the four corners in a mini-element for which uN is directed towards the exterior of the mini-element. Define the polynomial PN+1 (x) = (1 − 

2

  )LN ()(1 − 2 )LN (),

(46)

where LN is the Legendre polynomial of degree N. If we assume that PN+1 ≥ 0 in this mini-element then it is easy to see that PN+1 − Wei (u · ∇ )PN+1 ≤ 0. / 0, one can show that there is a point x∗ij along Now, for Wei = the streamline passing through xij where PN+1 attains a maximum and, at this point, we have PN+1 − Wei (u · ∇ )Pn+1 ≥ 0. By the intermediate value theorem there must exist a point zjk = xjk − hjk u(xij ) for which PN+1 (zjk ) + Wei (u · ∇ )PN+1 (zjk ) = 0.

13

(47)

By expanding PN+1 (zij ) as a Taylor series about xij and neglecting terms of third order and above, Owens et al. [17] showed that, for internal nodes, the shift or upwinding factor, hij , is the positive root

2N(N + 1) 3



u2ij 1 − 2

+

vij2 1 − 2

 h2ij + hij − 2 = 0.

(48)

Similarly, the shift factors corresponding to the edges of the elements can also be found. Note that for the multi-mode case we have different LUST shift factors for each mode as the Weissenberg number varies. 4.5. Iterative solution of the discrete equations Since the computation of velocity and pressure is decoupled from the computation of the extra-stress, a symmetric positive definite system for velocity and pressure is retained through the standard spectral element discretization of the mass and momentum equations. A Schur complement method has been used to reduce the size of the Helmholtz operator. This method eliminates the interior degrees of freedom associated with each spectral element to leave a system that serves to determine the unknowns at element interfaces. The Schur complement is stored as an LU decomposition and used as preconditioner for the conjugate gradient method. The preconditioner that was found to be the most efficient for the Uzawa operator is based on the finite element mesh comprised of the inner Gauss–Lobatto–Legendre nodes of each spectral element. The mesh is triangulated with no overlapping nodes. This preconditioner has been found to work very well in conjunction with the uncoupled method employed in this paper. This preconditioner is based on the finite element mass and stiffness matrices on the local finite elements, which are MkFE and EkFE and it is given as: PU−1 = FE−1 U =

 FE K  Mk T Rk

k=1

ˇ

−1

+ Ml EkFE

Rk

(49)

where Rk is the restriction operator which maps a global vector to a vector of length equal to the number of GLL nodes of the spectral element. The preconditioner is stored as an LU decomposition. The non-symmetric operator arising from the spectral element discretization of the constitutive equation is solved using the BiCGStab iterative method. The stress mass matrix is a simple but efficient preconditioner for this operator [24]. 5. Flow past a cylinder Consider the flow past a cylinder placed symmetrically in a channel of 2 units half-width and whose centre is located at the origin. The ratio of the radius of the cylinder to the half-width of the channel is 0.5. The geometry of the problem is shown in Fig. 4. The computational domain extends a distance 25 units upstream and 25 units downstream of the cylinder so that the assumption of fully developed flow conditions at entry and exit is valid for the range of Weissenberg numbers considered. A comprehensive discussion of this benchmark problem of flow around a cylinder may be found in the monograph of Owens and Phillips [20]. The effect of varying the individual parameters of the XPP model for flow around a cylinder has been discussed in previous work [24,25]. This paper is primarily concerned with the modelling of a commercial grade material in order to compare the numerical predictions with published experimental data. In this paper, the Weissenberg number is increased by increasing the flow rate rather than by changing the relaxation time since this effectively alters the fluid. An increase in flow rate is achieved in the numerical simulations by increasing the value of U. For the complex

14

N.J. Inkson et al. / J. Non-Newtonian Fluid Mech. 156 (2009) 7–20

Fig. 4. Spectral element mesh for flow past a solid cylinder for −5 ≤ x ≤ 5, K = 20 and N = 6.

flow of a polymer melt around a cylinder we impose a fully developed velocity profile at the inflow boundary. Unlike models such as the UCM and Oldroyd B [16] models, the XPP model does not possess an analytical solution for plane Poiseuille flow. To calculate the necessary inflow conditions, we simulated flow in a channel with periodic boundary conditions, for a range of characteristic velocities. The dimensionless parameters Re and We are defined with reference to U. This corresponds to a change in the average inflow velocity. The channel was of non-dimensional length of 64 to avoid any destabilizing effect of a short channel [23]. This provided the steady-state velocity and stress inflow conditions for predicting the flow around a cylinder. The mesh consisted of two spectral elements in the cross channel direction, with a line of symmetry in the centre. The polynomial order in the spectral approximation was varied to ensure that convergence with increasing polynomial order was obtained. The flow past a cylinder is calculated for several characteristic velocities, U. The characteristic length is equal to 1 (m). For this LDPE  = 0.92 g cm−3 giving a Reynolds number, Re = 0.2217U. Mesh convergence was investigated by changing N and K. For the characteristic velocity U = 0.1 ms−1 , the resulting dimensionless drag force, F, on the cylinder are shown in Table 3. The dimensionless stress component xx is plotted along the centreline of the channel and around the cylinder surface in Fig. 5(a). Fig. 5(b)

Table 3 Dependence of the dimensionless drag force, F and F˜ on N, U =0.1 ms−1 N

F F˜

4

5

6

7

116.78 0.883

118.84 0.898

118.51 0.896

118.64 0.897

shows the centreline stress component xx for N = 4, 6 and 7 with K =20. Zoomed areas of interest are shown in Fig. 6(a) and (b) for ˜ U = 0.1 and 0.5 ms−1 respectively. The normalised drag force, F, was obtained by dividing F by the drag corresponding to the characteristic velocity of U = 10−6 ms−1 . This latter value was 132.32. The results presented in the rest of this paper were obtained with N = 6 since this presents a sufficient level of refinement without excessive computational cost. Table 4 shows the drag force on the cylinder for each characteristic velocity along with the Weissenberg number corresponding to the slowest relaxing mode, Wemax . The highest characteristic velocity that yielded a converged numerical solution was U = 1.64 ms−1 . Simulations failed to converge for U > 1.64 ms−1 . The value U = 1.64 ms−1 corresponds to a Weissenberg number of 7.529 for the slowest relaxing mode. In order to achieve the necessary level of convergence for the highest characteristic velocity, the size of the time step was reduced from

Fig. 5. Mesh convergence of the dimensionless stress component xx for (a) U = 0.1 ms−1 and (b) U = 0.5 ms−1 .

Fig. 6. Zoom of the region near x = 1.5 of Fig. 5 showing the detail of the dimensionless stress component xx for (a) U = 0.1 ms−1 and (b) U = 0.5 ms−1 .

N.J. Inkson et al. / J. Non-Newtonian Fluid Mech. 156 (2009) 7–20

15

Table 4 ˜ on the cylinder, as a function of the characteristic velocity, U, together with the Weissenberg number corresponding Dimensionless drag force, F, and normalised drag force, F, to the slowest relaxing mode U

Wemax F F˜

0.005 ms−1

0.01 ms−1

0.1 ms−1

0.5 ms−1

0.75 ms−1

1.00 ms−1

1.64 ms−1

0.023 132.30 0.999

0.046 131.21 0.992

0.459 118.51 0.896

2.296 92.83 0.702

3.443 75.80 0.573

4.591 68.50 0.518

7.529 68.71 0.519

Fig. 7. Contours of the shear stress for U = 0.1 ms−1 .

Fig. 8. Contours of the shear stress for U = 0.5 ms−1 .

Fig. 9. Contours of the shear stress for U = 0.75 ms−1 .

Fig. 10. Contours of the shear stress for U = 1.64 ms−1 .

1 × 10−4 to 1 × 10−6 . The results reflect the nature of LDPE melt rheology. In the region of extensional flow between the cylinder and the wall, some strain hardening is predicted with an increase in normal-stress difference for lower values of U, which correspond to lower strain rates. At higher strain rates, beyond the peak in steady extensional viscosity, we expect more strain-softening. The shear viscosity response reflects the shear-thinning properties of the melt. This results in a drag reduction as U and therefore Wemax increases. We now analyze in more detail the contour plots of the shear stress, first normal-stress difference N1 and the stretch of the slowest relaxing mode,  for the characteristic velocities studied, U = 0.1, 0.5, 0.75 and 1.64 ms−1 . Contour plots for shear stress are shown in Figs. 7–10 respectively, for U = 0.1, 0.5, 0.75, 1.64 ms−1 . The minimum value of the shear stress occurs either on the channel wall opposite the cylinder or along the downstream bound-

ary of the cylinder. The maximum value of the shear stress occurs along the upstream boundary of the cylinder (see Fig. 7). The extreme values of the shear stress and their co-ordinates are given in Table 5. For U > 0.5 ms−1 , we find that the region of large positive shear stress moves towards the front of the cylinder while a

Table 5 Maximum and minimum values of the dimensionless shear stress on the cylinder and the minimum value on the channel wall as a function of U. U (ms−1 )

0.1 0.5 0.75 1.64

Cylinder

Wall

max

(x, y)

min

(x, y)

min

6.04 5.91 4.98 4.94

(−0.29, 0.98) (−0.34, 0.94) (−0.36, 0.93) (−0.36, 0.93)

−1.98 −3.26 −2.86 −2.97

(0.84, 0.54) (0.72, 0.70) (0.79, 0.61) (0.66, 0.77)

−4.51 −3.21 −2.29 −2.15

16

N.J. Inkson et al. / J. Non-Newtonian Fluid Mech. 156 (2009) 7–20

Fig. 11. Contours of the first normal-stress difference for U = 0.1 ms−1 .

Fig. 12. Contours of the first normal-stress difference for U = 0.5 ms−1 .

Fig. 13. Contours of the first normal-stress difference for U = 0.75 ms−1 .

Fig. 14. Contours of the first normal-stress difference for U = 1.64 ms−1 .

region of negative shear stress is located at the rear of the cylinder (see Figs. 8–10). It can be seen from Table 5 that the magnitude of the shear stress decreases with increasing shear rates due to the extreme shear-thinning nature of LDPE shear rheology. Contour plots of the first normal-stress difference are shown in Figs. 11–14 and are quite different from the shear stress. The minimum value of the first normal-stress difference occurs at the front of the cylinder, and local maxima occur near the top of the cylinder, the channel wall opposite the cylinder and downstream. A pronounced wake of first normal-stress difference extends downstream behind the cylinder (see Fig. 11), which grows in magnitude with increasing velocity (see Figs. 12–14). The largest values of the first normalstress difference are shown with position co-ordinates, in Table 6, along with the x-position of the maximum stress in the wake. The behaviour of the first normal-stress difference is interesting as there

Table 6 Maximum dimensionless values of the first normal-stress difference on the cylinder, in the wake and on the wall as a function of U U (ms−1 )

0.1 0.5 0.75 1.64

Cylinder

Wake

Wall

max

(x, y)

max

x

max

12.29 14.00 11.68 11.75

(0.29,0.96) (0.28,0.96) (0.14,0.99) (0.23,0.98)

2.61 5.96 6.79 6.41

1.47 1.55 1.50 1.55

5.84 8.97 7.04 7.36

is an increase in magnitude between U = 0.1 and 0.5 ms−1 which is due to the strain-hardening features of LDPE extensional rheology. This is due to the branched molecules being highly stretched in the region where the melt accelerates between the cylinder and the channel walls. After this initial increase in first normal-stress difference, the strain rates then attain values beyond those corresponding to the peak in extensional viscosity and the normal-stress difference then decreases in magnitude at the higher velocities. This behaviour can be attributed to the collapsing of the molecules into the tubes during the process of branch point withdrawal. Fig. 15 shows contours of the second normal-stress difference for U = 0.75 ms−1 . The second normal-stress difference is highest at the leading edge of the cylinder with a maximum at the front corner. There exists a circular isobars of second normal-stress difference downstream of the cylinder. There is a region of negative second normal-stress difference. For high characteristic velocities the symmetry fore and aft the cylinder is broken and the wake of stretched molecules extends well past the cylinder, beyond x = 5 (see Figs. 16 and 17). A boundary layer of stretched molecules is seen to develop at the front of the cylinder ( 2.2) which extends around and behind the cylinder. This boundary layer grows in thickness with increasing velocity. The maximum stretch, reached at U = 1.64 ms−1 , is  6.16 (see Fig. 17). The maximum possible stretch attainable for this slow mode is 10.8. We would only see stretching to this degree in a geometry with a stronger extensional component, for example in a

N.J. Inkson et al. / J. Non-Newtonian Fluid Mech. 156 (2009) 7–20

17

Fig. 15. Contours of the second normal-stress difference for U = 0.75 ms−1 .

Fig. 16. Contours of the stretch corresponding to the slowest relaxing mode for U = 0.75 ms−1 .

Fig. 17. Contours of the stretch corresponding to the slowest relaxing mode for U = 1.64 ms−1 .

Table 7 Maximum values of the dimensionless stretch of the slowest mode on the cylinder, in the wake and on the wall as a function of U U (ms−1 )

0.1 0.5 0.75 1.64

Cylinder

Wake

Wall

max

(x, y)

max

x

max

1.90 3.94 4.73 6.16

(0.00,1.00) (0.04,1.00) (0.21,0.98) (0.38,0.93)

1.12 3.28 4.41 6.07

1.50 1.62 1.60 1.67

1.70 3.37 3.83 5.67

contraction geometry. Detailed information on the maximum values of the molecular stretch for different velocities are given in Table 7. We can examine the stretch for each mode by plotting it against the distance around the cylinder wall, ranging from  = 0, upstream to  =  downstream. For the case of U = 0.5 ms−1 , this information is plotted in Fig. 18(a). The corresponding shear stress (dashed line) and first normal-stress difference (solid line) are displayed in Fig. 18(b). The stretch is computed using the algebraic Eq. (23)

which involves the trace of the extra-stress tensor. The derivatives of the components of the extra-stress tensor will, in general, be discontinuous across spectral element boundaries even though the extra-stress components are approximated by continuous spectral element approximations. This is the explanation for the kinks in Fig. 18(a). Finally, a quantitative comparison is performed between the numerical simulations performed here and the experimental results presented by Verbeeten et al. [27] for U = 1.64 ms−1 . The data from the paper of Verbeeten et al. [27] has been scanned. First, we compare with the velocity profiles at the axial locations x = −4, −2.5, −1.5, 1.5, 2.5 for U = 1.64 ms−1 , which corresponds to an average Weissenberg number of 2.9. The dimensionless velocity is represented along the horizontal co-ordinate as in the paper of Verbeeten et al. [27]. The cross-sectional velocities across the channel are plotted in Fig. 19(a) and good agreement is found between the simulations and experiments. The horizontal component of the velocity along the centreline is plotted in (y = 0) in Fig. 19(b) and agrees well with the data. There is a difference between the channel width/cylinder ratio used in the simulation and experiments, the

Fig. 18. (a) Stretch of three modes vs. distance along the cylinder, (b) shear stress (dashed line) and first normal-stress difference (solid line) vs. distance along the cylinder, for U = 0.5 ms−1 .

18

N.J. Inkson et al. / J. Non-Newtonian Fluid Mech. 156 (2009) 7–20

Fig. 19. (a) Dimensionless velocity profiles at various axial locations (b) horizontal component of the velocity along the centreline (dots are scanned data form Verbeeten et al. [27] and the lines are from simulation) when U = 1.64 ms−1 .

Fig. 20. Principal stress difference for U = 1.64 ms−1 compared with experimental data from Verbeeten et al. [27].

Fig. 21. Contours of shear stress with ˛i = 0 for U = 0.5 ms−1 .

ratio of channel half-width to cylinder radius is 2.0:1 in the numerical calculations and 2.084:1 in the experiments. Some allowance has been made for this in the comparisons. Experimentally, the stress in polymer melt flow can be measured using the technique of flow induced optical birefringence [9]. Linearly polarised monochromatic light is passed through quarterwave plates and received through an analyzer. The observed isochromate fringes represent integer differences in the principal refractive index, which in turn can be represented by contours of principal stress difference. The principal stress difference is

 PSD =

2 + ( −  )2 4xy xx yy

(50)

The stress-optical rule for 2D flows is used to find the stress difference between fringes:  =

 Cd

(51)

where C is the stress-optical coefficient for a given material, d is the length of the light path in the medium and  is the wavelength of the laser light. Fig. 20 shows contours of principal stress for U = 1.64 ms−1 compared with experimental data from Verbeeten et al. [27], using the same parameters. This corresponds to an average Weissenberg number of 2.9. There is good agreement between the numerical predictions and experimental data, particularly when noting the slight difference between the ratio of channel

Fig. 22. Contours of the first normal-stress difference with ˛i = 0 for U = 0.5 ms−1 .

N.J. Inkson et al. / J. Non-Newtonian Fluid Mech. 156 (2009) 7–20

19

a multi-mode XPP fluid past a cylinder has been simulated and the method is stable over an extended range of Weissenberg numbers (corresponding to the slowest relaxing mode) up to around 7.529. The drag on the cylinder has been determined as a function of the characteristic velocity of the problem. The drag decreases by a factor of two over the range of Weissenberg numbers considered and appears to have reached either a limiting or a minimum value. Qualitative and quantitative comparisons with experiments have been performed and good agreement is found with the numerical simulation confirming the ability of the XPP model to predict the flow of commercial polymer melts. The choice of an alternative constitutive model of pom-pom type should allow the spectral element method outlined in this paper to achieve flow simulations at higher Weissenberg numbers in the future. The construction of even more efficient preconditioners will also make the approach less computationally expensive. Fig. 23. Shear stress (black) and first normal-stress difference (grey) for DSM LDPE with non-zero ˛i (solid) and with ˛i = 0 (dashed) for U = 0.5 ms−1 .

half-width to cylinder radius which is 2.0:1 in the numerical calculations and 2.084:1 in the experiments. Similar flow features such as the small island in principal stress just downstream of the cylinder, the shape of the wake and the number of fringes upstream and near the channel wall are observed. This comparison demonstrates that the XPP model can be used to predict the flow of this particular LPDE. Finally, the effect of setting ˛i = 0, i = 1, 2, 3, for all modes, in order to produce a simulation corresponding to the DPP model, is investigated. This does not have a significant effect on the transient shear and elongational viscosity. In effect, it switches off the Giesekus-like term in the XPP equation. We have reason to believe that this parameter choice would not suffer from the effects of multiple solutions of the XPP model. It also has the effect of setting the second normal-stress difference to zero. Figs. 21 and 22 compare the shear stress and first normal-stress difference for U = 0.5 ms−1 when the DSM LDPE is modelled using the DPP model (˛i = 0) and / 0) XPP model. We can be compared to Figs. 8 and 12 for the (˛i = find that with this choice of parameters the calculated drag is now 80.45 for the DPP model, whereas with the XPP model the drag was 92.83. The reduction in drag is due to the stresses being of lower magnitude, on average, around the cylinder. The shear stress and first normal-stress difference on the surface of the cylinder are plotted in Fig. 23. 6. Conclusions An extension of a stabilized spectral element method to predict the flow of a branched polymer melt using a multi-mode extended pom-pom model has been described. A quantitative description of the rheology of LDPE has been used. The extended pom-pom model was first introduced to eradicate perceived deficiencies of the original pom-pom model such as the lack of a second normal-stress difference and a discontinuity in the derivative of the extensional viscosity. However, it is shown that the XPP model itself possesses some disconcerting attributes. In simple steady shear and uniaxial extensional flow, multiple solutions are found for many choices of the branching parameter, q, and the anisotropy parameter, ˛. A decoupled scheme is used in which the solution of the conservation equations is decoupled from the solution of the constitutive equation at each time step. The convection terms in the momentum and constitutive equations are treated using an operator-integration-factor-splitting technique. The conservation equations are solved using a nested preconditioned conjugate gradient method with efficient preconditioners. The complex flow of

Acknowledgements This research is funded by the Engineering and Physical Sciences Research Council of the United Kingdom through grant EP/C513037. We acknowledge the assistance of James Osborne for computational support in using Condor, a program which exploits spare computing capacity of open access work stations at Cardiff University, for performing some of the simulations in this paper. References [1] B. Bernstein, E.A. Kearsley, L.J. Zapas, A theory of stress relaxation with finite strain, Trans. Soc. Rheol. (1963) 391–410. [2] R.J. Blackwell, T.C.B. McLeish, O.G. Harlen, Molecular drag-strain coupling in branched polymer melts, J. Rheol. 44 (2000) 121–136. [3] N. Clemeur, R.P.G. Rutgers, B. Debbaut, On the evaluation of some differential formulations for the pom-pom constitutive model, Rheol. Acta 42 (2003) 217–231. [4] B. Debbaut, J. Dooley, Secondary motions in straight and tapered channels: experiments and three-dimensional finite element simulation with a multimode differential viscoelastic model, J. Rheol. 43 (6) (1999) 1525–1545. [5] P.G. de, Gennes, Reptation of a polymer chain in the presence of fixed obstacles, J. Chem. Phys. 5 (2) (1971) 572–579. [6] M. Doi, S.F. Edwards, The Theory of Polymer Dynamics, Oxford University Press, Oxford, 1986. [7] H. Giesekus, Die Elastizität von Flüssigkeiten, Rheol. Acta (1996) 29–35. [8] H. Giesekus, A simple constitutive equation for polymer fluids based on the concept of deformation-dependent tensorial mobility, J. Non-Newtonian Fluid Mech. 11 (1982) 69–109. [9] H. Janeschitz-Kriegl, Polymer Melt Rheology and Flow Birefringence, SpringerVerlag, Heidelberg, 1983. [10] N.J. Inkson, T.C.B. McLeish, O.G. Harlen, F.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. [11] N.J. Inkson, T.N. Phillips, Unphysical phenomena associated with the extended pom-pom model in steady flow, J. Non-Newtonian Fluid Mech. 145 (2007) 92–101. [12] Y. Maday, A.T. Patera, Spectral element methods for the incompressible Navier–Stokes equations, in: A.K. Noor, J.T. Oden (Eds.), State of the Art Surveys in Computational Mechanics, ASME, New York, 1989, pp. 71–143. [13] Y. Maday, A.T. Patera, E.M. Rønquist, An operator-integration-factor-splitting method for time-dependent problems: application to incompressible flow, J. Sci. Comput. (1990) 263–292. [14] T.C.B. McLeish, R.G. Larson, Molecular constitutive equations for a class of branched polymers: the pom-pom polymer, J. Rheol. 42 (1998) 81–110. [15] S.T. Milner, T.C.B. McLeish, Parameter-free theory for stress relaxation in star polymer melts, Macromolecules 30 (1997) 2159–2166. [16] J.G. Oldroyd, On the formulation of a rheological equation of state, Proc. R. Soc. Lond. A 200 (1950) 523–541. [17] R.G. Owens, C. Chauvière, T.N. Phillips, A locally-upwinded spectral technique (LUST) for viscoelastic flow, J. Non-Newtonian Fluid Mech. 108 (2002) 49–71. [18] N. Phan-Thien, A nonlinear network viscoelastic model, Trans. Soc. Rheol. 22 (1978) 259–283. [19] N. Phan-Thien, R.I. Tanner, A new constitutive equation derived from network theory, J. Non-Newtonian Fluid Mech. 2 (1977) 353–365. [20] R.G. Owens, T.N. Phillips, Computational Rheology., Imperial College Press, London, 2002. [21] J. van Meerfeld, Note on the thermodynamic consistency of the integral pompom model, J. Non-Newtonian Fluid Mech. 108 (2002) 291–299.

20

N.J. Inkson et al. / J. Non-Newtonian Fluid Mech. 156 (2009) 7–20

[22] R.G.M. van Os, T.N. Phillips, The prediction of complex flows of polymer melts using spectral elements, J. Non-Newtonian Fluid Mech. 122 (2004) 287–301. [23] R.G.M. van Os, T.N. Phillips, Spectral element methods for transient viscoelastic flow problems, J. Comput. Phys. 201 (2004) 286–314. [24] R.G.M. van Os, T.N. Phillips, Efficient and stable spectral element methods for predicting the flow of an XPP fluid past a cylinder, J. Non-Newtonian Fluid Mech. 129 (2005) 143–162. [25] M. Aboubacar, J.P. Aguayo, P.M. Phillips, T.N. Phillips, H.R. Tamaddon-Jahromi, B.A. Snigrev, M.F. Webster, Modelling pom-pom type models with highorder finite volume schemes, J. Non-Newtonian Fluid Mech. 126 (2005) 207–220.

[26] W.M.H. Verbeeten, G.W.M. Peters, F.T.P. Baaijens, Differential constitutive equations for polymer melts: the extended pom-pom model, J. Rheol. 45 (4) (2001) 823–843. [27] W.M.H. Verbeeten, G.W.M. Peters, F.T.P. Baaijens, Viscoelastic analysis of complex melt flows using the extended pom-pom model, J. Non-Newtonian Fluid Mech. 108 (2002) 301–326. [28] W.M.H. Verbeeten, G.W.M. Peters, F.T.P. Baaijens, Numerical simulations of the planar contraction flow for a polyethylene melt using the XPP model, J. NonNewtonian Fluid Mech. 117 (2004) 73–84.