Extensional response of the pom-pom model through planar contraction flows for branched polymer melts

J. Non-Newtonian Fluid Mech. 134 (2006) 105–126

Extensional response of the pom-pom model through planar contraction flows for branched polymer melts夽 J.P. Aguayo, H.R. Tamaddon-Jahromi, M.F. Webster∗ Institute of Non-Newtonian Fluid Mechanics, Department of Computer Science, University of Wales, Swansea SA2 8PP, United Kingdom Received 28 June 2005; received in revised form 27 November 2005; accepted 1 December 2005

Abstract One objective of this study is to assess the influence of the number of dangling arms at each end of the pom-pom molecule on flow characteristics in a 4:1 planar rounded-corner contraction, where varying the number of arms affects the level of entanglement in the system. For these complex flows, we evaluate the major influence of extensional viscosity and Trouton ratio when comparing kinetic-based as opposed to phenomenological networkbased models. The stability of the proposed extended pom-pom model (XPP) is investigated for a range of Weissenberg numbers utilising a hybrid finite-element/volume scheme. For this scheme, we solve the momentum-continuity equations by a fractional-staged Taylor–Galerkin/pressurecorrection procedure and invoke a cell-vertex fluctuation distribution finite volume scheme for the constitutive law. Distinction may be drawn between fluid response in the extension-hardening regime and those displaying a degree of softening. In particular, trends in salient-corner vortexintensity response can be associated with the extensional properties of the extended version of the pom-pom model. Critical levels of solution elasticity attainable depend on precise extensional viscosity characteristics. The influence of various model parameters upon field response is described through rates of deformation, stress and stretch. This would include consideration for parameters governing anisotropy (α), and the ratio of stretch to orientation relaxation times (). © 2006 Elsevier B.V. All rights reserved. Keywords: Pom-pom; Hybrid finite element/volume; Rounded-corner contraction; Viscoelasticity; Strain-hardening; Branched molecular structure

1. Introduction This article is intended to provide an extension to our earlier studies [1] for steady-state solutions with 4:1 abrupt contraction flows using the single-equation representation of the pompom model. The selected constitutive models must be capable of reproducing the rheological properties of the fluid in wellestablished rheometrical flows, at least in a qualitatively sense. Such models are candidates for use in the investigation of more complex flows, where the challenge shifts towards representing realistic fluid response, and suitable boundary conditions becomes a major issue. In developing the proposed hybrid finite-element/volume algorithm (hy fV), we have made extensive use of the wellestablished UCM-like constitutive models. This includes use of the Oldroyd-B model (constant shear-viscosity), for benchmark-

夽 Paper presented at the AERC 2005 Conference held in Grenoble, France. ∗

Corresponding author. Tel.: +44 1792 295656; fax: +44 1792 295708. E-mail address: [email protected] (M.F. Webster).

0377-0257/$ – see front matter © 2006 Elsevier B.V. All rights reserved. doi:10.1016/j.jnnfm.2005.12.006

ing and to resolve highly-elastic flow configurations, and a range of Phan-Thien–Tanner variants, that emulate the rheology of low and high density polymer melts ([2,3]). More recently, attention has turned to the study of transient flows, and the pursuit of highorder transient schemes [4,5]. Key aspects to achieving this goal are: a consistent treatment of the time-derivative in the stress constitutive equation, coupled with a novel weighting function for fv-cell and median-dual-cell-based contributions. To date, application of these ideas to a start-up flow of an Oldroyd-B fluid through a planar contraction has exposed some provocative transient flow features. This is in advance of its asymptotic steady-state, which is comparable to that reported in the literature [6]. Application to contraction flows, in both planar and axisymmetric configurations, has demonstrated the capacity of the proposed scheme to reproduce qualitatively many distinctive experimental flow features. For example, with constant viscosity models, vortex inhibition has been observed in planar contraction flows, as opposed to vortex enhancement in their axisymmetric counterparts [2]. The size and shape of upstreamvortices for shear-thinning fluids, is mainly governed by levels of strain-hardening and high Trouton ratio, following closely

106

J.P. Aguayo et al. / J. Non-Newtonian Fluid Mech. 134 (2006) 105–126

that observed experimentally for a range of LDPE and LLDPE blends [7]. As a consequence, we are lead naturally to a quantitative investigation for contraction flows, which may play a significant role within many industrial processes. Quantitative comparison against experiments requires, first and foremost, the use of constitutive models that can reproduce the steady rheological properties of the fluid in both pure shear and extensional flows (uniaxial and planar). For polymer melts, phenomenological equations (PTT, Giesekus, K-BKZ) are unable to achieve this over a significant range of deformation-rates. Although PTT-like models can account for both uniaxial and planar strain-hardening, shear-viscosity also reflects a dependence on the same parameters that govern extensional properties, see Matallah et al. [8]. This renders difficulty in fitting both viscosity functions independently. Giesekus models exhibit sustained strain-hardening response with increasing strain-rates (reaching a plateau), which is particularly appropriate experimentally for polymer melts. Typically, strain-softening is manifest only after a region of pronounced strain-hardening. Recently, a modified version of the integral K-BKZ model using a PSM damping function was employed by Mitsoulis et al. [9] to fit the rheological data for two commercial LDPE melts, that included a match to planar strain-hardening properties. Using eight modes, these authors reported quantitative agreement between numerical and experimental observations, with decreasing vortex-lengths for a melt of low Trouton-ratio. Nevertheless, and for the same melt, qualitative agreement was not achieved for the Bagley entry pressurecorrection. In contrast and for a strain-hardening melt with larger Trouton-ratio, computations reproduced experimental trends in both vortex-size and entry correction. This was an impressive achievement indeed! New constitutive equations have been specifically derived to describe the complex rheological behaviour of polymer melts with long side-branches, such as applicable for low density polyethylene. These are termed pom-pom equations and have been extracted in integral–differential form, see McLeish and Larson [10]. Verbeeten et al. [11] provided an extension to the original pom-pom model (XPP) on the basis of reptation dynamics and a simplified branched molecular structure. The pompom constitutive equation predicts the rheological behaviour well of LDPE-melts in both shear and extension [12]. Inkson et al. [13] reported predictions whilst representing a hypothetical (LDPE) melt via a pom-pom model, contrasting the use of different numbers of arms. These authors found that multi-mode pom-pom equations with a physically reasonable distribution of branching, could predict LDPE rheology over four decades in deformation-rate, covering three different flow problems. In Zatloukal [14], fitting and prediction of the proposed Leonov model was compared with the XPP and modified White–Metzner models, for steady shear and uniaxial extensional and three polymer melts (LDPE, mLLDPE, PVB, used mainly in melt-film production). There, the XPP model was shown to provide excellent predictions for steady shear and elongational viscosity, particularly at lower levels of deformation-rate. In addition, Bogaerds et al. [15] investigated the stability of the proposed pom-pom model of Verbeeten et al. [11] for simple viscometric shear flow

using both a one-dimensional spectral eigenvalue analysis and two-dimensional finite element computations. Bogaerds et al. found that the inclusion of a relatively small second normal stress-difference contribution (α = 0.1) had a stablizing effect on steady simple shear flow. They also reported that despite the similarity between the material functions of the EPTT and the XPP model, the dynamic response to perturbation with the XPP fluid was considerably different. Furthermore, Tanner and Nasseri [16] investigated the viscometric response of several constitutive models in both shear and extension. There, a new form of single-mode PTT model, named PTT–XPP, was shown to practically reproduce the steady elongational response of the XPP model of Verbeeten et. al [11,17], see also our Fluid-(III) below. This PTT–XPP model neglected the Giesekus term (with α = 0) and the (f (λ, τ) − 1)I term of the XPP-model, generalising the PTT{g}-term to the f (λ, τ)-term of the XPP-form (see below for terms and notation). It was noted that this PTT–XPP model differed from the XPP model through its shear viscosity at larger shear-rates. There, the XPP-form displayed a greater tendency to shear-thin when compared to the functionality of the PTT–XPP model. A number of other computational studies have contributed within this general area utilising a variety of different approaches and versions of pom-pom models. Bishko et al. [18] performed transient flow simulations for a 4:1 planar contraction with a single-mode differential pom-pom model. A dependence of vortex-size was identified on the degree of branched molecular structure with increasing Weissenberg number. In a similar manner, Clemeur et al. [19] demonstrated the performance of the double convected pom-pom model (DCPP), comparing simulation results for the abrupt contraction flow with experimental data. These authors observed a deviation of 10–15% between experimental data and numerical predictions. In addition, Verbeeten et al. [17,20] investigated the performance of the multimode XPP model (single and double equation) for transient flow around a cylinder, through a cross-slot device [17] and in a planar 3.29:1 contraction flow [20]. There numerically, a Discrete Elastic Viscous Stress Splitting technique, in combination with a Discontinues Galerkin method (DEVSS/DG) was employed. They observed the occurrence of negative molecular backbone stretch at the front and back stagnation point stations in flow around a cylinder and near the re-entrant corner of the contraction. Numerical convergence was secured by modifying the XPP model in which the stretch dynamics were adjusted [20], in agreement with a modification proposed by Van Meerveld [21] based on non-equilibrium thermodynamics (see on). Recently, a steady-state solution for a three-dimensional contraction geometry has been reported by Sirakov et al. [22], using a multi-mode pom-pom model and including comparison against experimental data. Findings confirmed that numerical predictions on upstream growth of vortex size agreed with experimental evidence to within a discrepancy of some 15%. Wapperom and Keunings [23] switched flow problem to study the behaviour of the pom-pom model through a planar 4:1:4 constriction flow with rounded-corners. These authors reported non-physical solution behaviour at high Weissenberg number, where the value of the stretch parameter (λ) became smaller

J.P. Aguayo et al. / J. Non-Newtonian Fluid Mech. 134 (2006) 105–126

than unity downstream of the constriction (see Verbeeten et al. [17,20] likewise, noting a resolution on this issue as proposed in [20]). In the present study, we focus upon the influence of the number of dangling arms at each end of the pom-pom molecule and on the corresponding flow characteristics generated in a 4:1 planar rounded-corner contraction. Varying the number of arms in this way affects the degree of entanglement associated with the system as a whole and significantly adjusts extensional properties in such a complex flow setting. In addition, we are able to compare and contrast flow response for these kinetic-based pompom models against their network-based PTT counterparts. By matching on extensional viscosity and Trouton ratio peaks across these two classes of models, we have been able to differentiate between the precise role and influence that these quantities have upon solution response in non-viscometric flows. 2. Governing equations and constitutive modelling We consider isothermal, incompressible creeping flow. Under such conditions, mass and momentum conservation equations can be expressed, through conventional operator notation and independent space–time variables (x, t), as: ∇ · u = 0,

(1)

and ρ

∂u = −ρu · ∇u + ∇ · (2µs d + τ) − ∇p. ∂t

(2)

Here, u, d, and p represent fluid velocity, rate-of-deformation (∇u + (∇u)T )/2, and hydrodynamic pressure, respectively. Stress is split into a solvent component of viscosity (µs ), and a polymeric contribution τ, governed by its polymeric viscosity (µe ). Fluid density is denoted by ρ. The constitutive equation for the single equation version of the XPP model may be expressed, using the identity tensor (I) in the form, ∇

f (λ, τ)τ + λ0b τ +G0 (f (λ, τ) − 1)I +

α τ · τ = 2G0 λ0b d, G0 (3)

with the backbone stretch approximated through the function f (λ, τ), utilising the trace operator (I),  
 λ0b ν(λ−1) 1 αIτ·τ 1 f (λ, τ) = 2 1− e , (4) + 2 1− λ0s λ λ 3G20 ∇

and where τ denotes the upper-convected material derivative of τ, viz. ∇

τ=

∂τ + u · ∇τ − (∇u)T · τ − τ · (∇u). ∂t

(5)

In this notational form, λ0b and λ0s represent the orientation and backbone stretch relaxation time-scales, respectively, and G0 is the linear relaxation modulus. For the Single eXtended PomPom (SXPP) model, backbone stretch is coupled directly to the

trace of the extra-stress (Iτ ), viz.  |Iτ | λ= 1+ . 3G0

107

(6)

For general categorisation, we appeal to conventional definitions of group numbers, in Reynolds and Weissenberg numbers, defined as: λ0b U ρUL , We = . Re = µ L We also find it convenient to refer to specific material parameter ratios on viscosity fractions and relaxation time-scales through parameters β and , identified as: β=

µs , µs + µ e

=

λ0s . λ0b

Above, (L) and (U) are characteristic scales on length and velocity, respectively; µ = µs + µe represents the total viscosity. To preserve similarity between the various forms of constitutive laws employed, such as Oldroyd-B, PTT, and SXPP models, we define, µe = G0 λ0b .  Iτ Note that negative values for the argument of 1 + 3G may 0 be avoided by adopting Eq. (6), based on the absolute value of (Iτ ), see Verbeeten et. al [17]. The parameter ν in (4) was incorporated into the model by Blackwell et al. [24] to remove discontinuity in the derivative of extensional viscosity, present in the differential approximation of the original pom-pom model. Its value was estimated by data-fitting and found to be inversely proportional to the number of arms (q). More precisely following Blackwell et al., we have taken ν=

2 . q

Further modifications to the pom-pom model have been advocated by Van Meerveld [21] and used by Verbeeten et al. [20]. The mathematical nature of these adjustments is to alter the multiplicative factor on the exponential term of Eq. (4) from {1 − λ−1 } (XPP) to {1 − λ−2 } (mXPP). This has some impact on the extensional rheology for the material parameters chosen, (see mXPP in Fig. 1) over the decade of extension rates [100 , 101 ]. For q = 2, there is initial increase in strain-softening, then reversal of this behaviour, before returning to similar softening response. Equivalently, Trouton ratios are lessened. For q = 15, the extensional viscosity indicates first a drop prior to extension rates of O(1), then increases to the same level of hardening as for the {1 − λ−1 } case. This initial decrease in extensional viscosities is not discernible in the {1 − λ−1 } instance. To get some appreciation for the relative significant of these model modifications, we can comment that at the level of elasticity of We = 20 (q = 5, Re = 0), stress components and backbone stretch can be reduced in peak corner values, by as much as O(40%) in stress and O(20%) in (λ). Ultimately, this may lead to enhanced levels of Wecrit steady-state solutions, see Verbeeten et al. [20]. Nevertheless, a significant issue in this regards is the influence of α = 0, introduced to gather some second normal stress-difference effects. We note, as others, that such parameterisation can introduce analytical singularity into the solution for

108

J.P. Aguayo et al. / J. Non-Newtonian Fluid Mech. 134 (2006) 105–126

sufficiently large α and with large enough number of branchedarms (q), irrespective of f (λ, τ)-form [25] (see on, Fig. 16). Such singularity appears at large shear-rates, being a direct consequence of the incorporation of the Giesekus-type parameter (α). In steady elongation for (α ≥ 0.5), we may observe an intermediate range of deformation-rates around [101 , 102 ] where a plateau emerges, prior to ultimate softening. We include some (α = 0) solutions in the present work to shed some light upon this aspect. To provoke variation in extensional properties, we give consideration to a particular version of the PTT network-based model: the exponential form (EPTT) [26,27], for which the constitutive equation may be expressed as: ∇

gτ + λ1 τ = 2µe d, where

(7)



 εptt λ1 g = exp trace(τ) , µe

(8)

with λ1 the relaxation time for the PTT model. The model parameter (εptt ≥ 0) is a unitless material parameter that can be evaluated by fitting to experimental data. For this study as common elsewhere, two values for parameter εptt , 0.02 and 0.25 are considered, see below. All variants of the EPTT model exhibit shear-thinning properties. Steady elongational viscosities remain finite for these PTT models, displaying alternative exposure to severity of hardening at moderate deformation-rates, whilst ultimately softening in the limit of large deformationrates. The intensity of the hardening is governed by the εptt -parameter, with stronger hardening corresponding to lower

values of εptt . In the asymptotic limit of εptt tending to zero, the infinite extensional dynamics of the Oldroyd model is recovered. 3. Material functions In previous articles, we have investigated the rheological response for the SXPP model in planar channel flows [28] and in a planar 4:1 abrupt contractions [1]. Parameters that influence the material functions, and in consequence, the flow response are: the ratio of relaxation times (stretch/orientation) of the backbone, represented by the non-dimensional parameter (); solvent to total viscosities ratio (β); the number of dangling arms at both ends of the molecule (q); the anisotropy-parameter (α). The SXPP model represents a shear-thinning fluid of hardening-softening form under extensional deformation. Fig. 1 illustrates the adjustment in rheological properties of the model when the number of arms (q) are varied. At this selected ratio of solvent:total viscosity, (β = 1/9), increasing the number of arms has an insignificant impact upon the shear viscosity. In contrast, this has a major influence upon planar extensional viscosmetric response. By increasing the number of arms from q = 2 to 15, the peak-value of the extensional viscosity rises almost one decade in the log–log plot. This peak occurs at a strain-rate of around 10–20 units for q ≥ 5. With respect to Trouton ratio (also planar), such an increment in arms-q produces an increase of one decade in maximum value. This can be explained in terms of the reduction in shear viscosity and the increasing influence of the extensional viscosity. The influence of q-variation on the first normal stress coefficient (ψ1 ) is negligible at low to moderate

Fig. 1. Rheological properties of XPP model, variation in q:  = 1/3, β = 1/9 and α = 0.15.

J.P. Aguayo et al. / J. Non-Newtonian Fluid Mech. 134 (2006) 105–126

109

Fig. 2. Transient response in planar extension and shear flow of XPP model, variation in q:  = 1/3, β = 1/9 and α = 0.15.

shear-rates (so up to five units in shear). At higher shear-rates, greater numbers of arms on the pom-pom molecule result in a decrement in rate-of-decline in ψ1 (Fig. 1c). In addition, adjusting the number of arms from 2 to 15 causes a significant increment in backbone stretch, λ, governing the length of the molecule. This becomes apparent in Fig. 1d when deformationrates exceed 0.5 units. To illustrate start-up and trends towards a terminating configuration, the SXPP transient response is presented in Fig. 2 for planar extensional and simple shear flows with corresponding backbone stretch. Time is scaled with the relaxation time for orientation of the backbone (λ0b ); transient viscosities (planar and shear) are scaled with the steady shear viscosity at zero shearrate (η0 ). These variables (ηPlanar , ηShear , λ) reach their steady state values more rapidly for larger fixed values of deformationdγ dγ d rate ( d dt or dt ). Additionally, (λ)-results for dt and dt at levels of 0.01 and 1, do not show any significant overshoot. For larger deformation-rates this feature is clearly present, in particular, for larger numbers of arms in the molecule. Note that, differences in stress and stretch fields, caused by unsteady Lagrangian influences, are not particularly significant here, because of the short residence time in the contraction flow zone. This is illustrated in Fig. 9, where the relaxation of stress and stretch is shown to be similar for different values of (q) at the plane-of-symmetry. One notes, however that elevation in (q) stimulates elevation in {N1 , λ} at the vicinity of the corner, which governs the respective plateau-levels achieved along the downstream wall, as displayed in Fig. 11. Hence, in this respect, we can detect significant impact of the larger numbers of arms.

4. Hybrid finite element/volume method The hybrid finite element/volume scheme is a semi-implicit, time-splitting, fractional-staged formulation, that invokes finite element discretisation for velocity–pressure parts of the system and finite volume for stress, see [5,29]. Under the feconstruction, a two-step Lax–Wendroff splitting is adopted, a Taylor–Galerkin variant [30,31], alongside an incremental pressure-correction procedure (with 0 ≤ θ1 ≤ 1, see Guermond and Quartapelle [32]). With a forward time increment factor θ2 = 0.5, such a pressure-correction implementation takes the scheme through to second-order temporal accuracy under incompressible conditions (see [33,34]). The incremental form of pressure-correction (θ1 > 0) proves to be superior in uniform temporal error bounds over its non-incremental counterpart (θ1 = 0). Appealing to semi-discrete representation on the time step [t n , t n+1 ], with starting values [un , τ n , pn , pn−1 ] for the present creeping flow setting, the three-stage structure can be expressed as [35]: • Stage 1a: 2Re n+(1/2) − un ) t (u − ∇(pn + θ1 {pn

n+(1/2) n = [∇τ]n + ∇ · 2β (d 2 +d )

− pn−1 }) + Fgn ,  2We n+(1/2) n ) = − Weu · ∇τ + 2(1 − β)d (τ − τ ta 1−β We [f (λ, τ) − 1]I n · ∇u + τ · ∇ut } ,

− f (λ, τ)τ − + We{τ

−

αWe 1−β τ

·τ

110

J.P. Aguayo et al. / J. Non-Newtonian Fluid Mech. 134 (2006) 105–126

• Stage 1b: Re ∗ t (u



− un ) = [∇τ]n+(1/2) + ∇ · 2β (d

∗ +d n )

Q=



2 n+(1/2) −∇(pn + θ1 {pn − pn−1 }) + Fg ,

We n+(1/2) n ) = −Weu · ∇τ + 2(1 − β)d (τ − τ ta αWe −f (λ, τ)τ − 1−β We [f (λ, τ) − 1]I− 1−β τ · τ  t n+(1/2)

+ We{τ · ∇u + τ · ∇u }

,

• Stage 2: ∇ 2 (pn+1 − pn ) =

Re ∇ · u∗ , θ2 t

• Stage 3: 2Re n+1 − u∗ ) = −θ2 ∇(pn+1 − pn ). (u t

(9)

Galerkin discretisation may be applied to the Stokesian components of the system; the momentum equation at Stage 1, the pressure-correction step at Stage 2 and the incompressibility correction constraint at Stage 3. For reasons of accuracy, the resulting Galerkin Mass matrix-vector equations at Stages 1 and 3 are solved using an efficient element-by-element Jacobi scheme, requiring only a handful of iterations (here, five chosen). With only a single matrix reduction phase necessary at the outset (t = 0), a direct Choleski decomposition procedure is invoked to handle Stage 2. Finally, diffusion terms at Stages 1a and 1b are treated in a semi-implicit manner in order to enhance stability. Pressure temporal treatment calls upon multi-step reference across three successive levels [t n−1 , t n , t n+1 ]. Stress is computed at the vertices of the fv-triangular subcells of the parent fe-triangle, constructed by connecting midside nodes of the parent triangle (see Fig. 3a). To outline our fv-discretisation, first let us reorganise the stress constitutive Eq. (3) into flux-source components, ∂τ = −R + Q, ∂t

(10)

1−β 1 [2(1 − β)d − f (λ, τ)τ] − [f (λ, τ) − 1]I We (We)2 α τ · τ + (∇u)T · τ + τ · (∇u), − (11) 1−β

where (R = u · ∇τ) and (Q) are the corresponding flux and source terms, respectively. Next, we consider each scalar stress components, τ, acting on an arbitrary volume Ω. Variation in τ is controlled through fluctuation in these components of flux and source. One must evaluate flux and source integral contributions on each fv-subcell and its associated Median-Dual-Cell. Fluctuation distribution coefficients govern the proportion of combination from each fv-cell to each of its cell-vertices. The update for a typical node (l) is then constructed by summing contributions from  its control volume, Ωl , composed of such fv-triangles, {Ωl = ∀Tl ΩTl }, see Fig. 3b. An area-weighted time-consistent distributional approach results over both flux and source contributions that ensures second-order accurate steady solutions. The Median-Dual-Cell component is necessary to retain stability under complex flow settings, according to which the following generalised complete nodal-update scheme emerges (CT3 ), ⎡ ⎤ n+1   ˆ lT ⎦ τl ⎣ δT αTl ΩT + δTMDC Ω t ∀Tl

=

 ∀Tl

MDCl

δT αTl bT +



δTMDC blMDC .

(12)

∀MDCl

Here, bT = (−RT + β1 QT ), blMDC = (−RMDC + β2 QMDC )l ,  T ∀Tl represents all fv-cells surrounding node l and αl the fluctuation distribution coefficients, indexed per fv-cell T for node l. Boolean factors β1 and β2 are taken here as unity, and imbue a class of scheme variants. For contraction flows, in which source terms may dominate some sectors, the precise fluctuation scheme, which governs the choice of αl -factors, is that of low diffusion B (LDB). This option is consistent, linearity preserving and minimal in numerical diffusion. The last property lies in distinct contrast to a linear positive scheme such as the N-scheme [35,36]. The LDB distribution coefficients (αl ) are defined per cell-node l in each fv-cell T, and depend upon a, the

Fig. 3. hy fV spatial discretisation: (a) fe-cell with four fv- sub-cells and FD per T, (b) fv-control volume for node l, with Median-Dual-Cell (MDC-shaded).

J.P. Aguayo et al. / J. Non-Newtonian Fluid Mech. 134 (2006) 105–126

111

magnitude of the cell advection velocity a, and the angles it subtends at an inflow vertex (see [36]). The authors retained δT = κ, with κ = ξ/3 if |ξ| ≤ 3 and 1 otherwise, with δMDC = 1 − δT . Here, ξ = We(a/ h), with h the area square-root for the fv-cell in question. 5. Flow through a 4:1 planar contraction 5.1. Problem specification Here, the benchmark flow problem is creeping flow within a 4:1 planar, rounded-corner contraction geometry. Fig. 4 displays the set of unstructured triangular meshes employed, through increasing refinement in mesh density. Table 1 provides the corresponding mesh characteristics. The length of upstream and downstream sections are 27.5L and 49L, respectively. To solve the governing system of partial differential equations for the XPP model it is necessary to impose suitable boundary and initial conditions. However, unlike the situation for UCM/OldroydB models analytical expressions for fully-developed velocity and extra-stress XPP-profiles are not readily available. Instead, such profiles must be determined computationally by solving the equivalent planar entry-channel flow problem (see Aboubacar et al. [28] for details). Fully developed boundary conditions are established at the outflow ensuring no change with respect to velocity component Ux and vanishing component Uy . In addition, no-slip boundary conditions are imposed along the stationary walls. Once flow kinematics have been established at the inlet, stress can be gathered by solving a set of pointwise ODEs. The number of dangling molecular-arms is varied across the range 2 ≤ q ≤ 15, for which computations have been performed using one particular set of XPP parameters, namely: Re = 0,  = 1/3, α = 0.15 and β = 1/9. Over the number of arms employed (q), we may classify each fluid through notation, Fluid-{qi } where i = {2, 5, 10, 15}. This is our default setting with fixed α = 0.15. Beyond this, we also consider some parameterisation in α-setting. Furthermore, we contrast pom-pom solutions with those of PTT by considering adjustment in the {q, }-pairing. Time-step termination is ensured when the L2 -norm relative maximum difference between solution approximations over two successive time steps falls below a set threshold of (10−7 ). Simulations are initiated from a We = 0.1 solution state and continuation in {We} parameter-space is employed until a critical level (Wecrit ) is reached. We have found it convenient to catalogue results through measures of upstream vortex cell-size (X) and vertical-wall/contraction-plane vortex cell-size (L), following reference [1]. Here, in our plots, we recalibrate graphical scales for ready interpretation of data, zooming all values by twice the upstream channel half-width. Table 1 Mesh characteristics Meshes

Elements

Nodes

Degrees of freedom (u, p, τ)

Rmin

(4 : 1)a (4 : 1)b (4 : 1)c

1086 1626 2693

2325 3433 5652

14570 21502 35392

0.0296 0.0170 0.0097

Fig. 4. Unstructured fe-triangular meshing in contraction zone: (a) coarse, (b) medium, (c) refined.

6. Pom-pom solutions with q-variation (α = 0.15) One of the primary features of interest in our current solutions is the development of flow field structure under present rheology and with increasing elasticity. This is traced for each setting of arms-q, through vortex patterns for a series of steadystate solutions up to critical levels of Weissenberg number. In the light of this, we then proceed to analyse corresponding solution states in stress, deformation and stretch. With pom-pom modelling, the fresh new features are associated with the information on stretch of the polymer molecule within the flow domain. 6.1. Vortex behaviour and Wecrit Elasticity levels attainable in the numerical simulations have been found to decline with increasing number of dangling branched-arms (q), with Wecrit for q = {2, 5, 10 and 15 } being {60, 25, 13 and 12}, respectively. Previous investigations, for example, [37–39], have indicated that the highest level of elasticity attainable may depend on the level of strain-hardening experienced by the fluid. We observe that for (q = 2), there is almost no hardening, however it does not soften until a non-dimensional strain-rate of 2 units is reached. In this case, we gather the largest Wecrit of 60. For (q = 5), the maximum value in extensional viscosity is 50% larger than that for (q = 2) and this would leads us to expect a reduction in Wecrit compared to the non-hardening case (around Wecrit = 25). For the two further selected levels of {q = 10, Wecrit = 13} and {q = 15, Wecrit = 12}, the hardening continues to increase, and here, there is still further decrease in the level of elasticity attainable, confirming earlier observations. Salient-corner vortex cell-size adjustment against We, measured

112

J.P. Aguayo et al. / J. Non-Newtonian Fluid Mech. 134 (2006) 105–126

on the vertical wall (L) and horizontal wall (X), follow similar patterns to each other, see Fig. 5a. We note, the magnitude of L is larger than for X at any chosen We. In addition, note that similarity in growth-trend for X and L can be associated with all the strongly strain-hardening cases considered (q ≥ 5, Fig. 5a). For trends in material functions see Fig. 1a. The patterns are almost identical for elasticity levels up to 10; whilst for the (q = 2) case there is slight departure in X and L-trends. In Fig. 5b, salient-corner vortex intensity (ψsal ) grows continuously when increasing the elasticity level for q = 5, 10 and 15. This is also true for q = 2 with We ≤ 5; beyond We = 5, there is a decline in intensity as some degree of softening takes hold (see arguments below). Plots of growth-trends in X and L are based on the streamlines patterns in Fig. 6 (q = 2, 5) and Fig. 7 (q = 10, 15), which reflect a relatively linear shape of separation-line throughout. This outstanding feature, clear at {q = 2, Wecrit = 60}, cannot be associated with response in shear rheology, as the four fluids (different q-cases Fluid-{qi }) exhibit practically identical shear viscosity. To explain these differences we may look to the strain-rate field, dxx , see on to Fig. 15a–d. One can observe that dxx -maxima occur in the entry-flow region, specifically around the corner and on the centre-plane. Although these fields are practically identical, the extensional viscosity for each value of q is quite different. We observe for (q = 2), that the extensionsoftening regime begins at λ0b ˙ ≈ 0.2, whilst for q = 5, 10 and 15, the softening regime is reached for λ0b ˙ ≥ 9.0. We can discern that for (q = 2), we have a large area of flow where we would anticipate strain-softening. For (q = 5) or more, this area reduces considerably, if not disappears altogether. This is explained in Figs. 8a and b where the regions inside the selected dxx -contour-lines are those where strain-softening is expected from the rheology of each particular fluid. These contour levels capture the approximate strain-rates and above, around which the peak in extensional viscosity occurs. It is important to bear in mind that the flow here is not purely extensional, but of complex dynamical form. Hence, deformation regimes quoted for the on-set of the softening phenomena are simply a guidance when applied to this 4:1 contraction flow. Expected softening regions, in terms of the generalized flow invariant (Γ )1 , are presented in Fig. 8c (q = 2) and Fig. 8d (q = 5). There, it can be observed again that for (q = 2), the softening region is larger than that for (q = 5). Note that such Γ -fields also identify shear deformation as well. We present in Fig. 9 profiles along the plane-of-symmetry of first normal stress-difference (N1 ), stretch (λ), stress component (τyy ) and pressure (p) for We = 10. The contraction is located at x = 27.5. By increasing the number of arms-(q) from 2 to 15, the relative increase from the (q = 2)-level in peak-value of the first normal stress-difference and stretch are around 20% and 10%, respectively. This increase can be associated with the influence of extensional viscosity, where we observe a rise at a strainrate of around 10 units as q increases. The relative increase in upstream pressure at the plane-of-symmetry is about 9% from the reference base of Fluid-{q2 } to Fluid-{q15 }.

1

√ Definition, Γ = 2 I2 , where I2 =

1 2

trace(d 2 ).

Fig. 5. Salient-corner vortex: increasing We; q = 2, 5, 10 and 15: (a) cell-size; (b) intensity.

An alternative representation is to consider pressure correction, available to us in the form of total pressure drop (p) relative to its Newtonian counterpart (pNewt ). Here, (p/pNewt ) decreases significantly below the Newtonian threshold for We < 2, asymptoting to a limit rapidly thereafter (Fig. 10). To quantify this discrepancy between q-solutions, we may comment that at We = 1, (p/pNewt ) with q = 2 falls away by 58% from the Newtonian position, this becoming 56% for q = 15; substantiating a 2% difference. This relative discrepancy is upheld at the larger value of We = 10, with levels dropping further by O(30%). This position may be understood through the shear viscometric properties in that more excessive shear-thinning is experienced by Fluid-{q2 } over Fluid-{q15 }, which leads to larger pressure-drops in the former instance. Profiles of N1 , λ, τxy and τyy along the downstream wall (y = 3) are displayed in Fig. 11. Here, we observe a significant increase in the level of N1 and λ in comparison to plane-ofsymmetry results. The relative increase in N1 from the reference base of (q = 2) to (q = 15) is more than O(100%). Equivalently, λ increases by O(50%). A quantitative measure of mesh convergence is provided through tabulation of the salient-corner vortex cell-size in Table 2. This data is for (q = 5) with solutions across three meshes, covering the range of 0.1 ≤ We ≤ 25. The information in Table 2 demonstrates that convergence with

J.P. Aguayo et al. / J. Non-Newtonian Fluid Mech. 134 (2006) 105–126

113

Fig. 6. Stream function: increasing We; q = 2 and 5.

mesh refinement has been achieved. In addition, vortex enhancement is clearly apparent with increasing Weissenberg number. 6.2. Stretch, stress and deformation-rate The main solution differences that arise when adjusting the number of branched-arms are as follows. Since the components of the deformation tensor have their largest value around the corner, this stimulates corresponding response in all solution

components at the same locations. Also important differences can be anticipated along the centre-plane near the entry-zone. Once the fluid has progressed beyond the contraction-plane, a strong shear flow develops, with large values of stress (τxx ) and stretch (λ) along the downstream wall. Differences in this region are not related to shear viscosity, but with variation in coefficient of the first normal stress-difference (ψ1 ). Elsewhere and where moderate shear deformation is dominant, only minor differences are detected in the solution, taking into account that variation in arms-q has insignificant influence on shear viscos-

Table 2 Mesh convergence: salient-corner vortex cell size (X), 0.1 ≤ We ≤ 25 and q = 5 Meshes

We 0.1

1.0

5.0

10.0

13.0

15.0

18.0

20.0

25.0

(4 : 1)a (4 : 1)b (4 : 1)c

0.151 0.152 0.152

0.162 0.163 0.163

0.182 0.184 0.184

0.190 0.194 0.194

0.200 0.203 0.203

n/a 0.205 0.205

n/a 0.207 0.207

n/a n/a 0.212

n/a n/a 0.218

114

J.P. Aguayo et al. / J. Non-Newtonian Fluid Mech. 134 (2006) 105–126

Fig. 7. Stream function: increasing We; q = 10 and 15.

Fig. 8. Rate-of-strains (dxx and Γ ) fields: We = 10; q = 2 and 5.

J.P. Aguayo et al. / J. Non-Newtonian Fluid Mech. 134 (2006) 105–126

115

Fig. 9. Profiles along plane of symmetry: We = 10; increasing q.

ity. One observes that centre-plane stretch gradually recovers its equilibrium state having transcended the contraction. 6.2.1. Trends in λ For molecular stretch, λ in Fig. 12e–h, we can distinguish regions where there is relatively unstretched material. These regions correspond to inflow and recirculation zones, and along the centre-plane upstream/downstream of the contraction, see Bishko et al. [18]. Such regions of unstretched material correspond to zones where shear dominates over extensional deformation. In addition, there is a ‘banded entry-flow’ zone of stretched material (termed as in [18]), where the influence of increasing levels of extensional viscosity with the number of branched-arms is reflected in larger extension of the polymer

Fig. 10. Normalised pressure-drop vs. We; q = 2, 5, 10 and 15.

molecules. As discussed above, the four fluid cases for q-values investigated (Fluid-{qi }) present similar patterns throughout these zones. Maximum backbone-stretch lies around the corner and also along the downstream wall, see Figs. 11 and 12. For q = 2, stretch generates a value of ≈1.8 at the boundary wall, whilst for q = 15, this rises to ≈2.8. At larger levels of elasticity (We = 25), solution differences with q-increase become more apparent over (We = 10)-equivalents. Fig. 13c and d illustrates the corresponding position at We = 25, in contrast to that at We = 10 of Fig. 12. The relative increase in λ over the contraction zone, from (q = 2) to (q = 5), lies around 25% for We = 25, as opposed to only 10% increase at We = 10. 6.2.2. Trends in N1 , τyy and τxy Fig. 12a–d illustrates field representation for polymeric stress-difference, (τxx − τyy ). Upstream of the contraction, all four plots are essentially identical, with maxima in (τxx − τyy ) just beyond the contraction-plane, lying just above the corner. The presence of more arms in the molecular structure has resulted in a general increase in N1 over the contraction zone and along the downstream wall. Likewise, we note in Fig. 13a and b the relative and more significant increase in N1 , as we elevate We from 10 to 25. There is now O(35%) increase in N1 when comparing peak contraction zone values between (q = 2) to (q = 5)-level, in contrast to O(20%) at We = 10. The banded entry-flow zone for backbone stretch is also associated with the τyy polymeric stress contribution. The maxima is located around the corner, near the contractionplane, whilst the variation of τyy along the centre-plane (Fig. 9 and Fig. 14a–d) is not particularly significant. Along the downstream wall, τyy becomes negative, and for (q = 15), the magnitude decreases slightly when compared against the (q = 2)-case. Nevertheless, the most substantial change in τyy

116

J.P. Aguayo et al. / J. Non-Newtonian Fluid Mech. 134 (2006) 105–126

Fig. 11. Profiles along downstream wall, y = 3: We = 10; increasing q.

Fig. 12. First normal stress-difference (N1 ) and backbone stretch (λ) fields: We = 10; increasing q.

J.P. Aguayo et al. / J. Non-Newtonian Fluid Mech. 134 (2006) 105–126

Fig. 13. First normal stress-difference (N1 ) and backbone stretch (λ) fields: We = 25; q = 2 and 5.

Fig. 14. Normal stress (τyy ) and shear stress (τxy ) fields: We = 10; increasing q.

117

118

J.P. Aguayo et al. / J. Non-Newtonian Fluid Mech. 134 (2006) 105–126

is located around the corner in the banded zone, as described above. The magnitude of τxy along the downstream wall is observed to increase slightly with increment in arms-q. This holds also for the entry-flow zone, where the differences become more noticeable. There, it is possible to also identify a banded zone similar to that exhibited in λ (Fig. 11 and Fig. 14e–h). 6.2.3. Trends in dxx and dxy For dxx and dxy , the response to an increment in arms-q can be considered as insignificant when compared to that equivalent in stress. There is only a slight increase along the centre-plane for dxx near the contraction-plane. In regions where shear dominates extensional deformation, dxy displays similar patterns to those presented for τxy (see Fig. 15). Regions in the flow of high deformation-rate lead to extremes in material functions, which

themselves influence stress and stretch fields. So, large strainrate zones particularly influence (ηe ), with large shear-rate zones affecting (ηs ) and (ψ1 ). Stretch responds to both zones of large strain-rate and shear-rate, but to a greater extent the larger strainrate zones. Entry flow patterns observed in dxx -fields therefore reflect those observed in stretch, whilst those corresponding to dxy begin to take up this relationship beyond the contractionplane and just after the rounding of the corner. Hence, it is the superposition of both these fields that affects λ-fields, which we may observe by back-reference to Fig. 12. The story is somewhat similar but interchanged, when relating N1 , τyy and τxy -fields to stretch, through Figs. 12 and 14. Here, it is τyy - and τxy -fields that actively impact upon λ in the entry-zone, whilst it is N1 (through τxx ) that retains the dependency above the rounded-corner and beyond. There is clear evidence here that molecular stretch responds to the complex deformation that the fluid experiences, as

Fig. 15. Rate-of-strains (dxx ) and (dxy ) fields: We = 10; increasing q.

J.P. Aguayo et al. / J. Non-Newtonian Fluid Mech. 134 (2006) 105–126

119

Fig. 16. Rheological properties of XPP model, variation in α :  = 1/3, β = 1/9 and q = 15.

does stress. The viscometric functions may be taken as a guidance only of anticipated response in such complex dynamical configurations. 6.3. Pom-pom solutions with α-variation {0.15, 0.25, 0.5} Fig. 16 illustrates the adjustment in rheological properties of the XPP model when the values of α are varied. Singularity

in shear viscosity and first normal stress-difference coefficient becomes more pronounced upon α-increase. Moreover, in extensional viscosity, an unrealistic flattening effect is observed. This unexpected response has influence upon both the Trouton ratio and stretch. Solution profiles along the downstream wall are included to particularly detect any abnormal behaviour here (Fig. 17). For α = 0.15 and 0.25, these profiles are similar in trend and magnitude; note also that there is practically no N2 before

Fig. 17. Profiles along downstream wall, y = 3: increasing α; q = 15, We = 10.

120

J.P. Aguayo et al. / J. Non-Newtonian Fluid Mech. 134 (2006) 105–126

Fig. 18. Rheological properties of XPP and EPPT models: variation in q,  and εptt ; β = 1/9, α = 0.15.

or after the contraction-plane. The most noticeable change is for (α = 0.5)-setting, where there is a general decrease in stress and stretch. That is with the exception of that in τzz , which yields an increase in magnitude of N2 . The observed decay in peak-values near the entry region may be related to the flattening influence of anisotropy on extensional viscosity. This is more apparent in N1 and λ (extensional fields) than in shear stress, and naturally of short-time influence due to the short residence time across this zone. In the downstream shear flow, the observed effect of (N1 , τxx , λ)-solution decrease at α = 0.5, may be interpreted on the basis of more exposure to shear-thinning at the larger level of α = 0.5, most prominent in N1 . The deformation-rates encountered under these circumstances are insufficient to enter the shear regime where we might anticipate singularity. On the positive side, for the selected values of α-parameter introduced into our hybrid finite element-volume implementation, we have been able to successfully generate stable solutions for this problem without encountering the theoretical singularities alluded to above.

6.4. Extensional response across models: pom-pom and PTT One aspect of this study has been to compare and contrast the behaviour in complex flows of both pom-pom (kinetic-based) and Phan-Thien/Tanner (PTT network-based) models, seeking to calibrate results in terms of our prior PTT-solutions and vortex enhancement characteristics [2,3]. In this respect, we configure matching peak extensional rheological response for the pompom model against that for the exponential PTT model (EPTT) at two different choice levels of PTT-parameter εptt . This allows us to independently match over extensional viscosity (ηe ) and Trouton ratio (Tr ), covering two scenarios of severe strain-hardening (EPTT (ε = 0.02)) and modest strain-hardening (EPTT (ε = 0.25)) models, see Fig. 18. Hence, and in principal, one can associate flow response with either the extensional properties or otherwise those moderated by the influence of shear. In order to acquire material functions for the SXPP model that approximate the response of those exhibited by the PTT

Table 3 SXPP and corresponding EPTT parameters, β = 1/9 Fluid-(I)a Fluid-(I)b Fluid-(II) Fluid-(III) Fluid-(IV)

Corresponding PTT

q



α

Comment

ηe : EPTT(ε = 0.02) ηe : EPTT(ε = 0.02) Tr : EPTT(ε = 0.02) ηe : EPTT(ε = 0.25) Tr : EPTT(ε = 0.25)

8 8 5 2 2

0.999999 0.999999 0.5 0.333333 0.075

0.15 0.05 0.15 0.15 0.15

Severe-hardening Severe-hardening Strong-hardening Delayed-softening Softening

J.P. Aguayo et al. / J. Non-Newtonian Fluid Mech. 134 (2006) 105–126

models cited above [2,3], we utilise the same viscosity ratio (β). This ensures that shear and extensional viscosity plateaux, at low and high deformations-rates are anchored at the levels set within the earlier work when increasing (q). Crucially, there are two parameters that influence the level of hardening of the pompom model: the ratio of relaxation times () and the number of molecular side-arms (q). An increment in either of these two variables will produce an increase in the extensional viscosity peak-value. For the chosen (β)-ratio, only the relaxation times ratio () has an effect upon shear response. The adjustment in Tr and ηe may be finely improved by varying both parameters, (q) and (). In this fashion, we arrive at four sets of pompom model parameters, yielding four different fluids of Fluid{(I)a ,(II),(III),(IV)}. Distinction in Fluid-(I)a and -(I)b is drawn through α-adjustment, for which we refer to further detail be-

121

low. The parameters for SXPP and corresponding EPTT-variants are displayed in Table 3, with rheological material functions displayed in Fig. 18. Fluid-(I)a represents the case with the largest extensional viscosity that best matches ηe for EPTT (εptt = 0.02) of severe hardening. Although the shear viscosity for Fluid-(I)a is closer to that for the PTT-fluid than for the other pom-pom fluids, its Trouton ratio exceeds that of the PTTfluid. This is because the network model ‘softens’ more rapidly than does the SXPP. As above, Fluid-(II) best approximates Tr for EPTT (εptt = 0.02). Fluid-(III) displays moderate hardening and closely approximates ηe for EPTT (εptt = 0.25). This was the base case in our previous study on pom-pom models [1,28]. Fluid-(IV) matches Tr for the moderate hardening EPTT (εptt = 0.25) model. Note that, for this set of parameters the fluid does not present any hardening and softens at low extension-rates. The possibility of modifying ηe with a non-

Fig. 19. Stream function: variation in q and , increasing We; match to strong-hardening EPTT (ε = 0.02).

122

J.P. Aguayo et al. / J. Non-Newtonian Fluid Mech. 134 (2006) 105–126

substantial change in shear viscosity is an important feature of the pom-pom model that is absent in these PTT network model variants. We can identify particular aspects of the solutions generated in each of the four fluid cases above, delineated against increasing We. This includes trends in adjustment of vortex cell-size (L along upstream-wall and X along contraction-plane wall), vortex intensity and separation-line curvature. In this rounded-corner planar geometry, it is the salient-corner vortex activity that is of primary interest. We note, that for εptt = 0.02, when greater separation-line curvature changes arise with increasing We, that it becomes more difficult to accurately measure the L-factor on an equitable basis. Here, we may either take this graphically when we are dealing with published results, or otherwise, by considering the nearest nodal neighbour to the upstream-wall as the wall reattachment point. We comment upon the severe and

modest strain-hardening instances separately, taking the match over (ηe ) and (Tr ) in turn. 6.4.1. Match on severe strain-hardening: EPTT (ε = 0.02), ηe and Tr Fluid-(I)a provides the best match to solution response against that for EPTT (ε = 0.02). In Fig. 19, one generally observes larger vortices with pom-pom Fluid-(I)a over Fluid-(II). Yet, the increase in vortex intensity (Ψsal ) and vortex cell-size (L and X) is similar in trend for both pom-pom Fluids-(I)a and (II), as the level of elasticity is increased. This is displayed in Fig. 21. Notably, pom-pom Fluid-(II) with less hardening response, reaches a larger We-value (Wecrit = O(17)) when compared to pom-pom Fluid-(I)a (Wecrit = O(9)). The EPTT (ε = 0.02) model (extreme strain-hardening case) exhibits vortex enhancement provoked from around We = 5. Such a trend continues

Fig. 20. Stream function: variation in , increasing We; match to modest-hardening EPTT (ε = 0.25).

J.P. Aguayo et al. / J. Non-Newtonian Fluid Mech. 134 (2006) 105–126

up to Wecrit of 30, see original Fig. 6 in Aboubacar et al. [2]. Here, we note that the XPP model shows similar solution response to the EPTT (ε = 0.02) model for We ≥ 5. In Fig. 19, there is insignificant change in orientation or shape of the vortex separation-line (linear shape) as the level of elasticity increases. In contrast, the separation-line for the EPTT model commences from a linear shape, gradually becoming more concave. At larger elasticity levels of (We = O(20)), the separation-line adopts a convex shape, at which point we observe an extremely large vortex formation. 6.4.2. Match on modest strain-hardening: EPTT (ε = 0.25), ηe and Tr With pom-pom Fluid-(III), smaller vortices are generated than with pom-pom Fluid-(I)a or (II). Ultimately, this setting leads to vortex reduction. Pom-pom Fluid-(IV) displays no apparent vortex enhancement up to the level of (We = 30), beyond which point there is slight growth both in vortex intensity and cell-size, see Figs. 20 and 21. In contrast to Fluid-(I)a and (II), larger We-values of O(60) are achieved with both pompom Fluids-(III) and (IV). For Fluid-(III), we observe a trend similar to that for EPTT (ε = 0.25); that is, vortex intensity increases up to and around (We = 5), with a small reduction

123

thereafter for We ≥ 5. In contrast, vortex enhancement is absent up to (We = 20) for Fluid-(IV). Here, vortex intensity increases slightly when it reaches the limiting state of Wecrit = 60. At the largest We (We = 100), the EPTT (ε = 0.25) model shows significant reduction in vortex intensity and cell-size. This is regarded as a consequence of its strain-softening response and associated levels of Trouton ratio. The extensional viscosity of EPTT (ε = 0.25) is lower than that for pom-pom Fluid-(III) at strain-rates larger than O(two units) (see Fig. 18). This accounts for the higher levels of elasticity reached for the EPTT model of O(100), compared to O(60) for the XPP model. For the delayed-softening case (Fluid-(III)), the separation-line is linear, as observed with Fluid-(I)a and Fluid-(II). In contrast, with Fluid-(IV), the shape of the separation-line is relatively convex. The separation-line for the corresponding EPTT model changes from linear (We = 1) to convex shape (We = 5), and for elasticity levels larger than (We = 20) recovers linear shape once more, see Aboubacar et al. [2]. 6.4.3. Fluid-(I)a and Fluid-(I)b : similarities and differences In order to detect the influence of variation in the selected values of α, we introduce the fine distinction between {Fluid(I)a , α = 0.15} and {Fluid-(I)b , α = 0.05}. To this point, we

Fig. 21. Salient-corner vortex intensity and cell-size: increasing We; variation in q;  and εptt .

124

J.P. Aguayo et al. / J. Non-Newtonian Fluid Mech. 134 (2006) 105–126

Fig. 22. First and second normal stress-difference (N1 , N2 ) fields: α = 0.05 and 0.15; We = 8.

have presented results for our standard value of α = 0.15. Field contour plots for first and second normal stress-difference are provided in Fig. 22 for these two instances, Fluid-(I)a and Fluid(I)b . Here, we observe that for N1 , two zones are clearly affected by such parameterisation. The first zone is a small region above the end of the rounding of the corner, where a reduction in first normal stress-difference (N1 ) can be related to typical response in viscometric flow. The second zone lies in the neighbourhood of the downstream wall, where differences in ψ1 at the two selected α-values utilised here, produce a decrease in N1 from states at α = 0.05 to 0.15. The influence of greater anisotropy (larger α) is expected to be more prominent in N2 . Some difference can be appreciated in the entry-zone where increase in α is reflected in a slight increase in N2 . Notwithstanding this, the most substantial variation in the magnitude of N2 is in shearflow along the downstream wall (increasing approximately from 0.07 to 1.5 units). Curiously and in contrast, α-adjustment has had barely any impact on N2 in the all-important converging entry-flow zone. If anything, there is even a slight decrease in N2 in this zone with increase of α. This is an unexpected consequence and leads us to conclude that α-adjustment must be associated primarily with a shear-flow response. The critical level of elasticity, Wecrit , is nine for α = 0.15 and reduces to eight for α = 0.05. This change, though not particularly significant, may be attributed to the reduction in strain-hardening at α = 0.15 over that prevailing at α = 0.05. 7. Conclusions In this paper, we have presented the influence of the number of dangling branched-arms (q) at each end of the pom-pom molecule on flow characteristic in a 4:1 non-abrupt planar contraction. The response under extensional deformation of the ex-

tended pom-pom model is significantly influenced by such parametric variation. With the selected set of parameters, we have commenced with a fluid of (q = 2), Fluid-{q2 }, dominated by its strain-softening behaviour, yet practically devoid of any hardening. We close by considering a fluid of (q = 15), Fluid-{q15 }, where there is practically no tendency to soften over the pertaining kinematic flow conditions. Differences in solutions derived across variation in the number of arms-q are related principally to properties in extensional viscosity and first normal stressdifference, noting the insignificant influence of arms-q on shear viscosity. As a result, we observe that deformation-rates, appreciated through dxx and dxy -fields, are relatively insensitive to variation in q. The impact of extension-softening may be appreciated under current flow conditions for (q = 2) through a decrease in vortex-intensity. In addition, as seen previously in our earlier work [1,28], the degree of hardening can be associated with the levels of elasticity reached. For moderate levels of increase in extensional viscosity, Wecrit is often larger than that obtained for strongly-hardening fluids. Vortex cell-size follows similar trends in all four fluid cases studied, Fluid-{qi }, covering the variant numbers of arms-q presented. It is particularly notable that for fluids devoid of strain-softening, cell-sizes (X and L) are almost identical for pom-pom fluids with arms-q of 5, 10 and 15. Zones in the flow field where extension dominates shear deformation are anticipated to present significant differences at any fixed We, for alternative configurations in the number of branched-arms. Such zones, detected in dxx , lie about the centreplane near the contraction-plane, and around the corner. We can observe for arms (q = 2), there are two distinct zones within the entry-flow, where we might anticipate strain-softening. This is practically absent once the number of branched-arms exceeds five, so that there is barely any exposure to softening, bar within

J.P. Aguayo et al. / J. Non-Newtonian Fluid Mech. 134 (2006) 105–126

an extremely localised neighbourhood about the rounded-corner. One may observe increase in the magnitude of field variables in stress and stretch, although for τyy , this is not as significant as in other stress components. Molecular stretch itself, is observed to continually increase and double over the entry region, this being heightened particularly from configurations with (q = 2) to (q = 5), remaining fairly static in (q)-change thereafter. In addition, pure shear flow is established in the exit-channel, along the downstream wall where extension-rates essentially vanish. In this shear-flow zone, significant increase in solution variables is provoked, with the largest being in N1 and λ. The magnitude of each variable elevates with increase in (q), change being more rapid from configurations (q = 2) to (q = 5) than thereafter. At any chosen (q)-level and as characterised by maxima in N1 , variable levels are sustained along the downstream wall, see associated profile plots of Fig. 11. Adjustment in solution levels for {N1 , λ} with increase in arms-(q) may be attributed to the shear response in first normal stress-coefficient, ψ1 , at the prevailing shear-rates developed in the exit-channel. Beyond parametric variation for the class of pom-pom (kinetic-based) models alone, we have sought to compare and contrast extensional response and vortex characteristics across model types, from pom-pom to EPTT (network-based phenomenological) models. This has revealed that similar trends in solution response may be extracted in terms of vortex enhancement or inhibition, depending upon the selected choice of extensional rheology forged. Strong strain-hardening and matching on peak-ηe with Fluid-(I) equates closest to the ηe -properties for EPTT (ε = 0.02). This setting stimulates the largest vortex response, but also the lowest level of Wecrit . Here, matching on peak-Tr instead with Fluid-(II), lessens vortex activity but doubles Wecrit . Under the moderate strain-hardening regime with Fluid-(III) and Fluid-(IV), these more modest extensional properties have been observed to ultimately lead to vortex inhibition as Tr declines; so that overall, vortices significantly diminish from their hardening counterparts. In this context, now matching on peak-ηe with Fluid-(III) equates to equivalent EPTT (ε = 0.25) extensional properties and corresponding vortex response. Here, slight vortex enhancement is observed up to We ≈ O(5), prior to subsequent decline. Alternatively, matching on Tr with Fluid-(IV), leads to immediate and sustained vortex reduction. Once the relative contribution of solvent and polymeric components has been established by fixing the β-parameter, the material functions for PTT models can be modified by changing the parameter εptt . Then, as εptt approaches zero, the extensional viscosity approximates Oldroyd-B behaviour. Unfortunately, the shear viscosity is also influenced by this parameter. In contrast, for the XPP models it is possible to practically adjust the degree of hardening without affecting the shear response over a working range of β-values and shear-rates. In the current study, this has been accomplished by fixing  and allowing q to vary. We also point out the influence larger values of anisotropy-parameter (α) may have, discussing the possible on-set of viscometric singularities. In general, we may pass comment that for fluids with equivalent levels of hardening, the PTT-version reaches its peak extensional viscosity value prior to its pom-pom equivalent. At the same time, the PTT-version is delayed in entering its shear-

125

thinning regime in contrast to its pom-pom equivalent, at the same levels of deformation-rate. This has its consequential effect upon Trouton ratios for both models, causing additional discrepancies accordingly. With large vortex generation firmly in mind, the strongly strain-hardening model variants are to be preferred. In combination with the knowledge on molecular stretch, this may provide a way forward in practical applications where one seeks to take advantage of uncoiling and recovering structures in a complex flow, alongside the presence of vortices. Acknowledgements The authors acknowledge financial support for the present work from the Engineering and Physical Sciences Research Council, UK, under grant GR/N38619, and for beneficial collaboration thereupon with Professor T.N. Phillips, INNFM, Cardiff University. One author (JPA) is grateful for support received from the Ministry of Education (SEP), Mexico. References [1] J.P. Aguayo, P.M. Phillips, T.N. Phillips, H.R. Tamaddon-Jahromi B.A. Snigerev, M.F. Webster, The numerical prediction of planar viscoelastic flows using the pom-pom model and high-order finite volume schemes, J. Comput. Phys., Submitted for publication. [2] M. Aboubacar, H. Matallah, H.R. Tamaddon-Jahromi, M.F. Webster, Numerical prediction of extensional flows in contraction geometries: hybrid finite volume/element method, J. Non-Newtonian Fluid Mech. 104 (2002) 125–164. [3] M. Aboubacar, H. Matallah, M.F. Webster, Highly elastic solutions for Oldroyd-B and Phan-Thien/Tanner fluids with a finite volume/element method: planar contraction flows, J. Non-Newtonian Fluid Mech. 103 (2002) 65–103. [4] M.F. Webster, H.R. Tamaddon-Jahromi, M. Aboubacar, Transient viscoelastic flows in planar contractions, J. Non-Newtonian Fluid Mech. 118 (2004) 83–101. [5] M.F. Webster, H.R. Tamaddon-Jahromi, M. Aboubacar, Time-dependent algorithms for viscoelastic flow—finite element/volume schemes, Numer. Meth. Part. Diff. Eqns. 121 (2005) 272–296. [6] M. Aboubacar, H.R. Tamaddon-Jahromi, M.F. Webster, Numerical simulation of contraction flows for Boger fluids using finite volume methods: bridge between finite-volume and finite-element methodology. In: K.J. Bathe (Ed.), Proceedings of the Second MIT Conference on Computational Fluid Solid Mechanics, MIT, USA, Elsevier Science Ltd., 2003, pp.815–818. [7] B. Tremblay, Visualization of the flow of linear low density polyethylene/low density polyethylene blends through sudden contractions, J. NonNewtonian Fluid Mech. 43 (1992) 1–29. [8] H. Matallah, P. Townsend, M.F. Webster, Viscoelastic multi-mode simulations of wire-coating, J. Non-Newtonian Fluid Mech. 90 (2000) 217–241. [9] E. Mitsoulis, M. Schwetz, H. M¨unstedt, Entry flow of ldpe melts in a planar contraction, J. Non-Newtonian Fluid Mech. 111 (2003) 41–61. [10] 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. [11] 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 (2001) 823–843. [12] P. Rubio, M.H. Wagner, LDPE melt rheology and the pom-pom model, J. Non-Newtonian Fluid Mech. 92 (2000) 245–259. [13] 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 (1999) 873–896.

126

J.P. Aguayo et al. / J. Non-Newtonian Fluid Mech. 134 (2006) 105–126

[14] M. Zatloukal, Differential viscoelastic constitutive equations for polymer melts in steady shear and elongational flows, J. Non-Newtonian Fluid Mech. 113 (2003) 209–227. [15] 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 (2002) 187–208. [16] R.I. Tanner, S. Nasseri, Simple constitutive models for linear and branched polymers, J. Non-Newtonian Fluid Mech. 116 (2003) 1–17. [17] 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. [18] 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. [19] N. Clemeur, R.P.G. Rutgers, B. Debbaut, Numerical simulation of abrupt contraction flows using double convected pom-pom model, J. NonNewtonian Fluid Mech. 123 (2004) 105–120. [20] W.M.H. Verbeeten, G.W.M. Peters, F.P.T. Baaijens, Numerical simulations of the planar contraction flow for a polyethylene melt using the XPP model, J. Non-Newtonian Fluid Mech. 117 (2004) 73–84. [21] J. Van Meerveld, Note on thermodynamic consistency of the integral pompom model, J. Non-Newtonian Fluid Mech. 108 (2002) 291–299. [22] I. Sirakov, A. Ainser, M. Haouche, J. Guillet, Three-dimensional numerical simulation of viscoelastic contraction flows using the pom-pom differential constitutive model, J. Non-Newtonian Fluid Mech. 126 (2005) 163– 173. [23] 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. [24] R.J. Blackwell, T.C.B. McLeish, O.G. Harlen, Molecular drag-strain coupling in branched polymer melts, J. Rheol. 44 (2000) 121–136. [25] N. Clemeur, R.P.G. Rutgers, B. Debbaut, On the evaluation of some differential formulation for the pom-pom constitutive model, Rheol. Acta 42 (3) (2003) 217–231. [26] N. Phan-Thien, A non-linear network viscoelastic model, J. Rheol. 22 (1978) 259–283.

[27] N. Phan-Thien, R.I. Tanner, A new constitutive equation derived from network theory, J. Non-Newtonian Fluid Mech. 2 (1977) 353–365. [28] M. Aboubacar, J.P. Aguayo, P.M. Phillips, T.N. Phillips, H.R. TamaddonJahromi, B.A. Snigerev, M.F. Webster, Modelling pom-pom type models with high-order finite volume schemes, J. Non-Newtonian Fluid Mech. 126 (2005) 207–220. [29] H. Matallah, P. Townsend, M.F. Webster, Recovery and stress-splitting schemes for viscoelastic flows, J. Non-Newtonian Fluid Mech. 75 (1998) 139–166. [30] J. Donea, A Taylor–Galerkin method for convective transport problems, Int. J. Num. Meth. Eng. 20 (1984) 101–119. [31] O.C. Zienkiewicz, K. Morgan, J. Peraire, M. Vandati, R. L¨ohner, Finite elements for compressible gas flow and similar systems, Proceedings of the Seventh International Conference on Computational Methods of Applied Science and Engineering, Versailles, France, 1985. [32] J.L. Guermond, L. Quartapelle, On stability and convergence of projection methods based on pressure Poisson equation, Int. J. Num. Meth. Eng. 26 (1998) 1039–1054. [33] D.M. Hawken, H.R. Tamaddon-Jahromi, P. Townsend, M.F. Webster, A Taylor–Galerkin based algorithm for viscous incompressible flow, Int. J. Num. Meth. Fluids 10 (1990) 327–351. [34] J. Van Kan, A second-order accurate pressure-correction scheme for viscous incompressible flow, SIAM J. Sci. Stat. Comput. 7 3 (1986) 870–891. [35] P. Wapperom, M.F. Webster, Simulation for viscoelastic flow by a finite volume/element method, Comput. Meth. Appl. Mech. Eng. 180 (1999) 281–304. [36] P. Wapperom, M.F. Webster, A second-order hybrid finite-element/volume method for viscoelastic flows, J. Non-Newtonian Fluid Mech. 79 (1998) 405–431. [37] E.O.A. Carew, P. Townsend, M.F. Webster, A Taylor–Petrov–Galerkin algorithm for viscoelastic flow, J. Non-Newtonian Fluid Mech. 50 (1993) 253–287. [38] R.A. Keiller, Entry-flow calculations for the Oldroyd-B and FENE equations, J. Non-Newtonian Fluid Mech. 46 (1993) 143–178. [39] B. Purnode, M.J. Crochet, Flows of polymer solutions through contractions. Part I. Flows of polyacrylamide solutions through planar contractions, J. Non-Newtonian Fluid Mech. 65 (1996) 269–289.