Anisotropy parameter restrictions for the eXtended Pom-Pom model

J. Non-Newtonian Fluid Mech. 165 (2010) 1047–1054

Contents lists available at ScienceDirect

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

Anisotropy parameter restrictions for the eXtended Pom-Pom model Michiel G.H.M. Baltussen a , Wilco M.H. Verbeeten c , Arjen C.B. Bogaerds b , Martien A. Hulsen a , Gerrit W.M. Peters a,∗ a

Materials Technology1 , Faculty of Mechanical Engineering, Eindhoven University of Technology, P.O. Box 513, 5600 MB Eindhoven, The Netherlands Materials Technology, Faculty of Biomedical Engineering, Eindhoven University of Technology, P.O. Box 513, 5600 MB Eindhoven, The Netherlands c Departamento de Ingeniería Civil, Área de Mecánica de los Medios Continuos y Teoría de Estructuras, Universidad de Burgos, Avenida Cantabria s/n, E-09001 Burgos, Spain b

a r t i c l e

i n f o

Article history: Received 1 April 2010 Received in revised form 4 May 2010 Accepted 10 May 2010

Keywords: eXtended Pom-Pom model Giesekus model Shear viscosity Extensional viscosity Second normal stress coefficient Restricted anisotropy parameter

a b s t r a c t A significant step forward in modelling polymer melt rheology has been the introduction of the Pom-Pom constitutive model of McLeish and Larson [T.C.B. McLeish, R.G. Larson, Molecular constitutive equations for a class of branched polymers: the Pom-Pom polymer, J. Rheol. 42 (1) (1998) 81–110]. Various modifications of the Pom-Pom model have been published over the years in order to overcome several inconveniences of the original model. Amongst those modified models, the eXtended Pom-Pom (XPP) model of Verbeeten et. al. [W.M.H. Verbeeten, G.W.M. Peters, F.P.T. Baaijens, Differential constitutive equations for polymer melts: the extended Pom-Pom model, J. Rheol. 45 (4) (2001) 823–843] has received quite some attention. However, the XPP model has been criticized for the generation of multiple and unphysical solutions. This paper deals with two issues. First, in the XPP model, anisotropy is implemented in a Giesekus-like manner which is known to result in unphysical solutions for non-linear parameter values ˛ ≥ 0.5. Hence, we put forward the conjecture that a similar limitation holds for the XPP model. In the present paper, the limits for the anisotropy parameter are elaborated on and result to be most restraining at high deformation rates where the backbone tube is oriented and the backbone tube stretch approaches the number of arms q. By restricting the anisotropy parameter to a maximum critical value the XPP model produces only one solution, which is the correct physical rheology. In the second part we show that, contrary to the results published by Inkson and Phillips [N.J. Inkson, T.N. Phillips, Unphysical phenomena associated with the extended Pom-Pom model in steady flow, J. Non-Newton. Fluid 145 (2–3) (2007) 92–101], for the special case where the anisotropy parameter equals zero, only one physically relevant solution exists in unaxial extensional. In addition to this physically relevant solution, also solutions exist in the physically unattainable part of the conformation space. However, the existence of these physically unattainable solutions is not a unique feature of the XPP model but rather general for non-linear differential type rheological equations. © 2010 Elsevier B.V. All rights reserved.

1. Introduction The processing of polymer materials has a large influence on the dimensional, mechanical, and optical properties of the end product. The complex rheological behaviour typically encountered in macromolecular fluids is an important reason for that influence. In order to predict the viscoelastic behaviour of polymer melts, simulation tools have been developed, which need constitutive models that can adequately model the polymer dynamics. A significant step forward in modelling polymer melt rheology has been the introduction of the Pom-Pom constitutive model [7]. This model is able to quantitatively predict the correct nonlinear behaviour in both

∗ Corresponding author. Tel.: +31 402474840; fax: +31 402447355. E-mail address: [email protected] (G.W.M. Peters). 1 http://www.mate.tue.nl/. 0377-0257/$ – see front matter © 2010 Elsevier B.V. All rights reserved. doi:10.1016/j.jnnfm.2010.05.002

shear and extension simultaneously for branched materials, such as low density polyethylene melts. Various modifications of the Pom-Pom model have been published over the years in order to overcome several inconveniences of the original model [2,3,8,9,12]. Amongst those modified models, the eXtended Pom-Pom (XPP) model [12] has received quite some attention. This particular model was successfully implemented in a finite element code and was able to satisfactorily predict the behaviour of a commercial LDPE melt in complex flow geometries in a quantitative manner [13,14]. Several authors have criticized the XPP model for both mathematical defects [3] and unphysical solutions [3,5]. Both seemingly alarming issues have each their own specific, yet simple solution. On the one hand, the mathematical defects will mostly be critical in numerical computations in the vicinity of geometric singularities. These can be circumvented by choosing the double-equation version of the XPP model, referred to as DXPP, as mentioned by

1048

M.G.H.M. Baltussen et al. / J. Non-Newtonian Fluid Mech. 165 (2010) 1047–1054

Verbeeten et al. [12] and Clemeur et al. [3]. On the other hand, the unphysical solutions, such as turning points [3] and bifurcation and multiple solutions [3,5], were shown to be related to the anisotropy parameter ˛. This parameter was introduced to produce a non-vanishing second normal stress difference. By restricting this parameter to a maximum critical value, the correct physical solutions do exist and will be encountered starting from an admissible initial conformation tensor. The restriction still leaves enough freedom to fit the second normal stress difference. Clemeur et al. [3] propose to use the double-equation version of the eXtended Pom-Pom model and setting the anisotropy parameter ˛ = 0. In this way, they avoid the mathematical defects of the single-equation version and by choosing ˛ = 0, bifurcation and multiple solutions are also ommited. This unfortunately comes at the cost of losing the second normal stress difference. They suggest an alternative way of introducing a non-vanishing second normal stress coefficient by combining both an upper- and lower-convected time derivative, similar to Johnson and Segalman [6]. However, such a combination of an upper-convected and lower-convected time derivative does not fit within the thermodynamic framework GENERIC [10], i.e. combinations are thermodynamically not allowed and a purely upper-convected or purely lower-convected time derivative is preferred. Since the anisotropy parameter in the XPP model is Giesekuslike and the anisotropy parameter ˛ in the Giesekus model is restricted [11], it is expected that some restrictions are also present for the ˛-parameter in the XPP model. The objective of the present paper is to indicate the limits for the ˛-parameter of the XPP model in order to avoid unphysical and multiple solutions. 2. Modelling To realistically describe the viscoelastic stresses of polymer fluids over a broad range of deformation rates, a multi-mode approximation of the extra-stress tensor  is defined as  =

M 

Gi (ci − I) .

(1)

i=1

Here M is the total number of different relaxation times, Gi is the shear modulus of the ith relaxation mode, ci is the conformation tensor, and I is the unit tensor. The conformation tensor of the ith relaxation mode is defined as ci = 32i Si ,

(2)

with i the backbone stretch and Si the orientation tensor of the backbone tube. In the remainder of this paper we will restrict ourselves to a single mode description of the constitutive behavior and omit the subscript i. For the eXtended Pom-Pom (XPP) model, time evolution of the conformation tensor follows from ∇

b c + [f (c) − 2˛] c + ˛c2 + (˛ − 1) I = 0, in which the function f (c) is given by



f (c) = 2r e(−1) 1 −

1 



+



(3)



1 ˛ 1 − ˛ − tr(c2 − 2c) . 3 2

(4)

Here b is the relaxation time of the backbone tube orientation, taken from the linear relaxation spectrum. ˛ is the anisotropy parameter that influences the second normal stress difference, r = b /s with s the relaxation time for the tube stretch, while  is a parameter determining the influence of the surrounding polymer chains on the backbone tube stretch and is defined as  = 2/q, where q denotes the amount of arms at the end of a backbone. Since the trace operator acting on the orientation tensor  yields 1 tr(c)/3. by definition, the backbone stretch is defined as  =

Eqs. (1)–(4) give the same XPP model as given in Clemeur et al. [3] and Inkson and Phillips [5]. However, instead of being written in terms of the orientation tensor S and backbone tube stretch  or extra-stress tensor  , it is written in terms of the conformation tensor c. A more appropiate equation for the function f (c), consistent with the thermodynamical framework GENERIC [8,14], reads



f (c) = 2r e(−1) 1 −

1 2



+





1 ˛ 1 − ˛ − tr(c2 − 2c) . 3 2

(5)

Since the introduction of the second normal stress difference is Giesekus-like by means of the anisotropy parameter ˛, and that parameter in the Giesekus model is restricted [11], the evolution equation of the Giesekus conformation tensor is also given for comparison ∇

 c + [1 − 2˛] c + ˛c2 + (˛ − 1) I = 0,

(6)

with  the linear relaxation time. 3. Anisotropy parameter restrictions Our starting point is the limitation on the parameter ˛ in the Giesekus model, 0 ≤ ˛ ≤ (1/2), as suggested by Bird et al. [1] and later by Schleininger and Weinacht [11], studying the linear stability of Couette flow. The restriction ensures that the Giesekus model does not give solutions with a maximum in the shear or elongational stress, leading to unstable, non-physical solutions. The restriction also ensures that the linear term in Eq. (6) ([1 − 2˛]c), is positive. We put forward the conjecture that a similar limitation holds for the XPP model, but we are not able to give a formal prove for this and maybe this is not even possible. However, all numerical experiments that we have performed as well as all relevant data obtained from literature support this conjecture. Also a violation of our proposed limitation on ˛ may lead to the flow condition that no steady state solution exists beyond some finite value of 2 (see Appendix A). With this as a starting point and considering the corresponding term in the XPP model [f (c) − 2˛]c, the restriction we put forward reads [f (c) − 2˛] ≥ 0. This leads to a positive linear term in the XPP model. In this way, the maximum anisotropy parameter ˛ allowed becomes a function of the other material parameters and the conformation tensor c, and thus depends on the stretch and orientation. For low shear and elongational rates, stretch and orientation are limited which implies that f (c) ≈ 1 and the behavior of the XPP model reduces to the behavior of the Giesekus model. However, at high deformation rates, f (c) changes significantly and, for given ratio r and number of arms q, the critical value of ˛ becomes a function of the applied deformation rate. Computationally, this can be confirmed by computing the values of the function g(c) = [f (c) − 2˛] for a relevant series of deformation rates and values of ˛. Contours of g = 0, as shown in Fig. 1, reveal the dependence of the maximum value for ˛ as a function of the applied shear or extensional rate. The figure shows that, for fixed ratio of relaxation times r and increasing number of arms q, the minimum allowable value for ˛ appears at high shear or elongational rate. In addition, it shows that the minimum allowable value for ˛ is different in shear and uniaxial extension and appears at different deformation rates. In order to find an expression for the maximum allowable value of ˛ in general complex flows (denoted by ˛max ) we write the function g(c) in the following way: g(c) =

1  2r e(−1) (2 − ) + 1 − ˛(1 + 34 tr(S2 )) . 2 

(7)

The possible values for c (and thus also of  and tr(S2 )) of the XPP model lie on a surface in (c1 , c2 , c3 )-space, where (c1 , c2 , c3 )

M.G.H.M. Baltussen et al. / J. Non-Newtonian Fluid Mech. 165 (2010) 1047–1054

1049

Fig. 1. Maximum allowable ˛ as a function of shear rate (left) and strain rate (right) for a single-mode XPP model. q = 2, 5, 10, 25, r = 3.

Fig. 2. Approximation of the maximum allowable ˛ (˛ ˆ max ) in (r, q)-space. √ Four √ regions can be observed, I and III where q ≥ 2 + 3, II and IV where q < 2 + 3. In III and IV, ˛ ˆ max is limited by the upper bound 0.5, whereas in I and II the approximation specified in the second argument of Eq. (10) is used. Fig. 4. Maximum admissible ˛ according to Eqs. (10) and (11) as a function of the amount of arms q and simple guidelines ˛ = 0.1/q and ˛ = 0.3/q with r = 1.

are the principal values of c (see Appendix A). In addition, this surface depends on the parameters (r, q, ˛) in the model. The value of ˛max is the maximum value of ˛ for which g(c) ≥ 0 on the surface in (c1 , c2 , c3 )-space. In Appendix A we have shown that for a given stretch  the maximum value of tr(S2 ) is obtained for uniaxial elongational flow. Therefore, it follows from the expression for g(c) in Eq. (7) that the minimum value for ˛max is obtained in uniaxial elongational flow. This means we only have to consider a curve instead of a surface in (c1 , c2 , c3 )-space to find ˛max . Note, that in equilibrium we have  = 1 and tr(S2 ) = (1/3) and from Eq. (7) we find g = 1 − 2˛ for this case. This leads to an upper bound ˛max ≤ (1/2) for any parameter set (r, q).

Even though we can restrict the possible values for  and tr(S2 ) to the curve for uniaxial elongational flow, the expressions are still complicated and finding an analytical solution is difficult. Therefore, we try to find an approximation for ˛max . First we approximate the uniaxial elongational flow curve by tr(S2 ) = 1, overpredicting tr(S2 ). If we denote the maximum allowed value for ˛ using this approximation by ˛∗max , we find ˛∗max < ˛max . By this choice g in Eq. (7) is only a function of  and linear in ˛, from which we easily find the value for ˛ for which g = 0 as a function of : ˛∗ () =

2r e(−1) (2 − ) + 1 , 1 + 34

(8)

Fig. 3. Relative error (˛max − ˛ ˆ max )/˛max . The left graph shows results for the original XPP model while the right graph shows results for the thermodynamically consistent XPP model using the same ˛ ˆ max as for the XPP model. The blocky behavior comes from the low amount of grey levels in the colorbar.

1050

M.G.H.M. Baltussen et al. / J. Non-Newtonian Fluid Mech. 165 (2010) 1047–1054

Fig. 5. Contour lines of the function H on a part of the physically attainable cxx , cyy -plane for (G = 1 [Pa], b = 3 [s], s = 1 [s], q = 10 [–], ˛ = 0 [–]) and ε˙ = 0.01 [1/s](left) and ε˙ = 0.1 [1/s](right).

Fig. 6. Contour lines of H on a part of the full cxx , cyy -plane for ε˙ = 0.1 [1/s](left) and corresponding absolute value of the uniaxial elongational viscosity (right). The dashed line is the negative elongational viscosity. The constitutive data corresponds to the data given in Fig. 5.

and ˛∗max is the minimum of this function over the range  ≥ 1. This function is still too complicated to find an analytical expression for ˛∗max . For large values of  the function ˛∗ () can be approximated by ˛∗ () ≈

2r e(−1) (2 − ) 2 = r e(−1) 3 34

 1 2

−

1 3



,

(9)

and the value for  where the right-hand side is minimal can be easily found by solving a quadratic equation. Based on that, we propose the following approximate value for ˛max , denoted by ˛ ˆ max :



˛ ˆ max = min

ˆ

ˆ 2 − ) ˆ +1 1 2r e(−1) ( , 2 4 ˆ 1 + 3



,

(10)

ˆ defined by with 

ˆ = 

⎧ 1+q ⎨ 2

⎩ 1+q+



for q < 2 + (q − 2)2 − 3 2

for q ≥ 2 +

√ 3 √

(11) 3

ˆ is found by minimizing the right-hand The second case for  side√of Eq. (9). Since this would become a complex value for q < 2 + 3, we use the real part only for that range of q. Note, that we included the upper-bound for ˛max of (1/2) found earlier. We can now distinguish four regions in (r, q)-space for √ the value shown in Fig. 2: I and III where q ≥ 2 + 3, II and IV of ˛ ˆ max which are √ where q < 2 + 3. In III and IV, ˛ ˆ max is limited by the upper bound 0.5, whereas in I and II the approximation specified in the second argument of Eq. (10) is used. The introduction of the approximation tr(S2 ) = 1 leads to an underprediction of ˛max . However, the sign of the error intro-

ˆ is undetermined. To confirm that the approximation duced by  given above, and shown in Fig. 2, is actually a conservative one, i.e. ˛ ˆ max ≤ ˛max , we consider the relevant part of the (r, q)-space and compare ˛ ˆ max with the “exact” ˛max using the numerical algorithm as used for producing Fig. 1. As has been shown before, ˛max can be found by considering uniaxial extensional flows only, which simplifies the problem of finding the ˛ ˆ max significantly. In the left graph of Fig. 3, the relative error (˛max − ˛ ˆ max )/˛max is shown for the original XPP model. It can be observed that this error is always positive in the entire (r, q)-space and thus the approximation is a conservative one. Also, it can be observed that the largest error occurs for small r and small q.2 This, however, is never more than 10% for q ≥ 3 and decreases fast in more relevant regions of the (r, q)-space. For the thermodynamically consistent version of the XPP model, one could argue to follow the same strategy to arrive at an estimate ˛ ˆ max . Unfortunately, this does not result in a conservative approxiation. However, as is shown in the right graph of Fig. 3, the approximation given by (10) and (11) also provides a conservative estimate for ˛max for the thermodynamically consistent XPP model albeit at the cost of accuracy. It should be noted that the estimate ˛ ˆ max provides nothing more than a guideline for the consistent implementation of the XPP models and one may always resolve the full set of equations for uniaxial extensional flow to determine the “exact” value ˛max .

2 Combinations of low r and low q are not very physical. Typically, r → 1 for large q-values and r ≥ 3 for q = 2. See [12–14] for a number of illustrative parameters sets.

M.G.H.M. Baltussen et al. / J. Non-Newtonian Fluid Mech. 165 (2010) 1047–1054

1051

Fig. 7. Contour lines of the corresponding H-function for the linear PTT model and ε˙ = 0.1[1/s](topleft) and corresponding value of the elongational viscosity (topright). The bottom graphs show the absolute value of u for the exponential PTT model (bottomleft) and the Giesekus model (bottomright). For all cases G = 1 and  = 3, only the upper convective derivative is applied and the value of the single remaining non-linear parameter is fixed at 0.1. The dashed lines represent the negative elongational viscosities.

In previous papers [12–14], a different guideline for choosing the anisotropy parameter ˛ was given as ˛ = 0.1/q or ˛ = 0.3/q in case no second normal stress difference or second planar viscosity data is available. That simple rule of thumb will generally work fine. However, especially for high q and low r it is recommended to check the admissibility of the ˛-parameter using Eqs. (10) and (11) in order to avoid possible numerical problems. In general, ˛ = 0.1/q will not violate the maximum allowable value for q < 45. Fig. 4 visually illustrates these statements. Inkson and Phillips [5] show numerous plots in their paper where, starting from different initial values, multiple solutions of the XPP model are found. In all of these plots, also a correct physical solution is found, with exception of the cases where the anisotropy

parameter ˛ is above its maximum allowable value. Table 1 shows the cases represented in their paper, the maximum allowable ˛parameter and if a correct physical solution was found or not. As can be seen from Table 1, most studied cases do encounter a correct physical solution, with the exceptions of all cases with ˛ = 1 and the q = 25, ˛ = 0.1 combination. This last case is also the case where Inkson and Phillips [5] found most problems, amongst others a positive second normal stress coefficient. Therefore, the unphysical solutions found by Inkson and Phillips [5] can be explained by the choice of an inadmissible ˛-parameter. Turning points were found in the steady shear viscosity of the XPP model by Clemeur et al. [3], as illustrated by Fig. 1 in their paper. It is speculated that this effect is also due to an anisotropy

Table 1 Non-linear parameters used to calculate bifurcations for the single-mode XPP model with G = 1, r = 3. ˛ ˆ max is calculated according to Eqs. (10) and (11). q

˛

˛ ˆ max

Correct physical solution

f (c)

Reference

2 2 2 2 10 10 10 10 25 25 25 25 10 20

0.00 0.10 0.15 1.00 0.00 0.03 0.10 1.00 0.00 0.012 0.10 1.00 0.03 unknown

0.5 0.5 0.5 0.5 0.1086 0.1086 0.1086 0.1086 0.0209 0.0209 0.0209 0.0209 0.1086 0.0317

Yes Yes Yes No Yes Yes Yes No Yes Yes No No Yes No

Eq. (4) Eq. (4) Eq. (4) Eq. (4) Eq. (4) Eq. (4) Eq. (4) Eq. (4) Eq. (4) Eq. (4) Eq. (4) Eq. (4) Eq. (5) Eq. (4)

[5] [5] [5] [5] [5] [5] [5] [5] [5] [5] [5] [5] [5] [3]

1052

M.G.H.M. Baltussen et al. / J. Non-Newtonian Fluid Mech. 165 (2010) 1047–1054

parameter ˛ above its limiting value. However, since the value for ˛ is not given in the paper, unfortunately, this can not be confirmed. 4. Multiple solutions Inkson and Phillips [5] also find that in uniaxial extensional flows, even for ˛ = 0, multiple solutions exist for low deformation rates. However, we were unable to reproduce these additional solutions. It is beyond the scope of this article to provide a full mathematical proof of the existence of only a single solution for ˛ = 0 but it can be readily shown that the multiple solutions found by Inkson and Phillips cannot be correct. In terms of the conformation tensor, the steady state Eq. (3) in uniaxial extensional flow for ˛ = 0 reduce to





I(cxx , cyy ) = cxx f (cxx , cyy ) − 2b ε˙ − 1 = 0,

(12)

J(cxx , cyy ) = cyy f (cxx , cyy ) + b ε˙ − 1 = 0,

(13)





with cxx and cyy the conformation in the extensional and the perpendicular directions respectively and f (cxx , cyy ) as defined by Eq. (4). If, depending on a fixed value of the deformation rate, the norm of the residues (I and J) on the cxx , cyy -plane is defined as



I(cxx , cyy )2 + J(cxx , cyy )2 ,

H(cxx , cyy ) =

(14)

solutions of (12) and (13) are found when H(cxx , cyy ) = 0. Numerically the function H is easily constructed on a 2-dimensional grid. For comparison, we take Fig. 10(a) of the paper of Inkson and Phillips [5] and the material parameters that it represents (G = 1 [Pa], b = 3 [s], s = 1 [s], q = 10 [–], ˛ = 0 [–]) as a reference. Fig. 5 now shows the contourlines of H in a physically allowable part of the cxx , cyy -plane (where c is positive definite) for ε˙ = 0.01 [1/s] and ε˙ = 0.1 [1/s]. For the dataset chosen here, the XPP model obviously predicts only one solution which is the one we may expect. Of course, this does not show yet that there may not exist other solutions in other parts of the cxx , cyy -plane. However, for the remainder of the physically attainable part of the conformation  space where tr(c) 3, the function f reduces to f ≈ 2r exp( tr(c)/3) 1. For the limit of small deformation rates (ε˙ 1) a linear combination of (12) and (13) yields

 

2 rtr(c) exp 3





tr(c) 3

= 1,

normal stress coefficient in shear flows as well as a second planar elongational viscosity in a Giesekus-like manner. Similar to the Giesekus model [11], we pose that constraints may exist for the anisotropy parameter. However, these are much more strict for the XPP model than for the Giesekus model. The maximum admissible value for ˛ depends on the applied flow but it is shown that it is most restrictive in uniaxial elongational. An estimate for the anisotropy parameter is given by Eqs. (10) and (11). If a too high value for ˛ is chosen, a solution for the XPP model either does not exist or is not stable and admissible for large values of the deformation rate, similar as for the Giesekus model. Both Inkson and Phillips [5] and Clemeur et al. [3] have encountered such incorrect solutions. Contrary to the results obtained by Inkson and Phillips [5], it is also shown that for the special case where ˛ equals zero (which by our definition should be admissible), only one physically relevant solution is obtained in unaxial extensional flows. In addition to this physically relevant solution, also solutions can exist in the physically unattainable part of the conformation space. It is shown that this is not a characteristic of the XPP model but rather a common feature of non-linear differential type rheological equations. In transient flow problems, it was shown by Hulsen [4] that these non-physical solutions can never occur (for admissible values of the non-linear parameter) provided that the initial conformation is positive definite. The XPP model also conforms to the conditions given by Hulsen [4] and the additional solution as shown in Fig. 6 should not lead to any problems when a time dependent flow problem is considered. However, when the steady state solution of a flow characterized by homogeneous elongational deformation is sought iteratively rather than by time integration (as is customary for the linear stability analyses of these flows) it is possible to converge to the physically incorrect solution instead of to the correct one.

Appendix A. The maximum value of tr(S2 ) for given stretch  in the XPP model In this appendix we will show that for given stretch  the value of tr(S2 ) is is maximal in uniaxial elongational flow, i.e. the highest orientation for given stretch is found in uniaxial elongation. In the XPP model the following tensor function has been defined

(15)

which, since b > s , does not have a solution in the conformation space outside the region shown in Fig. 5. Therefore only a single solution is found in the physically attainable conformation space as opposed to the results given in Inkson and Phillips [5]. Still, as is shown in Fig. 6, at least one additional solution exists for indefinite conformation tensors. This additional solution results in a ‘negative elongational viscosity’ which is shown in the right graph of Fig. 6. The same approach, applied to homogeneous shear flow, did not result in any additional (non-physical) solution at low deformation rates. For uniaxial deformation, the obtained non-physical solutions are not a unique feature of the XPP model. A similar analysis of both the linear and the exponential Phan-Thien Tanner (PTT) model (without mixed derivatives) and the Giesekus model shows that also for these models solutions exist for non-physical values of the conformation tensor (Fig. 7). 5. Conclusions The restrictions concerning the anisotropy parameter ˛ of the eXtended Pom-Pom (XPP) model have been investigated. This parameter is responsible for the introduction of a non-zero second

f (c) = f1 () +

1 ˛ [1 − ˛ − tr(c2 − 2c)], 3 2

(16)

where c is the conformation tensor, the stretch  =

  ⎧ ⎨ 2r e(−1) 1 − 1  f1 () =   ⎩ 2r e(−1) 1 − 1 2



tr(c)/3 and

original model,

(17) thermodynamically consistent model.

Note, that c is a positive definite tensor and therefore  is always defined. The function f (c) can be written as follows: f (c)

1 2 1 = f1 () + 2  1 = f1 () + 2  1 = f1 () + 2  = f1 () +





˛ tr(c2 − 2c) 3   ˛ 2˛ 1 − ˛ − tr(c2 ) + tr(c) 3 3   ˛ 1 − ˛ − tr(c2 ) + 2˛2 3   ˛ 1 − ˛ − tr(c2 ) + 2˛. 3 1−˛−

(18)

M.G.H.M. Baltussen et al. / J. Non-Newtonian Fluid Mech. 165 (2010) 1047–1054

1053

And thus g(c)

= f (c) − 2˛





1 ˛ 1 − ˛ − tr(c2 ) 3 2 1  = f1 () + 2 1 − ˛ − 3˛4 tr(S2 )  = f1 () +

where the orientation tensor S has been defined as c c . = S= tr(c) 32

(19)

(20)

Note, that tr(S) = 1 and that S is a positive definite tensor. Based on our approach of Section 3, we propose the following restriction for existence and stability: g(c) > 0.

(21)

For given material parameters (˛, r, ) and specified stretch  the function g(c) is minimal for tr(S2 ) obtaining the maximum value. Next we will show that for steady flow (in a Lagrangian sense) the maximum value of tr(S2 ) is obtained for uniaxial elongation flows. This means we only need to test uniaxial elongation flows for evaluating the minimum value of g(c). We write the XPP model in the following form ∇

b c = h1 (c)I + h2 (c)c + h3 (c)c2

(22)

with h1 (c) = 1 − ˛ h2 (c) = −f (c) + 2˛ = −g(c) h3 (c) = −˛ Note, that h1 and h2 are constants and h2 (c) (like g(c)) is a function of  and tr(S2 ). Following Hulsen [4] we find that ˙ = h1 I2 + (3h2 + h3 I1 )I3 I˙ 3 = I3 tr(c−1 c)

(23)

where the invariants of c are defined by I1 = tr(c) = 32 I2 =

1 2 9 [I − tr(c2 )] = 4 [1 − tr(S2 )] 2 1 2

I3 = det(c) = 276 det(S)

G(tr(S ), det(S)) = h1 I2 + (3h2 + h3 I1 )I3 = 0

which are both positive. In Fig. A.1 we have plotted contours of tr(S2 ) and det(S). Note, that the contours of tr(S2 ) are concentric circles with the center of the circles at the centroid of the triangle. The contour values increase with the radius. The contours of det(S) are curves with a 120 degrees rotational symmetry and a mirror symmetry with respect to the middle line towards the vertices (the uniaxial elongation lines). The contour values are maximal near the centroid and minimal near the edges of the triangle. If we start at a point on the uniaxial elongation lines and stay on the circles of constant tr(S2 ) the value of det(S) will always decrease. We will use this corollory in the following. The function G in Eq. (27) can be written as follows

(24) G = h4 + h5 I3

(30)

(25) (26)

Note, that I1 > 0, I2 > 0 and I3 > 0 since c is a positive definite tensor. Considering only steady solutions (in the Lagrangian sense) we find from Eq. (23) with I˙ 3 = 0 2

Fig. A.1. Contour lines of tr(S2 ) (thick solid lines) and det(S) (thin solid lines) on an equilateral triangle. Also shown are the lines of uniaxial elongation (dashed lines) with s2 = s3 = (1 − s1 )/2, s1 > (1/3) (cyclic). Note, that on the triangle (1/3) ≤ tr(S2 ) < 1, where the lowest value is in the centroid ((1/3), (1/3), (1/3)) and the highest value in the corners (cannot be reached by the model). Also, on the triangle 0 < det(S) ≤ (1/27), where the lowest value is on the edges (cannot be reached by the model) and the highest value is in the centroid.

(27)

describing a “surface” in solution space for c on which all steady state solutions can be found. Using (c1 , c2 , c3 ), the principal values of c, Eq. (27) is a surface in (c1 , c2 , c3 )-space, which consists of the first octant only. For a given stretch , i.e. c1 + c2 + c3 = 32 , the surface G = 0 is intersected by a plane with normal vector (1, 1, 1) and reduces to a curve. Also the intersection of that plane with the solution space of (c1 , c2 , c3 ) (first octant) is a equilateral triangle. The values of c on that triangle are most conveniently described using the orientation tensor S = c/(32 ) and its principal values (s1 , s2 , s3 ). Since s1 + s2 + s3 = 1, the principal values of S represent the barycentric coordinates on the triangle. For a given stretch , the two remaining invariants are tr(S2 ) and det(S) tr(S2 ) = s12 + s22 + s32

(28)

det(S) = s1 s2 s3

(29)

where h4 = h1 I2 = (1 − ˛)I2 h5 = 3h2 + h3 I1 = −3g(c) − ˛I1 We see that for specified  and tr(S2 ) the coefficients h4 and h5 are constant and that h4 > 0 (assuming ˛ < 1) and h5 < 0. Here we have used that g(c) > 0 because of existence and stability requirements we have imposed (see Section 3). Since h4 > 0 it is required that h5 < 0 otherwise the curve G = 0 does not exist. If the criterium g > 0 is violated it is very well possible that h5 can become positive and no solution will exist for the specified stretch . Note, that g < 0 corresponds to ˛ > (1/2) in the Giesekus model, denoting the range for ˛ where stress maxima in shear are observed. In Fig. A.2 we have plotted a curve G = 0 and a circular contour of tr(S2 ) that starts from the intersection of G = 0 with the uniaxial elongation line. Note, that on the edges of the triangle we have G = h4 > 0. Hence the value of G will be positive on the “outside” and negative on the “inside”. It we ‘walk’ along the circle starting from the uniaxial elongation line the value of det(S) will decrease (and thus also I3 ). Since the coefficient h5 is negative in Eq. (30), the value of G will increase (become positive). This means that the curve G = 0 will be on the inside of the circle. Consequently, the starting point on the uniaxial elongation line has the maximum value for tr(S2 ).

1054

M.G.H.M. Baltussen et al. / J. Non-Newtonian Fluid Mech. 165 (2010) 1047–1054

Fig. A.2. Curve G = 0 (thick solid line) for  = 1.1, ˛ = 0.1, r = 3 and  = 0.2. Contour of tr(S2 ) (thin solid line) starting from the intersection of G = 0 with the uniaxial elongation line (dashed line).

References [1] R.B. Bird, C.F. Curtiss, R.C. Armstrong, O. Hassager, Dynamics of Polymeric Liquids, vol. 2, 2nd ed., Wiley, New York, 1987. [2] R.J. Blackwell, T.C.B. McLeish, O.G. Harlen, Molecular drag-strain coupling in branched polymer melts, J. Rheol. 44 (1) (2000) 121–136.

[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 (3) (2003) 217–231. [4] M.A. Hulsen, A sufficient condition for a positive definite configuration tensor in differential models, J. Non-Newtonian Fluid Mech. 38 (1) (1990) 93–100. [5] N.J. Inkson, T.N. Phillips, Unphysical phenomena associated with the extended Pom-Pom model in steady flow, J. Non-Newtonian Fluid Mech. 145 (2–3) (2007) 92–101. [6] M.W. Johnson, D. Segalman Jr., A model for viscoelastic fluid behavior which allows non-affine deformation, J. Non-Newtonian Fluid Mech. 2 (3) (1977) 255–270. [7] T.C.B. McLeish, R.G. Larson, Molecular constitutive equations for a class of branched polymers: the Pom-Pom polymer, J. Rheol. 42 (1) (1998) 81–110. [8] J. van Meerveld, Note on the thermodynamic consistency of the integral Pom-Pom model, J. Non-Newtonian Fluid Mech. 108 (1–3) (2002) 291–299. [9] H.C. Öttinger, Thermodynamic admissibility of the Pompon model for branched polymers, Rheol. Acta 40 (4) (2001) 317–321. [10] H.C. Öttinger, Beyond Equilibrium Thermodynamics Chapter 2, John Wiley & Sons, Hoboken, New Jersey, 2005. [11] G. Schleiniger, R.J. Weinacht, A remark on the Giesekus viscoelastic fluid, J. Rheol. 35 (6) (1991) 1157–1170. [12] W.M.H. Verbeeten, G.W.M. Peters, F.P.T. Baaijens, Differential constitutive equations for polymer melts: the extended Pom-Pom model, J. Rheol. 45 (4) (2001) 823–843. [13] 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. NonNewtonian Fluid Mech. 108 (1–3) (2002) 301–326. [14] 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. NonNewtonian Fluid Mech. 117 (2–3) (2004) 73–84.