Viscoelastic fluids with pressure-dependent viscosity; exact analytical solutions and their singularities in Poiseuille flows

International Journal of Engineering Science 147 (2020) 103207

Contents lists available at ScienceDirect

International Journal of Engineering Science journal homepage: www.elsevier.com/locate/ijengsci

Viscoelastic fluids with pressure-dependent viscosity; exact analytical solutions and their singularities in Poiseuille flows Kostas D. Housiadas Department of Mathematics, University of the Aegean, Karlovassi, Samos 83200, Greece

a r t i c l e

i n f o

Article history: Received 11 May 2018 Revised 24 April 2019 Accepted 9 December 2019

a b s t r a c t We study the pressure-driven unidirectional flow of a viscoelastic fluid with pressuredependent viscosity in a straight channel and a circular tube. The upper convected Maxwell (UCM) constitutive equation is utilized to model the response of the viscoelastic fluid under deformation. The starting point of the analysis is based on the second-order, symmetric conformation tensor. As such, the final model is slightly different from previous models in terms of the extra-stress tensor as the primary variable. The solution of the governing equations is found analytically for all the dependent variables except for a constant which is introduced by the separation of variables method. A physical singular point of the solution is also revealed explicitly. The singularity is taken into account in order to construct very accurate high-order asymptotic solutions for the unknown constant. Further reprocessing of the asymptotic solutions by means of non-linear techniques which accelerate the convergence of series, result in more accurate analytical expressions even close to the singularity of the solutions. © 2019 Elsevier Ltd. All rights reserved.

1. Introduction This paper is a contribution on the subject of fluids with pressure-dependent viscosity. There is early and ample evidence in the literature that the shear viscosity of fluids increases with pressure (Stokes, 1845; Barus, 1891; 1893; Bridgman, 1948). Viscosity pressure-dependent effects become important in engineering applications where large pressures are applied or developed, e.g. in processing of polymer melts, in lubrication, in granular flows, in porous media and others (Denn, 2008; Fusi, Farina & Rosso, 2015; Rajagopal, 2006; Renardy, 2003). Solving the governing equations in such flows, either numerically or analytically, is more difficult than those with constant viscosity due to the additional non-linearities introduced in the equations. In polymer processing, such as injection moulding or flow through dies and extruders (Denn, 2008; Housiadas, 2015; Rajagopal, Saccomandi & Vergori, 2012; Renardy, 2003; and also references therein), the flowing material is usually a high-viscosity viscoelastic fluid which undergoes large external pressures, forces or torques. Thus, the effect of a pressuredependent viscosity becomes even more important than that for Newtonian or generalized Newtonian fluids. However, viscoelastic fluids exhibit complex behavior under flow deformation posing substantial difficulties in their modeling. Fortunately, nowadays, there is a big variety of constitutive equations which describe the dynamics of viscoelastic fluids (Beris & Edwards, 1994). Note that most viscoelastic models are non-linear even under creeping and isothermal flow conditions; the modeling effects become even more demanding for viscoelastic fluids with pressure-dependent viscosity. We also mention E-mail address: [email protected] https://doi.org/10.1016/j.ijengsci.2019.103207 0020-7225/© 2019 Elsevier Ltd. All rights reserved.

2

K.D. Housiadas / International Journal of Engineering Science 147 (2020) 103207

r*

r* z*

2R*

z*

R*

L*

L*

Circular tube (ξ=1)

Straight channel (ξ=0)

Fig. 1. Flow configurations and coordinate systems.

that experimental data for the dependence of the shear viscosity and the relaxation time of viscoelastic fluids have just appeared in the literature (Reynolds, Thompson & McLeish, 2018). In this paper, we investigate pressure-driven viscoelastic flows in a circular tube and a straight channel. We use the Maxwell constitutive equation to determine the extra-stress tensor which we modify suitably to take into account the effect of a pressure-dependent shear viscosity η = η ( p) where p is the pressure. Note that here, we will be referring to the Lagrange multiplier p as pressure. For the constitutive model that we use, the Lagrange multiplier is not the same as the mean normal stress, and if one defines the mechanical pressure as being the mean normal stress, then p is not the mechanical pressure (Karra, Prusa & Rajagopal, 2011; Rajagopal et al., 2012). For the shear-viscosity η two equations are usually reported in the literature. An exponential type equation, knowing as the Barus formula, η = η0 exp(β0 ( p − pre f )) (Barus, 1893; Huilgol & You, 2006; Hron, Malek & Rajagopal, 2001; Karra et al., 2011) where η0 is the shear viscosity at the reference pressure pref , and β 0 is a dimensional constant which can be determined by fitting experimental data. However, the existence of solutions of the governing equations for various formulas η = η ( p) is questionable. For instance, in the case of the Barus formula, Hron et al. (2001) showed that some unidirectional flows are not possible. Thus, its linearized form is preferred, and used here as well, i.e.

η = η0 (1 + β0 ( p − pre f ))

(1)

We emphasize that the viscoelastic model used here, i.e. the UCM model, is based on the conformation tensor instead of the viscoelastic extra-stress tensor (Housiadas, 2015a, 2015b; Karra et al., 2011; Siginer, Akyildiz & Boutaous, 2018). The constitutive model in terms of the conformation tensor is considered as a more reliable and accurate way of describing the polymer dynamics under flow deformation compared to models based on the polymer extra-stress tensor (Beris & Edwards, 1990, 1994; Wapperom & Hulsen, 1998). Thus, the final governing equations are different, with more non-linear terms. Under the assumption of unidirectional flow, the solution of the final equations is found analytically except from an unknown constant  introduced by the separation of variables method (we will refer to that as the “constant of separation”). This method of solution has been used before for various Poiseuille-type flows of fluids with pressure-dependent viscosity (Akyildiz & Siginer, 2016; Housiadas, 2015; Housiadas & Georgiou, 2018; Kalogirou, Poyiadji & Georgiou, 2011; Siginer et al., 2018). The unknown constant  is determined with the aid of the integral form of the total mass balance both numerically and asymptotically. It is shown that  controls the existence of the analytical solutions as well as that  diverges to infinity as the dimensionless pressure-viscosity coefficient β approaches a critical values β c which is found explicitly in terms of the geometrical aspect ratio a and the Weissenberg number Wi (for the definitions of a, β and Wi see Section 2). Note that in Housiadas (2015) a similar solution procedure had been used to determine the critical point numerically only. Moreover, we exploit the fact that the singular point of the solution has been found analytically and we derive highorder asymptotic solutions. The latter are reprocessed using non-linear techniques that accelerate the convergence of series representations of functions in order to increase their accuracy and extend their domain of validity. The acceleration techniques were investigated recently regarding their consistency, accuracy and efficiency for fundamental flow problems with viscoelastic fluids (Housiadas, 2017). It has been demonstrated that when implemented on high-order asymptotic solutions derived for simple shear and steady uniaxial elongational flows produce very accurate transformed analytical solutions over the whole range of the parameters which appear in the governing equations. These techniques have also been used successfully by Housiadas and Georgiou (2018) for the unidirectional two-dimensional flow of a Newtonian fluid with pressuredependent viscosity in a rectangular duct. 2. The mathematical model and its exact solution We consider the isothermal, steady, pressure-driven flow of an incompressible viscoelastic fluid with constant mass density ρ , in two geometrical configurations. In a straight channel of length L, height 2R and width W (where R < L < < W) and in a circular tube of length L and constant radius R (where R < L). For the channel we use a Cartesian coordinate system and for the tube a cylindrical coordinate system; both are shown in Fig. 1. The velocity vector and the total isotropic pressure associated with the constraint of incompressibility are denoted by v and p, respectively. The shear viscosity η of the fluid is assumed to depend linearly on p, following Eq. (1). As far as the single relaxation time λ of the fluid is concerned,

K.D. Housiadas / International Journal of Engineering Science 147 (2020) 103207

3

recent experimental data indicate that it changes only slightly (perhaps within the experimental error) with respect to p (see Fig. 12 in Reynolds et al., 2018), and therefore, as a first approximation, it can safely be assumed constant, i.e. λ = λ0 . In the absence of external forces, the kinematic constraint on the velocity field due to the incompressibility of the fluid, and the equation of balance of linear momentum are, respectively

∇ · v = 0, ρ (v · ∇ )v = −∇ p + ∇ · τ

(2)

In Eq. (2) τ is the extra-stress tensor which is determined by utilizing a suitable constitutive equation. We use the wellknown non-linear differential viscoelastic model which is the Upper Convected Maxwell (UCM) equation. The latter is given in terms of the symmetric, second order conformation tensor C (Beris & Edwards, 1990, 1994; Wapperom & Hulsen, 1998)

δ C/δt := ∂ C/∂ t + (v · ∇ )C − C · ∇ v − (∇ v )T · C = −(C − I )/λ

(3)

In Eq. (3) I is the unit tensor and δ C/δ t represents the upper convected derivative of C; note that C is already a dimensionless quantity. The extra-stress tensor τ is given by

τ = G (C − I )

(4)

where G is the elastic modulus of the fluid. From Eq. (4) we find C = τ /G + I which we substitute in Eq. (3) and using that δ I/δt = −(∇ v + (∇ v )T ), we get

1 δ 1  τ − ( ∇ v + ( ∇ v )T ) = − τ δt G Gλ

(5)

Using G = η/λ, expanding the upper convected derivative and rearranging, we find

    D ln(η ) Dλ δτ T τ +λ −τ +τ = η ∇ v + (∇ v ) δt Dt Dt

(6)

where D/Dt := ∂ /∂ t + v · ∇ is the material derivative. Compared to the classic UCM model which has been derived for constant viscosity and constant relaxation time, Eq. (6) contains two more terms, namely −λ(D ln(η )/Dt )τ and (Dλ/Dt)τ . At steady state and for constant relaxation time Eq. (6) reduces to

    τ + λ (v · ∇ )τ − τ · ∇ v − (∇ v )T · τ − τ (v · ∇ )ln(η ) = η ∇ v + (∇ v )T

(7)

We look for the simplest possible physical solution of the governing equations. This corresponds to unidirectional flow for which vr = 0 and vz = vz (r ), while the remaining flow variables, i.e. the pressure, the viscosity, and the components of the extra-stress tensor, are allowed to depend on both r and z. We scale the axial coordinate z by L, the radial coordinate r by R, the axial velocity component vz by the average velocity U at the exit plane, the pressure difference p − pre f by (3 + 5ξ )η0UL/R2 , and the components of the extra-stress tensor by η0 U/L, where ξ =0 is used for the planar configuration, and ξ =1 for the axisymmetric one. The average velocity at the exit plane is determined with the aid of the volumetric flow-rate Q˙ . For the planar case U = Q˙ /(2W R ), while for the axisymmetric U = Q˙ /(π R2 ). Noting that hereafter all variables and parameters are dimensionless, the final governing equations are

η = 1 + βp

(8)

  ∂p ∂ τzz ∂ τrz τrz − ( 3 + 5ξ ) +a a + +ξ =0 ∂z ∂z ∂r r − ( 3 + 5ξ )

∂p ∂ τrz + a3 =0 ∂r ∂z

(9)

(10)



2τ d u ∂τ v τ ∂η τzz + W i vz zz − rz z − z zz ∂z a dr η ∂z   η d vz ∂ τrz τrz ∂η τrz + W i − v = ∂z η ∂ z z a dr

 =0

(11)

(12)

Notice that the last term in the parenthesis on the left-hand-side of Eqs. (11) and (12) are the additional non-linear terms which appears in the constitutive equation compared to that used by Housiadas (2015a, 2015b) and Siginer et al. (2018). In Eqs. (8)–(12) three dimensionless numbers appear. The aspect ratio of the tube or channel, a, the viscosity pressuredependent coefficient, β , and the Weissenberg number, Wi, defined as

a≡

R , L

β ≡ ( 3 + 5ξ )

β0 η0 LU R2

,

Wi ≡

λ0U L

(13)

4

K.D. Housiadas / International Journal of Engineering Science 147 (2020) 103207

The system of Eqs. (8)–(12) closes with appropriate auxiliary conditions. Symmetry with respect to the axis of symmetry of the circular tube or the midplane of the channel is applied, no-slip of the fluid is imposed at the wall(s), and the pressure datum is set at the exit plane

v z  ( r = 0 ) = v z ( r = 1 ) = p( r = 1 , z = 1 ) = 0

(14)

Last, the dimensionless mass flow-rate at the exit plane, at z = 1, is unity (recall however that vz = vz (r ))



0

1

(2r )ξ vz dr = 1

(15)

2.1. Analytical solution The analytical solution of the partial differential Eqs. (8)–(12) and the accompanying auxiliary conditions (14) and (15) is found by using separation of variables. For a Newtonian fluid with pressure-dependent viscosity and under the assumption of a unidirectional flow, this procedure is described by Kalogirou et al. (2011) for a straight channel and a circular pipe, and also used by Akyildiz and Siginer (2016) and Housiadas and Georgiou (2018) for a rectangular duct. This method was modified suitably for the viscoelastic case by Housiadas (2015); the same is followed here. 2.1.1. Axisymmetric case For the axisymmetric case we set ξ =1 in Eqs. (9) and (10) and we use separation of variables for η, τ rz and τ zz , i.e. η = H (r )eβ (1−z) , τrz = Trz (r )eβ (1−z) , τzz = Tzz (r )eβ (1−z) where H, Trz and Tzz are unknowns functions, and  is a constant of separation that must be determined as part of the solution. Assumptions for the pressure, p, are not required because it depends linearly on the pressure (see Eq. (8)), and the same holds for the velocity which is a function of r only, i.e. vz = vz (r ). Substituting these expressions in Eqs. (8)–(12) we derive a set of ordinary differential equations which we solve analytically. The final solution is



η=

I0 (ca r ) I0 (ca )

1/μa

eβ (1−z ) ,

τzz =

128ηW i a2 ca2



I1 (ca r ) I0 (ca r )

2

,

τrz = −

8η I1 (ca r ) 8 , vz = ln aca I0 (ca r ) ca2



I0 (ca ) I0 (ca r )



(16)

where I1 , I0 are the modified Bessel functions of first and zero order, respectively, and

μa ≡ 1 + 16W i/(a β ), ca ≡



2

(aβ )2 + 16(W iβ )

(17)

The constant  is calculated with the aid of the total mass balance, Eq. (15). Setting ξ =1 and substituting the velocity profile given in Eq. (16) into Eq. (11) gives



1 ln (I0 ( ca ) ) − 2

1

0

ln(I0 ( ca r ) )rdr =

ca2  16

(18)

Eq. (18) can be solved numerically, or asymptotically for small values of the ca parameter (see in Section 3). 2.1.2. Planar case For the planar case we set ξ =0 in Eqs. (9) and (10) and we follow the same procedure as for the axisymmetric case. The final solution is

 η=

cosh(c p r ) cosh(c p )

1/μ p

eβ (1−z ) , τzz =

3η 18ηW i 3 2 tanh (c p r ), τrz = − tanh(c p r ), vz = ln ac p a2 c2p c2p



cosh(c p ) cosh(c p r )

 (19)

where

μ p ≡ 1 + 6W i/(a β ), c p ≡ 2



(aβ )2 + 6(W iβ )

(20)

and the constant of separation  is calculated with the aid of the total mass balance. Setting ξ =0 and using the velocity profile which is given in Eq. (19) into Eq. (15) we get

ln (cosh(c p ) ) −



1 0

ln (cosh( c p r ) )dr =

c2p  3

(21)

Again, Eq. (21) is solved numerically, and asymptotically for small values of cp . It is worth noting here that the constant of separation  depends only on the dimensionless groups cp and ca , as Eqs. (28) and (21) clearly show. In addition, we notice that if we set−ca in place of ca in Eq. (18), and −c p in place of cp in Eq. (21), the equations remain the same due to the fact that all terms in these equations are even functions of ca or cp .Thus, we conclude that  is an even function of ca and cp for the axisymmetric and planar case, respectively.

K.D. Housiadas / International Journal of Engineering Science 147 (2020) 103207

5

Before proceeding with the solution of Eqs. (18) and (21), we mention that as Wi goes to zero, μp and μa go to unity, while the τ zz component of the viscoelastic extra-stress tensor goes to zero. Thus, the Newtonian solution reported previously by Kalogirou et al. (2011) is recovered. For Wi > 0, solutions (16) and (19) are different from the solutions derived by Housiadas (2015) because the equations solved here contain more non-linear terms. Some qualitative differences on the results are discussed in Section 3. 2.2. The singular points of the solution For similar flow problems with a Newtonian fluid Housiadas & Georgiou, 2018; Kalogirou et al., 2011) and with a viscoelastic fluid (Housiadas, 2015), solved with the same method as used here, it has been observed numerically that  increases monotonically with β and eventually goes to infinity at a critical value, β c . Thus, the solution becomes singular at β = βc . We proceed here by investigating if the analytical solutions ((16) and (19) have singular point(s). 2.2.1. Axisymmetric case Assuming that  → ∞ as ca → ca,c ≡

1





ln

I0 (ca ) I0 (ca r )





→

c→ca,c (→∞ )

(aβc )2 + 16(W iβc ), we observe that

ca,c (1 − r )

(22)

Thus, the velocity profile given in Eq. (16) reduces to

vz

→

ca →ca,c (→∞ )

8(1 − r )/ca,c

(23)

Substituting in the integral form of the total mass balance, Eq. (15), performing the integration and solving for ca, c , gives ca,c = 8/3, and hence, vz → 3(1 − r ) as ca → ca, c . Using the definition for ca, c and solving the resulting quadratic equation with respect to β c , we find two solutions



8W i 8 =− 2 ± 3a a

βc1,2

 1+

3W i a

2 (24)

Note that β c1 is positive while β c2 is negative. Thus, the solution has a physical singularity at β = βc1 and a non-physical one at β = βc2 . We will use the symbol β c to refer to the physical singularity, i.e. β c ≡ β c1 . 2.2.2. Planar case

2 As for the axisymmetric case, we assume that  → ∞ as c p → c p,c ≡ (aβc ) + 6(W iβc ). We observe that

1





ln

cosh(c p ) cosh(c p r )



→

c p →c p,c (→∞ )

c p,c (1 − r )

(25)

Thus, the velocity profile given in Eq. (19) becomes

vz

→

c p →c p,c (→∞ )

3(1 − r )/c p,c

(26)

Substituting in the integral form of the total mass balance, Eq. (15), performing the integration and solving for cp, c gives c p,c = 3/2, and hence vz → 2(1 − r ) as cp → cp, c . Using the definition for cp, c and solving with respect to β c we find two solutions



βc1,2

3W i 3 =− 2 ± 2a a

 1+

2W i a

2 (27)

In Eq. (27) β c1 is positive, β c2 is negative, and we will use the symbol β c to refer to the positive root. The above analysis reveals that as the dimensionless group ca approaches ca, c (cp and cp, c , for the planar case, respectively) the constant of separation diverges to infinity and the solution is lost. In terms of the pressure-viscosity coefficient there are two singular points, only one of which is physical. Consequently, Eq. (16) for the axisymmetric case and Eq. (19) for the planar case are valid in the range 0 ≤ ca < ca, c (0 ≤ cp < cp, c ), or in terms of β , in the range 0 ≤ β < β c . Note that for a Newtonian fluid, i.e. for Wi → 0, β c → 8/(3a) for the circular pipe, and β c → 3/(2a) for the straight channel. These limits have also been reported before by Kalogirou et al. (2011). Expressions (24) and (27) also reveal that β c for the viscoelastic fluid is less than its Newtonian counterpart. Since aβ is a positive quantity, rearrangement of Eqs. (24) and (27) gives that the analytical solution is valid in the window 0 ≤ Wiβ < 4/9 for the axisymmetric configuration, and in the window 0 ≤ Wiβ < 3/8 for the planar one; thus the higher the β parameter the lower the Wi number for which the solution exists.

6

K.D. Housiadas / International Journal of Engineering Science 147 (2020) 103207

Fig. 2. The constant of separation  as a function of c/cc .

3. Results In order to present both the axisymmetric and the planar cases succinctly, hereafter we use the symbols



c=

ca , c p,



ξ =1 , ξ =0

cc =

ca,c = 8/3, c p,c = 3/2,

ξ =1 ξ =0

(28)

3.1. Numerical solution of the integral mass balance We complete the solution for both geometric configurations by determining the constant of separation  as a function of the dimensionless group c. Recalling that  is an even function of c, and in order to exploit the fact that  goes to infinity as c → cc , we express the constant of separation  as follows

=

ˆ  2 1 − (c/cc )

(29)

Then, we solve Eqs. (18) and (21) numerically using a composite Gauss–Legendre quadrature for the integrals, and a ˆ . We impose the strictest possible convergence criteria to compute  ˆ with the maximum Newton iterative scheme for ˆ ) are considered exact. The results are shown in Fig. 2 for (machine) accuracy. Hence, the numerical values for  (and  the planar case (solid line) along with the axisymmetric case (line with dots) as a function of c/cc . Note that the numerical ˆ is finite as c → cc (~1.5 for the axisymmetric case, and ~2.0 for the planar one); this implies simulations showed that  that with Eq. (29) the singularities of the solution have been fully removed. 3.2. Asymptotic solution of the integral mass balance ˆ We proceed by determining  approximately using analytic methods. We assume a regular perturbation scheme for  with the small parameter being c2 , namely

 ≈ n =

1

1+

1 − (c/cc )

2

n 



2k c

2k

(30)

k=1

where 2k are constants that must be determined, and in principle n can go to infinity. Notice that the first term in the asymptotic expression (30) is unity, i.e. 0 = 1, since as β goes to zero, c goes to zero and  goes to unity. Then, the velocity profile reduces to the well-known parabolic profile (Poiseuille solution). Substituting Eq. (30) in Eq. (18) for the axisymmetric case (or in Eq. (21) for the planar one), expanding all quantities suitably, and collecting together the terms with regard to c, allows us to find 2k . Except from 0 we have found four more terms in expression (30), i.e. we have found the solution up to O(c8 )



11c2 c4 c6 1 31c8  ≈ 4 = 1− − − + 192 768 30720 13 271 040 1 − 9c2 /64



(31)

K.D. Housiadas / International Journal of Engineering Science 147 (2020) 103207

7

Fig. 3. The error for the approximate analytical solutions for  as a function of c/cc (a) axisymmetric configuration (cc =8/3), (b) planar configuration (cc =3/2).

Fig. 4. (a) The average pressure drop required to drive the flow, (b) the mean Darcy friction factor. Results are shown for the axisymmetric case (ξ =1) with an aspect ratio α =0.1.

For the planar configuration the solution is



 ≈ 4 =

11 2 41c4 127c6 22927c8 1 1− c − − − 2 45 1575 23625 16372125 1 − 4c /9

 (32)

We see that for both configurations the magnitude of the coefficients of c decreases very fast which is an indication that the asymptotic series converges. For this reason, the expression (30) with five terms, i.e. the approximate analytical solution 4 is expected to be a good approximation to the exact solution. Because the coefficients given in (31) are smaller in magnitude than those given in (32), it is also expected that the asymptotic results must be more accurate in the axisymmetric case. Then, we use Shanks’ non-linear transformation Shanks, 1955) and Padé-type diagonal approximants (Padé, 1892) on the numerator only of expression ((30) and we derive new transformed expressions. For instance, the one-level Shanks’ transform and the diagonal (1,1) Padé-type approximant when applied on the first three terms of the numerator in (30) give the same expression (Housiadas, 2017). For the axisymmetric configuration the transformed solution is



 ≈ 1/P = 1/S =

Similarly, for the planar configuration the result is



 ≈ 1/P = 1/S



1 121c2 1+ , 2 1 − 9c /64 48(c2 − 44 )

0 ≤ c < 8/3

(33)



1 847c2 = 1+ , 1 − 4c 2 /9 9(41c2 − 385 )

0 ≤ c < 3/2

(34)

8

K.D. Housiadas / International Journal of Engineering Science 147 (2020) 103207

If we use the first five terms of the numerator in Eq. (30) and apply the two-level Shanks transform and the Padé-type (2,2) diagonal transformation in terms of c2 , we derive 2/S and 2/P , respectively (Housiadas, 2017). These new transformed solutions are long and are not given here. Last, it is worth noting that the asymptotic solutions, and their transformed expressions, are valid for both Newtonian and UCM fluids because they have been derived in terms of c. In Fig. 3, we compare the percentage absolute relative error, Error%, for some of the approximate solutions app that we have derived, where app = 2 , 4 , 1/S ≡ 1/P , 2/S , and 2/P . The error is based on the numerical value N which, as mentioned above, it is considered exact, i.e. Er ror % = 100|1 − app /N |.The results are presented as a function of c/cc for both geometrical configurations in Fig. 3(a) and (b). The horizontal line corresponds to 5% error which we consider as an acceptable tolerance in the calculations because in this case the differences between the exact (numerical) values for  and the analytical results are not even visible. For the axisymmetric configuration, and for all the analytical approximate solutions, we observe that the error, shown in Fig. 3(a), is extremely small over the whole range of c/cc except from the region very close to the singularity, i.e. as c/cc approaches unity. We also see that 4 (given by Eq. (32)) is consistently more accurate than 2 (given by Eq. (32) simply by neglecting the O(c6 ) and O(c8 ) terms) except again from the region close to the singularity where 4 is slightly worse than 2 . The latter can be used with an error less that 5% up to c/cc ≈ 0.979, while the former can be used with the same error up to c/cc ≈ 0.967. Regarding the error of 1/S ≡ 1/P (given by Eq. (33)), it is 2–4 orders of magnitude lower than the corresponding error for 2 except from the region close to the singularity. It can used with a tolerance less than 5% up to c/cc ≈ 0.95. Last, we mention that the errors of 2/S and 2/P are very similar with the error of 4 , i.e. an improvement is not observed. The results for the planar case are presented in Fig. 3(b). First, we see that the error for the analytical solutions increases consistently with c/cc , and also that 4 is more accurate than 2 over the whole range of c/cc . When 2 is transformed to 1/S ≡ 1/P (given by Eq. (34)) the error reduces but not as much. Last, we see when 4 is transformed to the expressions 2/S and 2/P , 1–3 orders of magnitude reduction of the error is achieved, with 2/S being slightly better than 2/P . The slight superiority of the 2-levels Shanks transformation over the diagonal (2,2) Padé approximant has also been observed before by Housiadas (2017). 3.3. Quantities of interest Interesting quantities that are reported for internal, pressure-driven, flows are the mean Darcy friction factor, Cf , and the average pressure difference,  p , required to drive the flow

Cf = −

a 3+ξ

 p =

 0

1



0

1

τrz (r = 1, z )dz

(35)

( p(r, z = 0 ) − p(r, z = 1 ))(2r )ξ dr

where in general for any function φ , φ := φ (z = 0 ) − φ (z = 1 ) and φ := tion can be derived (Housiadas, 2015)

C f =  p −

a2  τzz

3 + 5ξ

(36) 1 0

(2r )ξ φ dr. Using Eq. (9), the following equa(37)

from where it is seen that Cf and  p are not the same due to the presence of the viscoelastic extra-stress τ zz ; recall that for a Newtonian fluid τzz = 0, and the same holds for a viscoelastic fluid with constant viscosity. Using the analytical solutions, the average pressure drop is found as follows

⎧ 2(eβ − 1 )  1 1 / μa ⎪ ⎪ rdr, ξ = 1 ⎨ 0 (I0 (ca r ) ) 1/μ β I0 (ca ) ) a (  p = 1 (eβ − 1 ) ⎪ ⎪ (cosh(c p r ) )1/μ p rdr, ξ = 0 ⎩ 1/μ p 0 β (cosh(c p ) )

(38)

Similarly, the mean Darcy friction factor is found as follows

⎧ β 2(e − 1 ) I1 (ca ) ⎪ , ξ =1 ⎨ β ca I0 (ca ) Cf = ⎪ (eβ − 1 ) ⎩ tanh(c p ), ξ = 0 β c p

(39)

Note that we have also calculated  τ zz and confirmed that Eq. (37) is satisfied within machine accuracy. Results for  p as a function of Wi β are shown in Fig. 4(a) for the axisymmetric case. Recall that due to the singularity, the analytical solution is valid in the window 0 ≤ Wi β < 4/9. The aspect ratio is a = 0.1, and the pressure-viscosity coefficient is β = 0.004 and 0.1. A substantial monotonic increase of the required pressure drop is observed as the Weissenberg number increases. We also see that the results are not very sensitive to the β value. Regarding Cf , which is shown

K.D. Housiadas / International Journal of Engineering Science 147 (2020) 103207

9

in Fig. 4(b), it is seen that this decreases with the Weissenberg number. This is in qualitative difference with the corresponding results previously derived using the constitutive model based on the polymer extra-stress tensor as a primary variable for which Cf decreases with the increase in Wi (Housiadas, 2015). Experimental validation of these predictions is definitely required to confirm which model describes better the effect of fluid elasticity for a fluid with pressure-dependent viscosity. 4. Conclusions We have studied internal pressure-driven flows of viscoelastic fluids for which the shear viscosity depends linearly on the pressure. The extra-stress tensor is described in terms of the conformation tensor resulting in the constitutive equation with more non-linear terms than the constitutive equation used previously in the literature in terms of the extra-stress tensor as a primary variable. The solution of the governing equations has been found analytically for all the flow variables, except for a constant which has been introduced due to the technique of separation of variables. A physical singular point of the solution, cc , has also been revealed explicitly. This has been used to determine β c and thus the analytical solution is defined and is valid in the range 0 ≤ β < β c . We have also determined that the solution is valid in the window 0 ≤ Wiβ < 4/9 for the axisymmetric case, and in the window 0 ≤ Wiβ < 3/8 for the planar configuration. The singularity, cc , has been taken into account in order to calculate the constant of separation by solving the integral form of the mass balance numerically, and asymptotically for small values of the dimensionless group c. The asymptotic expressions are very good approximations of the exact (calculated numerically) solution over the whole range of c/cc . When these solutions are reprocessed further by means of Shanks’ non-linear transformation and Padé-type diagonal approximants more accurate analytical solutions are derived, even close to the singularity. The analytical results derived here can be used as a reference for finding analytical solutions for more complex problems, or for developing numerical methods for viscoelastic flows with pressure-dependent viscosity. References Akyildiz, F. T., & Siginer, D. (2016). A note on the steady flow of Newtonian fluids with pressure dependent viscosity in a rectangular duct. International Journal of Engineering Science, 104, 1–4. Barus, C. (1891). Note on the dependence of viscosity on pressure and temperature. Proceedings of the American Academy of Arts and Sciences, 27, 13–18. Barus, C. (1893). Isothermals, isopiestics and isometrics relative to viscosity. American Journal Sciences, 45, 87–96. Beris, A. N., & Edwards, B. J. (1990). Poisson bracket formulation of viscoelastic flow equations of differential type: A unified approach. Journal of Rheology, 34(4), 503–538. Beris, A. N., & Edwards, B. J. (1994). Thermodynamics of flowing systems. New York: Oxford University Press. Bridgman, P. W. (1948). Viscosities to 30,0 0 0 kg/cm2 . Proceedings of the American Academy of Arts and Sciences, 22, 1471–1475. Denn, M. M. (2008). Polymer melt processing. New York: Cambridge University Press. Fusi, L., Farina, A., & Rosso, F. (2015). Mathematical models for fluids with pressure-dependent viscosity flowing in porous media. International Journal of Engineering Science, 87, 110–118. Housiadas, K. D. (2015a). An exact analytical solution for viscoelastic fluids with pressure-dependent viscosity. Journal of Non-Newtonian Fluid Mechanics, 223, 147–156. Housiadas, K. D. (2015b). Internal viscoelastic flows for fluids with exponential type pressure-dependent viscosity and relaxation time. Journal of Rheology, 59(3), 769–791. Housiadas, K. D. (2017). Improved convergence based on linear and non-linear transformations at low and high Weissenberg asymptotic analysis. Journal of Non-Newtonian Fluid Mechanics, 247, 1–14. Housiadas, K. D., & Georgiou, G. C. (2018). Analytical solution for the flow of a Newtonian fluid with pressure-dependent viscosity in a rectangular duct. Applied Mathematics and Computation, 322, 123–128. Hron, J., Malek, J., & Rajagopal, K. R. (2001). Simple flows of fluids with pressure-dependent viscosities. Proceedings of the Royal Society of London A, 457, 1603–1622. Huilgol, R. R., & You, Z. (2006). On the importance of the pressure dependence of viscosity in steady non-isothermal shearing flows of compressible and incompressible fluids and in the isothermal fountain flow. Journal Non-Newtonian Fluid Mechanics, 136, 106–117. Kalogirou, A., Poyiadji, S., & Georgiou, G. C. (2011). Incompressible Poiseuille flows of Newtonian liquids with a pressure-dependent viscosity. Journal of Non-Newtonian Fluid Mechanics, 166, 413–419. Karra, S., Prusa, V., & Rajagopal, K. R. (2011). On Maxwell fluids with relaxation time and viscosity depending on the pressure. International Journal of Non-Linear Mechanics, 46, 819–827. Padé, H. (1892). Sur la representation approachee d’une function pour des functions rationeless Thesis. École Normale Supérieure. Rajagopal, K. R. (2006). On implicit constitutive theories for fluids. Journal of Fluid Mechanics, 550, 243–249. Rajagopal, K. R., Saccomandi, G., & Vergori, L. (2012). Flow of fluids with pressure- and shear-dependent viscosity down an inclined plane. Journal of Fluid Mechanics, 706, 173–189. Renardy, M. (2003). Parallel shear flows of fluids with a pressure-dependent viscosity. Journal of Non-Newtonian Fluid Mechanics, 114, 229–236. Reynolds, C., Thompson, R., & McLeish, T. (2018). Pressure and shear rate dependence of the viscosity and stress relaxation of polymer melts. Journal of Rheology, 62(2), 631–642. Shanks, D. (1955). Non-linear transformations of divergent and slowly convergent sequences. Journal of Mathematics and Physics, 34, 1–42. Siginer, D. A., Akyildiz, F. T., & Boutaous, M. (2018). Thermally developing heat transfer with non-linear viscoelastic and Newtonian fluids with pressure dependent viscosity. ASME Journal of Heat Transfer, 140(10), 101701 (7 pages). Stokes, G. G. (1845). On the theories of the internal friction of fluids in motion, and of the equilibrium and motion of elastic solids. Trans. Cambridge Philos. Soc., 8, 287–305. Wapperom, P., & Hulsen, M. A. (1998). Thermodynamics of viscoelastic fluids: The temperature equation. Journal of Rheology, 42(5), 999–1019.