Phys. Fluids 29, 121604 (2017); https://doi.org/10.1063/1.4993782 29, 121604 © 2017 Author(s). Stresses of PTT, Giesekus, and Oldroyd-B fluids in a Newtonian velocity field near the stick-slip singularity Cite as: Phys. Fluids 29, 121604 (2017); https://doi.org/10.1063/1.4993782 Submitted: 30 June 2017 . Accepted: 20 September 2017 . Published Online: 16 October 2017 J. D. Evans, I. L. Palhares Junior , and C. M. Oishi ARTICLES YOU MAY BE INTERESTED IN Normal stress differences from Oldroyd 8-constant framework: Exact analytical solution for large-amplitude oscillatory shear flow Physics of Fluids 29, 121601 (2017); https://doi.org/10.1063/1.4994866 Passive non-linear microrheology for determining extensional viscosity Physics of Fluids 29, 121603 (2017); https://doi.org/10.1063/1.4993736 Helical instability in film blowing process: Analogy to buckling instability Physics of Fluids 29, 121501 (2017); https://doi.org/10.1063/1.4992015 http://oasc12039.247realmedia.com/RealMedia/ads/click_lx.ads/test.int.aip.org/adtest/L16/1772678147/x01/AIP/HA_POF_PDF_AIPPAcademy_2019/HA_POF_PDF_AIPPAcademy_2019.jpg/4239516c6c4676687969774141667441?x https://doi.org/10.1063/1.4993782 https://doi.org/10.1063/1.4993782 https://aip.scitation.org/author/Evans%2C+J+D https://aip.scitation.org/author/Palhares+Junior%2C+I+L http://orcid.org/0000-0002-7046-5108 https://aip.scitation.org/author/Oishi%2C+C+M http://orcid.org/0000-0002-0904-6561 https://doi.org/10.1063/1.4993782 https://aip.scitation.org/action/showCitFormats?type=show&doi=10.1063/1.4993782 http://crossmark.crossref.org/dialog/?doi=10.1063%2F1.4993782&domain=aip.scitation.org&date_stamp=2017-10-16 https://aip.scitation.org/doi/10.1063/1.4994866 https://aip.scitation.org/doi/10.1063/1.4994866 https://doi.org/10.1063/1.4994866 https://aip.scitation.org/doi/10.1063/1.4993736 https://doi.org/10.1063/1.4993736 https://aip.scitation.org/doi/10.1063/1.4992015 https://doi.org/10.1063/1.4992015 PHYSICS OF FLUIDS 29, 121604 (2017) Stresses of PTT, Giesekus, and Oldroyd-B fluids in a Newtonian velocity field near the stick-slip singularity J. D. Evans,1,a) I. L. Palhares Junior,1,2 and C. M. Oishi3 1Department of Mathematical Sciences, University of Bath, Bath BA2 7AY, United Kingdom 2Instituto de Ciências Matemáticas e Computação, Universidade de São Paulo, 13566-590 São Carlos, SP, Brazil 3Departamento de Matemática e Computação, Faculdade de Ciências e Tecnologia, Universidade Estadual Paulista “Júlio de Mesquita Filho”, 19060-900 Presidente Prudente, SP, Brazil (Received 30 June 2017; accepted 20 September 2017; published online 16 October 2017) We characterise the stress singularity of the Oldroyd-B, Phan-Thien–Tanner (PTT), and Giesekus viscoelastic models in steady planar stick-slip flows. For both PTT and Giesekus models in the presence of a solvent viscosity, the asymptotics show that the velocity field is Newtonian dominated near to the singularity at the join of the stick and slip surfaces. Polymer stress boundary layers are present at both the stick and slip surfaces. By integrating along streamlines, we verify the polymer stress behavior of r�4/11 for PTT and r�5/16 for Giesekus, where r is the radial distance from the singularity. These asymptotic results for PTT and Giesekus do not hold in the limit of vanishing quadratic stress terms for Oldroyd-B. However, we can consider the Oldroyd-B model in the fixed kinematics of a prescribed Newtonian velocity field. In contrast to PTT and Giesekus, this is not the correct balance for the momentum equation but does allow insight into the behavior of the Oldroyd-B equations near the singularity. A three-region asymptotic structure is again apparent with now a polymer stress singularity of r�4/5. The high Weissenberg boundary layer equations are found to manifest themselves at the stick surface and are of thickness r3/2. At the slip surface, dominant balance between the upper convected stress and rate-of-strain terms gives a slip boundary layer of thickness r2. The solution of the slip boundary layer shows that the polymer stress is now singular along the slip surface. These results are supported through numerical integration along streamlines of the Oldroyd-B equations in a Newtonian velocity field. The Oldroyd-B model thus extends the point singularity at the join of the stick and slip surfaces to the whole of slip surface. As such, it does not have a physically meaningful solution in a Newtonian velocity field. We would expect a similar stress behavior for this model in the true viscoelastic velocity field. Published by AIP Publishing. https://doi.org/10.1063/1.4993782 I. INTRODUCTION Flows with singularities are particularly challenging for viscoelastic fluids. These may arise due to a non-smooth flow domain or a sudden change in flow conditions. An example of the former is sharp (re-entrant) corners in abrupt contraction flows, whilst a typical situation for the latter is a no-slip surface meeting a shear-free surface in die extrusion. The treatment of such singularities and associated high stress concentrations is one of the most challenging problems in computational non-Newtonian fluid mechanics. Characterising the mathe- matical behavior of a viscoelastic fluid near such singularities is important for the following important reasons: 1. The stress singularity is a test of the rheology: it is not clear that all constitutive equations behave properly at very high stresses. 2. It aids numerical schemes, where an analytical solution in the neighbourhood of the singularity can be used to guide discretisations and improve accuracy. a)E-mail: masjde@maths.bath.ac.uk. Telephone: +44 1225 386994. Fax: +44 1225 386492. Here we focus on a problem involving a sudden change of flow conditions relevant to die extrusion. The limit of large surface tension yields an undeformable free surface, giving the so-called stick-slip problem.1 In such a flow, the fluid enters upstream with the Poiseuille flow and upon exiting the die (which can be planar or cylindrical) transitions to plug flow. It is the sudden change at the die exit from a sticking bound- ary condition within the die to a slipping boundary condition outside it that causes the stress field to be singular. Under- standing this flow problem provides a better insight into the mechanics of extrudate swell.2 Because it is simple in geome- try but numerically difficult due to the singularity, the stick-slip flow problem has been cited as a benchmark problem for viscoelastic flow simulations.3 The viscoelastic fluids that we discuss are the Phan- Thien–Tanner (PTT),4,5 Giesekus,6,7 and Oldroyd-B8 (see also, for example, Ref. 9). The Oldroyd-B model describes well the behavior of polymeric liquids composed of a low con- centration high-molecular-weight polymer in a very viscous Newtonian solvent for moderate shear rates. Termed Boger fluids, they are well-characterised by experimental data (see, for example, Refs. 10 and 11), allowing a comparison of theo- retical predications with experimental results. The Oldroyd-B 1070-6631/2017/29(12)/121604/19/$30.00 29, 121604-1 Published by AIP Publishing. https://doi.org/10.1063/1.4993782 https://doi.org/10.1063/1.4993782 https://doi.org/10.1063/1.4993782 mailto:masjde@maths.bath.ac.uk http://crossmark.crossref.org/dialog/?doi=10.1063/1.4993782&domain=pdf&date_stamp=2017-10-16 121604-2 Evans, Palhares Junior, and Oishi Phys. Fluids 29, 121604 (2017) rheological equation is also referred to as the convected Jef- freys model and can be derived from a molecular model in which the polymer molecule is idealized as a Hookean spring connecting two Brownian beads.9 However, many polymeric fluids exhibit shear thinning, where the viscosity decreases with increasing shear rate as well as non-zero first and second normal stress differences. Two popular models that aim to cap- ture such effects are the PTT and Giesekus models. These models are partly motivated by molecular ideas and modify the Oldroyd-B equations by adding quadratic stress terms. One well known deficiency of the Oldroyd-B model (see Refs. 2, 9, and 12) is that it gives infinite stresses and elongational viscos- ity at finite elongational rates, a feature both PTT and Giesekus correct. However, since the elongational flow is relevant to the stick-slip problem after the fluid exits the die, we will see that the Oldroyd-B model gives similar unphysical behavior of infi- nite stress on the slip surface in contrast to the well-behaved PTT and Giesekus models. The low Reynolds number stick-slip flow of a Newto- nian fluid was originally solved analytically by Richardson1 in the planar case and Trogdon and Joseph13 in the axisym- metric case. Numerical solutions of the same Stokes flow problem have been obtained by Coleman14 using a contour integral formulation and by Georgiou et al.15 using a singu- lar finite element method. The latter uses knowledge of the form of the singularity at the die lip to improve the numeri- cal solution in the neighbourhood of the singularity by using special singular elements. This contrasts markedly with the extremely fine mesh needed for convergence, when Salamon et al.16 used ordinary Lagrangian-type elements. However, the type of singular elements used by Georgiou et al.15 is guided by the singularity of the stress, the local analysis for which was performed originally by Michael,17 and then later by Moffatt18 and Richardson.1 Elliotis et al.19 used a singular function boundary integral method to improve numerical con- vergence results, by incorporating several terms of the singular expansion around the stick-slip location. Due to the practical importance of viscoelastic extrusion, there has been much interest in steady viscoelastic stick- slip flow. Numerical simulations of the Oldroyd-B model are the most common, with PTT also having been considered, but not Giesekus. Early work includes that of Coleman20 who extended the boundary integral method to accommodate the Oldroyd-B model in planar flow and managed to obtain solutions for Weissenberg numbers up to 2.6. Marchal and Crochet21 achieved similar Weissenberg numbers using a mixed finite element method with an over diffusive and incon- sistent streamline upwind method for the same problem. Apelian et al.22 applied a Galerkin finite element method and obtained convergence up to Weissenberg 0.4 for the upper con- vected Maxwell (UCM) model (i.e., the Oldroyd-B model with no solvent viscosity). A questionable stress singularity of r�1 was obtained. For higher Weissenberg numbers, severe stress oscillations occurred near the singularity. Consequently, the authors proposed a modified UCM model that admits a New- tonian stress singularity of r�0.5 and obtained results up to Weissenberg 1.9. Rosenberg and Keunings23 compared sev- eral finite element techniques [Galerkin, streamline-upwind, and streamline-upwind Petrov-Galerkin (SUPG)] as well as a streamline integration scheme for the planar Oldroyd-B prob- lem but in fixed Newtonian kinematics. They only attained Weissenberg numbers of 0.5 and commented upon the pres- ence of steep stress boundary layers at the slip surface. Similar Weissenberg numbers were obtained by Owens and Phillips24 using a spectral domain decomposition method. However, they also investigated smoothing the singularity through regular- ising the sudden change in boundary conditions at the join of the stick and slip surfaces. This regularisation permitted higher Weissenberg numbers to be attained. Salamon et al.25 investigated analytically and numerically the modified situa- tion of allowing partial-slip along the wall. They considered both fixed Newtonian kinematics and the true fully coupled viscoelastic flow field. A logarithmic stress singularity was determined, which compared well with the numerical results of an elastic-viscous stress-splitting (EVSS-G)/SUPG finite element scheme on a highly refined mesh. Fortin et al.26 and Baaijens27 used discontinuous Galerkin finite element methods for planar stick-slip flow of Oldroyd-B and PTT models. Fortin et al. obtained similar results to Rosenberg and Keunings23 for Oldroyd-B in fixed Newtonian kinematics, where mesh convergence could not be obtained near to the slip surface. Relaxing the fixed kinematics, the authors man- aged to obtain Weissenberg numbers up to 4 for Oldroyd-B and no bound for the PTT model with a sufficiently large quadratic stress parameter. Baaijens managed significantly higher Weis- senberg numbers of 25 for UCM (Oldroyd-B with no solvent viscosity) and no limit for PTT. Ngamaramvaranggul and Web- ster28 considered the axisymmetric stick-slip case for Oldroyd- B and obtained convergence for Weissenberg numbers up to 2.2 using a semi-implicit Taylor-Galerkin/pressure-correction method. Karapetsas and Tsamopoulos29 used a mixed finite element EVSS/SUPG scheme for both planar and axisymmet- ric stick-slip flows of PTT fluids. They obtained steady state solutions for high Weissenberg numbers (practically without bound) and carefully extracted the velocity and stress behav- iors at the singularity. Almost all results they give are in the absence of a solvent viscosity. Recent progress has been made on the analytical behavior of the singularity for the PTT and Giesekus models in the work of Evans.30,31 The method of matched asymptotic expansions has been used to show that the solvent stress Ts dominates the polymer stress Tp with the radial behaviors, Ts ∼ r− 1 2 , PTT & Giesekus, Tp ∼   r− 4 11 , PTT, r− 5 16 , Giesekus. (1.1) Boundary layers are required at the stick and slip surfaces, of different character and thicknesses. At the stick surface, the boundary layer thickness is θ ∼   r 1 6 , PTT, r 1 4 , Giesekus, (1.2) whilst at the slip surface, we have boundary layer thicknesses π − θ ∼   r 3 20 , PTT, r 3 14 , Giesekus. (1.3) 121604-3 Evans, Palhares Junior, and Oishi Phys. Fluids 29, 121604 (2017) FIG. 1. Geometry of the stick-slip problem. Here (r, θ) are polar coordinates centred at the singularity as in Fig. 1. Since the solvent stress dominates the polymer stress for these models near the singularity, then the flow field is locally Newtonian as conjectured in Ref. 32 for PTT. As such, we may also investigate the behavior of the Oldroyd-B model in such a flow field and derive the following results: Tp ∼ r− 4 5 , (1.4) with boundary layer thicknesses at the stick and slip surfaces of θ ∼ r 1 2 , π − θ ∼ r, (1.5) respectively. The main results that we present here are as follows: 1. By integrating the PTT and Giesekus models along streamlines, we confirm numerically in Sec. III the asymptotic stress singularities in (1.1), originally derived in Refs. 30 and 31. Further, in the same section, we give an order of magnitude derivation of the bound- ary layer thicknesses (1.2) and (1.3) through balanc- ing components of the polymer stresses and velocity gradients. 2. The corresponding results (1.4) and (1.5) for the Oldroyd- B model in a Newtonian velocity field are also derived in Sec. III, with the details presented in the context of matched asymptotics in Sec. IV. Importantly, the asymp- totics show that the stresses of the Oldroyd-B model are unbounded at the slip layer, in contrast to both PTT and Giesekus models. This conclusion is supported by the numerical results in Sec. III. 3. Finally, in Sec. V, we confirm numerically the core behaviors of the natural stress variables determined for PTT and Giesekus in Refs. 30 and 31, as well as Oldroyd-B in Sec. IV. The analytical results derived here for the Oldroyd-B model in the fixed kinematics of a Newtonian velocity field are consistent with the numerical simulation results of Refs. 23 and 26, where excessively large stresses were obtained near the slip surface. We do not have the asymptotic results for Oldroyd- B when the assumption of a Newtonian velocity field is relaxed, i.e., allowing for a truly viscoelastic flow field. The numeri- cal simulation results for this case in the literature are mixed, with some researchers27 obtaining solutions (at least for mod- erate Weissenberg numbers) and others clearly failing.22 The careful singularity regularisation in the numerical approach in Ref. 24 suggests that this situation may also be badly behaved at the slip surface. Thus at this stage, it is difficult to know if the Oldroyd-B behavior in Newtonian kinematics has any indicative bearing for its behavior in the true viscoelastic flow. However, comparison may also be made with the behavior of the Oldroyd-B model at the re-entrant corner singularity. In the benchmark case of 270◦, Oldroyd-B has a polymer stress singularity of r�0.74 in Newtonian kinematics,33 whilst it is r�2/3 in the viscoelastic velocity field.34 The solvent stresses are less singular in both cases, being r�0.4555 in the Newto- nian kinematics and r�5/9 for viscoelastic kinematics. The wall boundary layers are narrower in the Newtonian kinematics. Taking θ to be the angle measured from the upstream wall, the boundary layer thickness is θ ∼ r0.4555 in the Newtonian kinematics, compared with θ ∼ r0.3333 in the viscoelastic flow field. A similar thickness boundary layer is present at the down- stream wall.35,36 Thus, the fixed Newtonian kinematics gives a more singular stress field and thinner wall boundary layers. It remains to be seen if this same trend in relative behavior for the stresses and boundary layers occurs for Oldroyd-B at the stick-slip singularity. II. PROBLEM FORMULATION We consider the planar stick-slip problem shown sche- matically in Fig. 1. The fluid flows in a channel {(x, y) : x ∈ (−∞,∞), 0 ≤ y ≤ 2H } of half-width H. The channel has a no-slip (stick) surface at y = 0 for x ≥ 0 and a shear-free (slip) surface y = 0 for x < 0, with the centre line y = H being a line of symmetry. The fluid flows from right to left, with an assumed fully developed Poiseuille flow within the channel of mean speed V. On exiting the channel, the flow transitions to a shear- free plug flow of speed V. The fluid develops singular stress fields at the points (0, 0) and (0, 2H) in Fig. 1. Experimentally, a bank of parallel plates can be set up at y = ±2H,±4H , and so on, so that die-swell is suppressed for the fluid exiting the chan- nel. This allows the stick-slip regime to be realised practically for viscoelastic fluids. The governing equations for steady, planar, and incom- pressible creeping flow are taken in the dimensional form ∇ · v = 0, 0 = −∇ p + ∇ · τ, (2.1) where v is the velocity field and p is the pressure. The total extra-stress tensor τ is rheologically split into solvent τs and polymer τp components as follows: τ = τs + τp, (2.2) with the Newtonian solvent stress τs = 2ηs D, (2.3) where ηs is the solvent viscosity and D = 1 2 (∇v + (∇v)T ) is the rate-of-strain tensor. The main constitutive equation for the 121604-4 Evans, Palhares Junior, and Oishi Phys. Fluids 29, 121604 (2017) polymer stress of interest will be the Oldroyd-B model in the form τp + λp 5 τp = 2ηp D, (2.4) with λp a constant polymer relaxation time, ηp a constant polymer viscosity, and the upper-convected stress derivative 5 τp = (v · ∇) τp − (∇v) τp − τp (∇v)T . (2.5) We will compare results to the quadratic stress models of PTT and Giesekus, where (2.4) is modified to τp + λp ( 5 τp + g(τp) ) = 2ηp D, (2.6) with g(τp) =   ε ηp tr(τp)τp, PTT, α ηp (τp)2, Giesekus. (2.7) Both these models introduce an additional parameter that con- trols the influence of the quadratic extra-stress terms, being designated the model parameter ε for PTT and the mobility fac- torα for Giesekus. The model parameter ε is taken positive and usually considered for small values, typically up to 0.05.4,29 The mobility factor can be taken in the range 0 ≤ α ≤ 1, although it is usually restricted to 0 ≤ α ≤ 0.5 (see Refs. 9 and 37). Both models reduce to Oldroyd-B, when these param- eters are set to zero. Boundary conditions of no-slip are taken on the channel walls and no shear-stress on the free surface. Using H and V as the characteristic length and flow speeds, we non-dimensionalise as follows: x = H x̄, v = V v̄, p = η V H p̄, τ = η V H T̄, τs = ηsV H T̄s, τp = ηpV H T̄p, where η = ηs + ηp is the total viscosity. Dropping the overbars, the governing equations become ∇ · v = 0, 0 = −∇ p + ∇ · T, T = βTs + (1 − β)Tp, (2.8) Ts = 2D, Tp + Wi ( 5 Tp + κg(Tp) ) = 2D, (2.9) where now g(Tp) =   tr(Tp)Tp, PTT, (Tp)2, Giesekus, (2.10) with dimensionless parameters of the Weissenberg number Wi, retardation parameter β (dimensionless retardation time or solvent viscosity fraction), and model parameter κ, given by Wi = λpV H , β = ηs η , κ =   ε , PTT, α, Giesekus. We examine these equations near the join of the stick and slip surfaces, where the change in the boundary conditions from no-slip to no-shear stress gives rise to a stress singularity. We take Cartesian and polar coordinates centred at the join of the stick and slip surfaces as shown in Fig. 1. The momentum equation may be written as 0 = −∇p + β∇2v + (1 − β)∇ · Tp (2.11) and since we are interested in the flow behavior near to the singularity, we consider this equation for small radial distances r � 1. Away from both the stick and slip surfaces, we postulate that the solvent stress dominates the polymer stress, i.e., (1 − β)Tp � βTs as r → 0. (2.12) Consequently, (2.11) reduces to the Stokes flow equation 0 = −∇p + β∇2v as r → 0. (2.13) A discussion on the separable self-similar solutions for Stokes flow is given in the work of Richardson,1 noting the earlier work of Michael17 and Moffatt.18 Introducing the stream func- tion ψ, the physically relevant dominant self-similar solution is ψ ∼ 2C0r 3 2 sin ( θ 2 ) sin θ, p ∼ 2βC0r− 1 2 sin ( θ 2 ) , as r → 0, (2.14) with the arbitrary constant C0 being determined by flow away from the singularity. In the Newtonian case, C0 = −(3/2π)1/2 ≈ −0.691,29,32 whereas its value has yet to be determined for the viscoelastic fluids of PTT and Giesekus. The negative value of C0 ensures that the flow is from right to left and exits the channel in Fig. 1. Since we are working on small length scales, the following scalings r =Wi2r̂, v =Wiv̂, Tp = T̂ p Wi , (2.15) remove the Weissenberg number from our system of equations. Without loss of generality, we thus set Wi = 1. The analysis of Secs. III–V may be thought of as being for the hat variables in (2.15), with the Weissenberg number easily being reintroduced using the transformation (2.15) to obtain the original variables. III. STRESS INTEGRATION ALONG STREAMLINES In this section, we investigate the polymer constitutive equation (2.9) in the given Newtonian velocity field (2.14). We will use polar coordinates with the stream function (2.14) for the Stokes problem, written in the form ψ = C0r 3 2 f (θ), f (θ) = 2 sin ( θ 2 ) sin θ. (3.1) The components of the velocity are given by vr = 1 r ∂ψ ∂θ = C0r 1 2 f ′(θ), vθ = − ∂ψ ∂r = −C0 3 2 r 1 2 f (θ). (3.2) In polar coordinates, (2.9) takes the form (v · ∇)Tp rr − 2 ∂vr ∂r Tp rr − 2 r ∂vr ∂θ Tp rθ + κgrr + Tp rr = 2 ∂vr ∂r , (v · ∇)Tp rθ + vθ r Tp rr − 1 r ∂vr ∂θ Tp θθ − ∂vθ ∂r Tp rr + κgrθ + Tp rθ = ( 1 r ∂vr ∂θ − vθ r + ∂vθ ∂r ) , (v · ∇)Tp θθ + 2 vθ r Tp rθ − 2 ∂vθ ∂r Tp rθ − 2 r ∂vθ ∂θ Tp θθ − 2 vr r Tp θθ + κgθθ + Tp θθ = 2 (1 r ∂vθ ∂θ + vr r ) , (3.3) 121604-5 Evans, Palhares Junior, and Oishi Phys. Fluids 29, 121604 (2017) where v · ∇ = vr ∂ ∂r + vθ r ∂ ∂θ and grr =   (Tp rr + Tp θθ )Tp rr , PTT, (Tp rr)2 + (Trθ )2, Giesekus, grθ = Tp rθ (Tp rr + Tp θθ ), gθθ =   (Tp rr + Tp θθ )Tp θθ , PTT, (Tp rθ )2 + (Tp θθ )2, Giesekus. (3.4) Since we consider the velocities as given, (3.3) is a system of ordinary differential equations (ODEs) along streamlines. The streamlines are the level curves of the function ψ = C0r 3 2 f (θ) and parameterising each streamline by θ, we have r = ( ψ C0 f (θ) ) 2 3 . (3.5) Since dr/(rdθ) = vr/vθ holds along streamlines, then v · ∇ = vθ r d dθ , (3.6) where d/dθ is the total derivative with respect to θ along the streamline. A. Asymptotic analysis A certain amount of analytical progress can be made for the streamlines that pass close to the singularity. Upstream of the singularity, the flow is weak with the microstructure little deformed. The streamlines are geometrically self-similar and tangential to the stick surface of the channel wall. The no- slip boundary condition means that the flow is predominately simple shear with viscometric polymer stress behavior. This behavior is confined to a narrow region at the wall, after which radial flow dominates. The thickness of the wall boundary layer will be determined by the analysis below and differs for each viscoelastic model. After leaving the viscometric region, the radial stretching becomes important, where the polymer advects and deforms affinely with the fluid. This occurs in a small region of high velocity gradients where the polymer is dominated by a commonly termed stretching solution. The sol- vent stresses genuinely dominate the polymer stresses for the PTT and Giesekus models in this region. For these two mod- els, this behavior persists until the polymer stresses become comparable to the solvent stresses, which occurs near to the slip surface as the flow becomes elongational. This gives rise to a boundary layer at the slip surface in which the polymer stress growth is arrested for the PTT and Giesekus models. We will see that no such restriction in the growth of the polymer stress occurs for the Oldroyd-B model. For the discussion below, two components of the velocity gradient are important, namely, the shear rate γ̇ and radial strain rate ė given by γ̇ = 1 r ∂vr ∂θ − vθ r , ė = ∂vr ∂r . (3.7) 1. Behavior near the stick surface We begin by discussing the behavior of the solution for small values of θ. At leading order in θ, the stream function and velocity are ψ ∼ C0r 3 2 θ2, vr ∼ 2C0r 1 2 θ, vθ ∼ − 3 2 C0r 1 2 θ2, (3.8) with the shear and radial strain rates γ̇ ∼ 2C0r− 1 2 , ė ∼ C0r− 1 2 θ. (3.9) Since θ � 1, the radial flow dominates the angular flow with also the shear rate dominating the radial strain rate. Using in (3.3) and retaining only the leading-order terms in θ gives 3 2 θ2 dTp rr dθ + 2θTp rr + 4Tp rθ − r 1 2 C0 ( κgrr + Tp rr ) = −2θ, 3 2 θ2 dTp rθ dθ + 3 4 θ2Tp rr + 2Tp θθ − r 1 2 C0 ( κgrθ + Tp rθ ) = −2, 3 2 θ2 dTp θθ dθ + 3 2 θ2Tp rθ − 2θTp θθ − r 1 2 C0 ( κgθθ + Tp θθ ) = 2θ. (3.10) Following the work of Renardy33 for the UCM model, it is convenient to rescale the stresses as Tp rr = Qrr , Tp rθ = θQrθ , and Tp θθ = θ 2Qθθ . Then (3.10) gives 3 2 θ2 *... , Q′rr Q′rθ Q′θθ +/// - + θ *... , 2 4 0 3 4 3 2 2 0 3 2 1 +/// - *... , Qrr Qrθ Qθθ +/// - − r 1 2 C0 *...... , κQ2 rr + Qrr κQrrQrθ + Qrθ κ   QrrQθθ Q2 rθ + Qθθ +////// - = 1 θ *... , 0 −2 2 +/// - , (3.11) ignoring lower order terms in θ in the quadratic stress and rate- of-strain terms. The behavior of (3.11) for small θ is dominated by the last two terms as long as θ � r 1 2 . As soon as θ � r 1 2 , then the first two terms in (3.11) are the most important. These first two terms have solutions proportional to θ2α/3, where α = 1,− 3 2 ,−4 are the eigenvalues of the matrix in the second term. Consequently, as θ increases, the solution grows accord- ing to the positive eigenvalue until it dominates the behavior. We now use these insights to analyse viscometric behavior and estimate when it breaks down. Returning to Eqs. (3.10), for θ � r 1 2 , viscometric polymer stresses satisfy κgrr + Tp rr = 2γ̇Tp rθ , κgrθ + Tp rθ = γ̇ ( Tp θθ + 1 ) , κgθθ + Tp θθ = 0, (3.12) where we have introduced the wall shear rate from (3.9). For the Oldroyd-B model, (3.12) with κ = 0 give the usual relationships Tp θθ = 0, Tp rθ = γ̇, Tp rr = 2γ̇2. (3.13) For the PTT model, we have Tp θθ = 0, Tp rθ + 2κ(Tp rθ )3 = γ̇, Tp rr = 2(Tp rθ )2, (3.14) 121604-6 Evans, Palhares Junior, and Oishi Phys. Fluids 29, 121604 (2017) whilst for the Giesekus model, the equations are κ ( (Tp rr)2 + (Tp rθ )2 ) + Tp rr = 2γ̇Tp rθ , κ ( Tp rr + Tp θθ ) Tp rθ + Tp rθ = γ̇ ( Tp θθ + 1 ) , κ ( (Tp rθ )2 + (Tp θθ )2 ) + Tp θθ = 0. (3.15) Although less simplification takes place in the Giesekus case, following Ref. 38, it is possible to obtain the expressions Tp θθ = 1 2κ ( −1 + φ(Tp rθ ) ) , Tp rr = γ̇(1 + Tp θθ ) κTp rθ − (1 + κTp θθ ) κ , γ̇ = 2κTp rθ ( 1 + (2κ − 1)φ(Tp rθ ) ) ( 2κ − 1 + φ(Tp rθ ) )2 , (3.16) where φ(Tp rθ ) = √ 1 − 4κ2(Tp rθ )2. The first and second expres- sions follow from the third and second equations in (3.15), whilst the third expression gives the shear rate in terms of the shear stress and follows from the first equation in (3.15). The sign choices of the terms involving φ(Tp rθ ) have been chosen so that the expressions reduce to those for the UCM/Oldroyd-B model when κ = 0. The above viscometric relationships for both PTT and Giesekus simplify for small radial distances r � 1, where the shear rate γ̇ is large. Consequently for large shear rates, the viscometric stresses for PTT from (3.14) become Tp θθ = 0, Tp rθ = ( γ̇ 2κ ) 1 3 , Tp rr = 2 ( γ̇ 2κ ) 2 3 , (3.17) whilst for Giesekus, we have from (3.15) or (3.16) that Tp θθ = −1 + ( 2κ γ̇ ) 1 2 ( 1 − κ κ ) 3 4 , Tp rθ = ( 1 − κ κ ) 1 2 , Tp rr = ( 2γ̇ κ ) 1 2 ( 1 − κ κ ) 1 4 . (3.18) We now determine the extent of the region in which the approximate equations (3.10) hold. This viscometric region breaks down when θ ≥ r 1 2 . The discussion in terms of the Q variables indicates that this occurs when the terms of the upper convective derivative become comparable in size. Since the largest stress components are the normal radial stress Tp rr and shear stress Tp rθ , this first occurs when the second and third terms in the first equation of (3.10) become the same size giving the relationship ėTp rr ∼ γ̇Tp rθ , (3.19) on using the shear and radial velocity gradient components from (3.9). Using the viscometric relations (3.13), (3.17), and (3.18), we thus have at viscometric breakdown the fol- lowing relationships between the radial and shear strain rates for the each of the models, Oldroyd-B : ė ∼ 1 2 , PTT : ė ∼ κ 1 3 ( γ̇ 2 ) 2 3 , Giesekus : ė ∼ κ 1 4 (1 − κ) 1 4 ( γ̇ 2 ) 1 2 . (3.20) In terms of polar coordinates, using (3.9), these give the esti- mates of the boundary thicknesses at the stick surface to be θ ∼ rn, n =   1 2 , Oldroyd-B, 1 6 , PTT, 1 4 , Giesekus. (3.21) In Cartesian coordinates, we have y ∼ x1+n using θ ∼ y/x and r ∼ x, which agree with the thicknesses for the PTT and Giesekus stick layers in Refs. 30 and 31 derived through matched asymptotics and dominant balance. Later in Sec. IV, we will use the matched asymptotics approach to give an independent derivation of this result for Oldroyd-B. 2. Behavior in the core region On small radial distances r � 1 near the singularity, but away from both surfaces, we have 1 � Tp � ∇v as r → 0. Consequently the polymer constitutive equation reduces to 5 Tp = 0, which has the exact so-called stretching solution Tp = λ(ψ)vvT as r → 0, (3.22) where the function λ(ψ) is constant along streamlines. We des- ignate the region in which this solution holds as the core region, following the original terminology introduced in Ref. 34 for the analysis at the re-entrant corner. To find the variation of the polymer stress across stream- lines, as described by λ(ψ), we may use the stress scaling in the viscometric region. On a given streamline, the angle θ is related to the radial distance by θ = ( ψ C0 ) 1 2 r− 3 4 , (3.23) whilst the viscometric approximation breaks down at θ ∼ rn, with n as given in (3.21) for the different models. Consequently, r ∼ ( ψ C0 ) 2 3+4n , (3.24) from which we can deduce vr ∼ ( ψ C0 ) 1+2n 3+4n , γ̇ ∼ ( ψ C0 )− 1 3+4n . (3.25) The radial component of the polymer stress from the vis- cometric approximations (3.13), (3.17), and (3.18) is given by Tp rr ∼ γ̇ m, m =   2, Oldroyd-B, 2 3 , PTT, 1 2 , Giesekus. (3.26) 121604-7 Evans, Palhares Junior, and Oishi Phys. Fluids 29, 121604 (2017) At the breakdown, we have Tp rr = λ(ψ)v2 r , which gives λ(ψ) = C1 C2 0 ( ψ C0 )n1 , where n1 =   − 6 5 , Oldroyd-B, − 10 11 , PTT, − 7 8 , Giesekus, (3.27) with C1 an arbitrary constant. Consequently, in component form, (3.22) gives Tp rr = C1r1+ 3n1 2 f ′(θ)2f (θ)n1 , Tp rθ = − 3 2 C1r1+ 3n1 2 f ′(θ) f (θ)n1+1, Tp θθ = 9 4 C1r1+ 3n1 2 f (θ)n1+2, (3.28) which yields the singular polymer stress behaviors as Tp ∼   r− 4 5 , Oldroyd-B, r− 4 11 , PTT, r− 5 16 , Giesekus. (3.29) This gives an alternative derivation of the results for PTT and Giesekus in Refs. 30 and 31. We note that the solvent stress has singular behavior Ts ∼ r− 1 2 , (3.30) on using (3.1), so that the main assumption (2.12) is sat- isfied for PTT and Giesekus but not Oldroyd-B. As such, we are fixing the kinematics for Oldroyd-B, the true veloc- ity field being non-Newtonian (i.e., viscoelastic) near the singularity. 3. Behavior near the slip surface Near the slip surface, θ ≈ π, the stream function and velocity components are ψ ∼ 2C0r 3 2 (π − θ), vr ∼ −2C0r 1 2 , vθ ∼ −3C0r 1 2 (π − θ). (3.31) The radial strain rate is now ė ∼ −C0r− 1 2 , (3.32) with the shear-rate being significantly smaller. Elongational flow thus dominates. Using these expressions in (3.3) and retaining only the leading-order terms in π � θ give 2r dTp rr dr − 2Tp rr − r 1 2 C0 ( κgrr + Tp rr ) = 2, 2r dTp rθ dr + 3 2 (π − θ)Tp rr − r 1 2 C0 ( κgrθ + Tp rθ ) = − 3 2 (π − θ), 2r dTp θθ dr + 3(π − θ)Tp rθ + 2Tp θθ − r 1 2 C0 ( κgθθ + Tp θθ ) = −2, (3.33) where we choose here to use the radial distance as the parameter along streamlines. These equations receive stress values from the core region, which have the streamline behaviors Tp rr = λv 2 r = 4C1 ( ψ C0 )n1 r, Tp rθ = λvrvθ = 3C1 ( ψ C0 )n1+1 r− 1 2 , Tp θθ = λv 2 θ = 9 4 C1 ( ψ C0 )n1+2 r−2, (3.34) on entering this slip surface region. Consequently, we see that the normal radial stress component Tp rr is the largest and growing linearly with the radial distance. The radial strain- rate is decreasing and can become comparable to the normal radial stress component. For the PTT and Giesekus models, the first equation in (3.33) suggests that this first occurs when κTp rr 2 ∼ ėTp rr , which simplifies to κTp rr ∼ ė = −C0r− 1 2 , (3.35) representing a balance between the quadratic stress and upper convective derivative terms. We may use this to give estimates of the slip boundary layer thickness. The polymer stress enters the slip layer with behavior (3.34), which on using with (3.31) in (3.35) gives π − θ ∼ r− 3(1+n1) 2n1 =   r 3 20 , PTT, r 3 14 , Giesekus, (3.36) or equivalently in Cartesian coordinates π − θ ∼ y/(−x) and r ∼ (−x), y ∼   (−x) 23 20 , PTT, (−x) 17 14 , Giesekus. (3.37) It is clear that no such balance is possible for the Oldroyd-B model since the quadratic stress terms are not present and the relaxation and rate-of-strain terms in the first equation in (3.33) remain subdominant for small r. However, a slip boundary layer still occurs, but its scaling now comes from balancing Tp θθ and the rate-of-strain terms in the third equation in (3.33). Using (3.34), this gives π − θ ∼ r, (3.38) or equivalently y ∼ (−x)2. In this case, for small r, there is thus no mechanism in the first equation in (3.33) to arrest the growth of the normal radial stress component. As a consequence, we will show in Sec. IV that the stress is actually infinite on the slip surface. B. Numerical results We now present numerical results for the system of Eq. (3.3). For a numerical solution, we fix the streamline, i.e., the value of ψ in (3.1), and solve the resulting sys- tem of ODEs. We start sufficiently far upstream, impos- ing viscometric stresses obtained from solving (3.12). The interval of integration for θ is taken as [10�6, π � 10�10] and use Matlab’s ode15s solver with AbsTol = RelTol = 10�6. Parameter values in all simulations were κ = 0.1 for PTT and Giesekus. For convenience, we take C0 = �1, so that the stream function takes values relative to this normalisation. 121604-8 Evans, Palhares Junior, and Oishi Phys. Fluids 29, 121604 (2017) FIG. 2. A plot of selected streamlines using (3.1) near the singularity, along with the boundary layer curves for all three models at both the stick and slip surfaces. Also shown is the line θ ≈ 1.91 where the radial dis- tance is minimum and radial velocity vanishes on each geometrically similar streamline. We note for the stream function (3.1) that the minimum radial distance along a given streamline occurs when f ′(θ) = 0 giving θ = 2 tan−1 √ 2 ≈ 1.9107 or 109.5◦. This is also the location at which the radial velocity component vanishes. To orient ourselves, Fig. 2 shows selected streamlines that pass close to the singularity. Plotted are the stick and slip boundary layer curves from (3.21), (3.36), and (3.38). Evi- dent is their cusp like nature, with PTT wider than Giesekus FIG. 3. Stress and velocity components along the streamline ψ = �10�6. The top row [(a)–(c)] shows the polar stress components for all three models. The middle row [(d)–(f)] gives the components the dyadic tensor vvT arising in the stretching solution (3.22). The bottom row [(g)–(i)] gives estimates of λ formed from the ratio of the appropriate components of stress and the dyadic product of the velocities. 121604-9 Evans, Palhares Junior, and Oishi Phys. Fluids 29, 121604 (2017) at both surfaces and Oldroyd-B giving the thinnest. The slip boundary layers are marginally thicker than the correspond- ing stick layers for PTT and Giesekus. The reverse is true for Oldroyd-B, with the slip layer being substantially thinner than the stick layer. Figure 3 gives the polymer stress profiles for all three models along the streamline ψ = �10�6, together with the FIG. 4. Verification of the singularity behavior (3.29) for the Oldroyd-B, PTT, and Giesekus models along the line θ = π/2. components of the dyadic product of velocity vvT and numer- ical estimates for λ from the computed ratio of the stresses and products of the velocity components. Noticeable is the difference in behavior between the models of the normal radial stress component Tp rr as the slip surface is approached (θ → π): FIG. 5. Behavior of Tp rr along the streamline ψ = �10�6 with the radial dis- tance. Also plotted are the stick and slip layer behaviors, the core stretching solution, and estimates of the demarcation between these behaviors given by the vertical lines obtained from the boundary layer scalings. 121604-10 Evans, Palhares Junior, and Oishi Phys. Fluids 29, 121604 (2017) in the Oldroyd-B model, it continues to grow following v2 r from the stretching solution, whilst in the PTT and Giesekus models its growth is arrested and reduced. The other stress components are also shown, with the values of the Oldroyd-B model having to be reduced 40 times to fit on the same plots in (b) and (c). FIG. 6. Behavior of Tp rr now along the streamline ψ = �10�12. Also plot- ted are the stick high shear rate viscometric behaviors as well as the slip layer behaviors and the core stretching solution. The vertical lines give esti- mates of the boundary layer edges and transition to and from core stretching behavior. This illustrates the significantly higher stresses given by the Oldroyd-B model. To confirm the stretching solution (3.22), estimates of λ from all three stress components are also shown in (g)–(i) and give a consistent constant value in the main part of the domain away from the stick θ = 0 and slip θ = π surfaces: deviation occurs in the plots (g) and (h) at θ ≈ 109.5o where the radial velocity vanishes, the ratios involving 3r not being reliable near this point. Figure 4 confirms the singularity behaviors (3.29) for each of the three models along the line θ = π/2. Finally, Figs. 5 and 6 examine the stress behaviors near the stick and slip surfaces. In Fig. 5, the normal radial stress component Tp rr is plotted with the radial distance along the streamline ψ = �10�6. In all figures, the lower curve corresponds to stress values at locations before the minimum radial distance, i.e., θ ≤ 1.9107 ≈ 109.5◦ and thus emanating close to the stick surface, whilst the upper curve represents locations after the minimum radial distance with θ ≥ 1.9107 ≈ 109.5◦ and thus close to the slip surface. Plotted for each model is the stretching solution along with their viscometric stick layer behavior and their slip layer behavior. The verti- cal lines give estimates of the location of the boundary layers at the stick surface from (3.21) and the slip surface from (3.36) and (3.38). In all three models, they give good estimates of when the stretching solution begins its dominance when leaving the stick surface region, with the transition between viscometric behavior and stretching behavior being clearly evident. They also give good estimates when the stretching solution no longer applies near the slip surface region for PTT and Giesekus. Noticeable is that the behavior of the mod- els near the slip surface is markedly different for PTT and Giesekus compared to Oldroyd-B. We see continued growth for Oldroyd-B in (a) with slope one for the green line, con- sistent with the radial index in the stretching behavior (3.34). This growth is not arrested within the slip surface layer, the location of which is influenced by the angular normal stress Tp θθ . For PTT and Giesekus, the slip region shows a clear transition to the behavior (3.35), the slope of the green lines being �1/2. Figure 6 shows similar information but along the streamline ψ = �10�12. Now the large shear rate viscometric behaviors are excellent approximations within the stick region (being less so on the streamline ψ = �10�6 and not shown). The continued stress growth of Oldroyd-B and its arrest for PTT and Giesekus are again strikingly evident in the slip region. IV. ASYMPTOTICS OF THE OLDROYD-B MODEL IN FIXED KINEMATICS We now focus solely upon the Oldroyd-B model and describe here the matched asymptotic solution near the sin- gularity. The analysis is performed in a Newtonian velocity field, so that the kinematics as given by (2.14) are fixed. We know, in contrast to PTT and Giesekus, that such a velocity field is not correct for the Oldroyd-B model near the singu- larity. However, the analysis does shed light on the behavior of the model, critical features of which we would expect to transfer over to the case of the true viscoelastic velocity field. Moreover, fixed Newtonian kinematics is a common test 121604-11 Evans, Palhares Junior, and Oishi Phys. Fluids 29, 121604 (2017) situation to consider first, as it amounts to simply omitting the contribution of the polymer stress term in the momen- tum equation. This is often a simpler problem for numer- ical schemes, particularly for those that modify Newtonian based solvers/strategies. Thus, analysing this situation gives further understanding to the numerical results of Rosenberg and Keunings23 and Fortin et al.26 and highlights that the Oldroyd-B model is not well-behaved in this often taken test situation. For later use, we state the Oldroyd-B polymer stress con- stitutive equation in (2.4) in both Cartesian and natural stress forms. The natural stress allows a complete solution to be obtained since the stress variables can be fully matched at leading order between the regions. The Cartesian formulation does not allow this without proceeding to higher terms in its expansion in the core region, which is prohibitive. Neverthe- less, we simultaneously perform the Cartesian formulation at leading order in the regions as it serves as a consistency check on the natural stress variables. In the Cartesian stress form, (2.4) is Tp 11 + ( u ∂Tp 11 ∂x + v ∂Tp 11 ∂y − 2 ∂u ∂x Tp 11 − 2 ∂u ∂y Tp 12 ) = 2 ∂u ∂x , Tp 22 + ( u ∂Tp 22 ∂x + v ∂Tp 22 ∂y − 2 ∂v ∂x Tp 12 − 2 ∂v ∂y Tp 22 ) = 2 ∂v ∂y , Tp 12 + ( u ∂Tp 12 ∂x + v ∂Tp 12 ∂y − ∂u ∂y T22 − ∂v ∂x Tp 11 ) = ∂u ∂y + ∂v ∂x . (4.1) Aligning the polymer stress field along streamlines39 gives the natural stress form λ + (v · ∇)λ + 2µ∇ · w = 1 |v|2 , µ + (v · ∇)µ + ν∇ · w = 0, ν + (v · ∇)ν = |v|2, (4.2) where v = ( u v ) , w = 1 |v|2 ( −v u ) , Tp + I = λvvT + µ(vwT + wvT ) + νwwT . (4.3) The polymer stress expression in (4.3) links the Cartesian (Tp 11, Tp 12, Tp 22) and natural stress components (λ, µ, ν) together. For later use, these are Tp 11 = −1 + λu2 − µ 2uv (u2 + v2) + ν v2 (u2 + v2)2 , Tp 12 = λuv + µ (u2 − v2) (u2 + v2) − ν uv (u2 + v2)2 , Tp 22 = −1 + λv2 + µ 2uv (u2 + v2) + ν u2 (u2 + v2)2 . (4.4) As described in Sec. III, we expect a three region structure near the singularity composed of a core region with narrow boundary layers at the stick and slip surfaces. We now address the details of these regions in turn. A schematic summary is shown in Fig. 7, where we use the parameter ε to represent the size of the small length scale on which this asymptotic structure pertains. A practical estimate of it can be obtained from the numerical results in Fig. 4, suggesting it should be 10�2 or smaller. A. The outer core solution Following Sec. III A 2, we expect the stretching solution Tp = λ(ψ)vvT , λ(ψ) = C1 C2 0 ( ψ C0 )− 6 5 as r → 0, (4.5) where C1 is an arbitrary constant, to apply near the singularity but away from the stick and slip surfaces. This leading order outer solution gives the estimates v = O(r 1 2 ), ∇v = O(r− 1 2 ), Ts = O(r− 1 2 ), Tp = O(r− 4 5 ) as r → 0. The assumption of the dominance of the solvent stress over the polymer stress, required for the velocity field (2.14), does not hold. Consequently, the momentum equation will not be satisfied, and we proceed on the understanding that we study the Oldroyd-B equations for the simpler problem in which the kinematics are fixed. The behavior of the natural stress variables in this region may be determined from (4.2), which at leading order are (v · ∇)λ = 0, (v · ∇)µ = 0, (v · ∇)ν = 0. (4.6) All three variables are thus constant along streamlines, with µ = C2 ( ψ C0 )− 1 5 , ν = C2 0C3 ( ψ C0 ) 4 5 , along with that for λ in (4.5) being determined by matching to the stick surface boundary layer (described next in Sec. IV B). FIG. 7. Asymptotic structure local to the singularity for small length scales ε . 121604-12 Evans, Palhares Junior, and Oishi Phys. Fluids 29, 121604 (2017) The three free constants C1, C2, and C3 are associated with each natural stress variable and communicate the neces- sary polymer stress information between the stick and slip boundary layers. The estimates λ = O(r− 9 5 ), µ = O(r− 3 10 ), ν = O(r 6 5 ) (4.7) may be used to confirm the dominance of the terms (4.6) within the constitutive equation (4.2). The Cartesian stress components for (4.5) are Tp 11 = λu2, Tp 12 = λuv , Tp 22 = λv 2. (4.8) We may now determine the limiting behaviors as the stick and slip surfaces are approached. This is required for matching between the regions and can also be used to determine where a change in dominant balance occurs leading to the required boundary layers at the stick and slip surfaces. Approaching the stick surface, we have as y → 0+, x > 0 : ψ ∼ C0x− 1 2 y2, u ∼ 2C0x− 1 2 y, v ∼ 1 2 C0x− 3 2 y2, T s 11 = −T s 22 ∼ −2C0x− 3 2 y, T s 12 ∼ 2C0x− 1 2 , Tp 11 ∼ 4C1x− 2 5 y− 2 5 , Tp 12 ∼ C1x− 7 5 y 3 5 , T22 ∼ 1 4 C1x− 12 5 y 8 5 , λ ∼ C1 C0 2 x 3 5 y− 12 5 , µ ∼ C2x 1 10 y− 2 5 , ν ∼ C2 0C3x− 2 5 y 8 5 . (4.9) The polymer and solvent stress components become the same size when Tp 12 = O(T s 12), giving y = O(x 3 2 ) as the thickness of the stick boundary layer. This is consistent with the scaling obtained in Sec. III A 1. For the slip surface, we have as y → 0+, x < 0 : ψ ∼ 2C0(−x) 1 2 y, u ∼ 2C0(−x) 1 2 , v ∼ C0(−x)− 1 2 y, T s 11 = −T s 22 ∼ −2C0(−x)− 1 2 , T s 12 ∼ 2C0(−x)− 3 2 y, Tp 11 ∼ 4C1(−x) 2 5 (2y)− 6 5 , Tp 12 ∼ C1(−x)− 3 5 (2y)− 1 5 , Tp 22 ∼ 1 4 C1(−x)− 8 5 (2y) 4 5 , λ ∼ C1 C2 0 (−x)− 3 5 (2y)− 6 5 , µ ∼ C2(−x)− 1 10 (2y)− 1 5 , ν ∼ C2 0C3(−x) 2 5 (2y) 4 5 . (4.10) We see that the polymer Tp 11 and Tp 12 components dominate their corresponding solvent values. A change in dominance is only possible when Tp 22 = O(T s 22), giving y = O((�x)2) as the thickness of the slip boundary layer, again agreeing with the estimate obtained in Sec. III A 3. B. The stick surface boundary layer To systematically compare terms, we introduce a small artificial parameter ε (not to be confused with the PTT model parameter), through the scaling x = εX̄. (4.11) The scalings for the other variables are then y = ε 3 2 Ȳ , ψ = ε 5 2 Ψ̄, u = ε ū, v = ε 3 2 v̄ , Tp 11 = ε −1T̄p 11, Tp 12 = ε − 1 2 T̄p 12, Tp 22 = T̄p 22, λ = ε−3 λ̄, µ = ε− 1 2 µ̄, ν = ε2 ν̄, which follow from (4.9) and dominant balance in the con- stitutive equation (4.1) or (4.2). Thus, for the given stream function, Ψ̄ = C0X̄− 1 2 Ȳ2, (4.12) the leading order (in ε) polymer Cartesian stress boundary layer equations are T̄p 11 + * , ∂Ψ̄ ∂Ȳ ∂T̄p 11 ∂X̄ − ∂Ψ̄ ∂X̄ ∂T̄p 11 ∂Ȳ − 2 ∂2Ψ̄ ∂X̄∂Ȳ T̄p 11 − 2 ∂2Ψ̄ ∂Ȳ2 T̄p 12 + - = 0, T̄p 22 + * , ∂Ψ̄ ∂Ȳ ∂T̄p 22 ∂X̄ − ∂Ψ̄ ∂X̄ ∂T̄p 22 ∂Ȳ + 2 ∂2Ψ̄ ∂X̄∂Ȳ T̄p 22 + 2 ∂2Ψ̄ ∂X̄2 T̄p 12 + - = −2 ∂2Ψ̄ ∂X̄∂Ȳ , T̄p 12 + * , ∂Ψ̄ ∂Ȳ ∂T̄p 12 ∂X̄ − ∂Ψ̄ ∂X̄ ∂T̄p 12 ∂Ȳ + ∂2Ψ̄ ∂X̄2 T̄p 11 − ∂2Ψ̄ ∂Ȳ2 T̄p 22 + - = ∂2Ψ̄ ∂Ȳ2 , (4.13) subject to as Ȳ → 0+, T̄p 11 ∼ 8C2 0 X̄−1, T̄p 12 ∼ 2C0X̄− 1 2 , T̄p 22 ∼ 2C0X̄− 3 2 Ȳ , (4.14) as Ȳ → +∞, T̄p 11 ∼ 4C1X̄− 2 5 Ȳ− 2 5 , T̄p 12 ∼ C1X̄− 7 5 Ȳ 3 5 , T̄p 22 ∼ 1 4 C1X̄− 12 5 Ȳ 8 5 . (4.15) We recognise Eq. (4.13) as the high Weissenberg boundary layer equations of Renardy,40 which clearly manifest them- selves at solid surfaces near a singularity even for Wi = 1. The conditions (4.14) give viscometric behavior, and (4.15) are the matching conditions with the outer core solution (4.9). These equations possess the similarity solution ξ = |C0 | Ȳ X̄ 3 2 , T̄p 11 = C2 0 X̄−1tp 11(ξ), T̄p 12 = |C0 |X̄ − 1 2 tp 12(ξ), T̄p 22 = tp 22(ξ), which allows for changes in the sign of C0 and hence con- veniently covers the cases of both flow directions, stick-slip C0 < 0 and slip-stick C0 > 0. We thus obtain 5 2 ξ2tp 11 ′ + 4t12 − |C0 | C0 t11 = 0, 5 2 ξ2tp 22 ′ + 2ξt22 − 3 2 ξ2t12 − |C0 | C0 t22 = −2ξ, 5 2 ξ2tp 12 ′ + ξt12 + 2t22 − 3 4 ξ2t11 − |C0 | C0 t12 = −2, (4.16) as ξ → 0+, tp 11 ∼ 8, tp 12 ∼ 2 C0 |C0 | , tp 22 ∼ 2 C0 |C0 | ξ, (4.17) as ξ → +∞, tp 11 ∼ 4C∗1ξ − 2 5 , tp 12 ∼ C∗1ξ 3 5 , tp 22 ∼ 1 4 C∗1ξ 8 5, (4.18) 121604-13 Evans, Palhares Junior, and Oishi Phys. Fluids 29, 121604 (2017) where ′ denotes d/dξ and C∗1 = C1 |C0 | − 8 5 . (4.19) In natural stress variables, the analogous statement for the boundary layer equation (4.13) with imposed viscometric wall behavior and outer matching conditions is (v̄ · ∇̄)λ̄ + 2 µ̄ ∂ ∂Ȳ (1 ū ) + λ̄ = 0, (v̄ · ∇̄)µ̄ + ν̄ ∂ ∂Ȳ (1 ū ) + µ̄ = 0, (v̄ · ∇̄)ν̄ + ν̄ − ū2 = 0, (4.20) as Ȳ → 0+, λ̄ ∼ 2 Ȳ2 , µ̄ ∼ 2C0X̄− 1 2 , ν̄ ∼ 4C2 0 X̄−1Ȳ2, (4.21) as Ȳ → +∞, λ̄ ∼ C1 C2 0 X̄ 3 5 Ȳ− 12 5 , µ̄ ∼ C2X̄ 1 10 Ȳ− 2 5 , ν̄ ∼ C2 0C3X̄− 2 5 Ȳ 8 5 . (4.22) Equations (4.20) follow immediately from (4.2) or can be deduced from (4.13) using T̄p 11 = λ̄ū2, T̄p 12 = λ̄ūv̄ + µ̄, T̄p 22 = −1 + λ̄v̄2 + 2 µ̄ v̄ ū + ν̄ ū2 , (4.23) where ū = ∂Ψ̄ ∂Ȳ = 2C0X̄− 1 2 Ȳ , v̄ = − ∂Ψ̄ ∂X̄ = 1 2 C0X̄− 3 2 Ȳ2. Equations (4.23) themselves follow from (4.4). In self-similar form, we have λ̄ = C2 0 X̄−3 λ̃(ξ), µ̄ = |C0 |X̄ − 1 2 µ̃(ξ), ν̄ = X̄2 ν̃(ξ), and the system (4.20)–(4.22) now becoming 5 2 ξ2 λ̃ ′ + 6ξλ̃ + µ̃ ξ2 − λ̃ |C0 | C0 = 0, 5 2 ξ2 µ̃′ + 3µ̃ + ν̃ 2ξ2 − µ̃ |C0 | C0 = 0, 5 2 ξ2 ν̃′ − 4ξν̃ − |C0 | C0 (ν̃ − 4ξ2) = 0, (4.24) as ξ → 0+, λ̃ ∼ 2 ξ2 , µ̃ ∼ 2 C0 |C0 | , ν̃ ∼ 4ξ2, (4.25) as ξ → +∞, λ̃ ∼ C∗1ξ − 12 5 , µ̃ ∼ C∗2ξ − 2 5 , ν̃ ∼ C∗3ξ 8 5 , (4.26) with C∗1 as given in (4.19) and introducing the scaled far-field parameters, C∗2 = C2 |C0 | − 3 5 , C∗3 = C3 |C0 | 2 5 . (4.27) The relationships (4.23) in the similarity form are tp 11 = 4ξ2 λ̃, tp 12 = ξ 3 λ̃+ µ̃, tp 22 = −1+ 1 4 ξ4 λ̃+ 1 2 ξ µ̃+ ν̃ 4ξ2 , which link the Cartesian and natural stress formulations together. These suggest that the leading order far-field behav- ior (4.18) in the Cartesian statement can be replaced with the more accurate expressions tp 11 ∼ 4C∗1ξ − 2 5 , tp 12 ∼ C∗1ξ 3 5 + C∗2ξ − 2 5 , tp 22 ∼ 1 4 C∗1ξ 8 5 + 1 2 C∗2ξ 3 5 − 1 + 1 4 C∗3ξ − 2 5 , and highlight how further terms are needed in the Cartesian formulation to furnish the same information. For the natural stress formulation, it is convenient to use the scaled variables `(ξ) = ξ2 λ̃, m(ξ) = µ̃, n(ξ) = ν̃ ξ2 , for which (4.24)–(4.26) become 5 2 ξ2`′ + ξ` + m − ` |C0 | C0 = 0, 5 2 ξ2m′ + ξm + 1 2 n − m |C0 | C0 = 0, 5 2 ξ2n′ + ξn − |C0 | C0 (n − 4) = 0, (4.28) as ξ → 0+, ` = 2, m = 2 C0 |C0 | , n = 4, (4.29) as ξ → +∞, ` ∼ C∗1ξ − 2 5 , m ∼ C∗2ξ − 2 5 , n ∼ C∗3ξ − 2 5 . (4.30) Figure 8 shows the profiles and convergence to the far-field behaviors for the Cartesian system (4.16)–(4.18) and natural stress system (4.28)–(4.30) solved as initial-value-problems (IVPs) with initial data (4.17) and (4.29), respectively. Here we take C0 = �1, relevant to the stick-slip regime and again use Matlab solver ode15s with AbsTol = RelTol = 10�6. Both formulations give consistent estimates of C∗1 , the far-field estimates of the constants being C∗1 = 0.990 89, C∗2 = −1.23 861, C∗3 = 4.128 71. The core solution constants C1, C2, and C3 are thus fully deter- mined with (4.19) and (4.27) giving their scalings with the stream function constant C0. We remark that we could in prin- ciple consider the reverse flow situation of slip-stick for which we would take C0 = 1. Again (4.28)–(4.30) are solved as an IVP, but now initial data are posed in the far-field through (4.30) and the equations integrated into the wall to obtain the behavior (4.29). As we will see in Sec. IV C, the Oldroyd-B model is not well-posed at the slip surface for this problem and so this regime is of little interest. C. The slip surface boundary layer The boundary layer variables at the slip surface are given by the scalings x = εX̄ , y = ε2Ȳ , ψ = ε 5 2 Ψ̄, u = ε 1 2 ū, v = ε 3 2 v̄ , Tp 11 = ε −2T̄p 11, Tp 12 = ε −1T̄p 12, Tp 22 = T̄p 22, λ = ε−3 λ̄, µ = ε− 1 2 µ̄, ν = ε 3 2 ν̄, (4.31) 121604-14 Evans, Palhares Junior, and Oishi Phys. Fluids 29, 121604 (2017) FIG. 8. Stick-slip boundary layer similarity solution solved as initial-value-problems with viscometric wall behavior as initial data. Cartesian profiles are shown in (a) with (b) illustrating convergence to their far-field matching behavior. The natural stress variables are shown in (c) with their far-field convergence in (d). where now x < 0 and ε > 0 again an artificial small parameter. At leading order in ε , the polymer Cartesian stress boundary layer equations are * , ∂Ψ̄ ∂Ȳ ∂T̄p 11 ∂X̄ − ∂Ψ̄ ∂X̄ ∂T̄p 11 ∂Ȳ − 2 ∂2Ψ̄ ∂X̄∂Ȳ T̄p 11 − 2 ∂2Ψ̄ ∂Ȳ2 T̄p 12 + - = 0, * , ∂Ψ̄ ∂Ȳ ∂T̄p 22 ∂X̄ − ∂Ψ̄ ∂X̄ ∂T̄p 22 ∂Ȳ + 2 ∂2Ψ̄ ∂X̄∂Ȳ T̄p 22 + 2 ∂2Ψ̄ ∂X̄2 T̄p 12 + - = −2 ∂2Ψ̄ ∂X̄∂Ȳ , * , ∂Ψ̄ ∂Ȳ ∂T̄p 12 ∂X̄ − ∂Ψ̄ ∂X̄ ∂T̄p 12 ∂Ȳ + ∂2Ψ̄ ∂X̄2 T̄p 11 − ∂2Ψ̄ ∂Ȳ2 T̄p 22 + - = 0, (4.32) with outer core matching conditions as Ȳ → +∞, T̄p 11 ∼ 4C1(−X̄) 2 5 (2Ȳ )− 6 5 , T̄p 12 ∼ C1(−X̄)− 3 5 (2Ȳ )− 1 5 , T̄p 22 ∼ 1 4 C1(−X̄)− 8 5 (2Ȳ ) 4 5 . (4.33) The leading order stream function here is given by Ψ̄ = 2C0(−X̄) 1 2 Ȳ , (4.34) for which (4.32) possesses an exact solution. This is most conveniently expressed in the self-similar form T̄p 11 = (−X̄)2tp 11(ξ), T̄p 12 = (−X̄)−1tp 12(ξ), T̄p 22 = tp 22(ξ), (4.35) where ξ = 2Ȳ (−X̄)2 , tp 11 = K1ξ − 6 5 , tp 12 = K1 4 ξ− 1 5 + K2ξ − 2 5 , tp 22 = K1 16 ξ 4 5 + K2 2 ξ 3 5 + K3ξ 2 5 − 1, and K1, K2, and K3 are arbitrary constants. The matching condition (4.33) determines K1 = 4C1. (4.36) 121604-15 Evans, Palhares Junior, and Oishi Phys. Fluids 29, 121604 (2017) The other two constants are taken to be K2 = K3 = 0, as suggested below using natural stress variables. Thus, as we approach the slip surface, we have the behavior as Ȳ → 0+, T̄p 11 ∼ 4C1(−X̄) 2 5 (2Ȳ )− 6 5 , T̄p 12 = C1(−X̄)− 3 5 (2Ȳ )− 1 5 , T̄p 22 ∼ −1. (4.37) Importantly, on the slip surface Ȳ = 0, both T̄p 11 and T̄p 12 are infinite. This is in complete contrast to the solution for the PTT and Giesekus models in Refs. 30 and 31, where the polymer stresses are finite (with the shear stress zero). In the Oldroyd- B model, we may interpret this as the stress singularity at the join of the slip and stick surfaces being extended along the slip surface. As such the model gives physically unrealistic polymer stress values on the free surface and does not seem to give a sensible practical solution in this flow situation. To complete the analysis, we determine the natural stress behavior in this region. In boundary layer variables (4.31), the natural stress equations are ε 1 2 λ̄ + (v̄ · ∇̄)λ̄ + ε 1 2 µ̄(∇̄ · w̄) = ε 5 2 (ū2 + ε2 v̄2) , ε 1 2 µ̄ + (v̄ · ∇̄)µ̄ + ν̄(∇̄ · w̄) = 0, ε 1 2 ν̄ + (v̄ · ∇̄)ν̄ = (ū2 + ε2 v̄2), (4.38) where ∇̄.w̄ = ∂ ∂Ȳ ( ū ū2 + ε2 v̄2 ) + ε2 ∂ ∂X̄ ( − v̄ (ū2 + ε2 v̄2) ) and matching conditions Ȳ → +∞, λ̄ ∼ C1 C2 0 (−X̄)− 3 5 (−2Ȳ )− 6 5 , µ̄ ∼ C2(−X̄)− 1 10 (2Ȳ )− 1 5 , ν̄ ∼ C2 0C3(−X̄) 2 5 (2Ȳ ) 4 5 . (4.39) Posing the expansions λ̄ ∼ λ0, µ̄ ∼ µ0, ν̄ ∼ ν0 + ε 1 2 ν1, ū ∼ u0, v̄ ∼ v0, as ε → 0, (4.38) give at leading order (v0 · ∇̄)λ0 = 0, (v0 · ∇̄)µ0 = 0, (v0 · ∇̄)ν0 = u2 0, (4.40) and at the next order, (v0 · ∇)ν1 = −ν0, (4.41) where v0 = ( u0 v0 ) , u0 = 2C0(−X̄) 1 2 , v0 = C0(−X̄)− 1 2 Ȳ . We note that ∇̄ · w̄ = o(1) since u0 is independent of Ȳ and care has been taken in deriving (4.41) to ensure that the first order velocity terms do not genuinely enter the equation. The solution for λ0 and µ0 is thus λ0 = C1 C2 0 ( Ψ̄ C0 )− 6 5 , µ0 = C2 ( Ψ̄ C0 )− 1 5 , where Ψ̄ is the leading order stream function (4.34). These natural stress variables are thus unchanged in this layer from their limiting outer core behaviors. Writing ν0 = C0(−X̄) 3 2 ν̃0(ξ), the third equation in (4.40) for ν0 gives 5ξ d ν̃0 dξ − 3ν̃0 = 4, with solution ν̃0 = − 4 3 + A0ξ 3 5 , for arbitrary constant A0. The matching condition (4.39) implies A0 = 0. Similarly for ν1 in (4.41), we have ν1 = (−X̄)2 ν̃1(ξ), with 5ξ d ν̃1 dξ − 4ν̃1 = 4 3 and the exact solution ν̃1 = − 1 3 + A1ξ 4 5 . (4.42) The matching condition (4.39) finally fixes the arbitrary constant A1 to be A1 = C2 0C3. (4.43) In summary, we thus have the solution λ̄ ∼ C1 C2 0 ( Ψ̄ C0 )− 6 5 , µ̄ ∼ C2 ( Ψ̄ C0 )− 1 5 , ν̄ ∼ − 4 3 C0(−X̄) 3 2 + ε 1 2  − 1 3 (−X̄)2 + C2 0C3 ( Ψ̄ C0 ) 4 5  . (4.44) The leading order forced solution for ν̄ is noteworthy, being a consequence of the slip velocity, with the core stretch- ing behavior only appearing at the next order. Further, the expressions (4.4) give T̄p 11 ∼ λ̄ū2 − ε2, T̄p 12 ∼ λ̄ūv̄ + ε 1 2 µ̄, T̄p 22 ∼ −1 + λ̄v̄2 + ε 1 2 ( 2 µ̄ v̄ ū + ν̄ ū2 ) , (4.45) and thus at leading order we have T̄p 11 = λ0u0 2, T̄p 12 = λ0u0v0, T̄p 22 = −1 + λ0v0 2, which recovers (4.35) with K2 = K3 = 0. The expressions (4.45) show how further expansion terms are needed in the Cartesian formulation to capture the information carried by the µ and ν natural stress variables. V. NATURAL STRESS ALONG STREAMLINES To verify the asymptotics of Sec. IV, we may consider the behavior of the natural stress variables along streamlines. This we do in very much the same manner of Sec. III, where for comparison we also present results for the PTT and Giesekus models. Using the construction (4.3), the extension of (4.2) to the PTT and Giesekus cases takes the form λ + (v · ∇)λ + 2µ∇ · w + κgλ = 1 |v|2 , µ + (v · ∇)µ + ν∇ · w + κgµ = 0, ν + (v · ∇)ν + κgν = |v|2, (5.1) 121604-16 Evans, Palhares Junior, and Oishi Phys. Fluids 29, 121604 (2017) FIG. 9. Profiles of natural stress variables with θ along the streamline ψ = −10−6. where gλ =   ( λ|v|2 − 2 + ν |v|2 ) ( λ − 1 |v|2 ) , PTT,( λ − 1 |v |2 )2 |v|2 + µ2 |v |2 , Giesekus, gµ = ( λ|v|2 − 2 + ν |v|2 ) µ, gν =   ( λ|v|2 − 2 + ν |v|2 ) ( ν − |v|2 ) , PTT,( ν − |v|2 )2 1 |v|2 + µ2 |v|2, Giesekus. (5.2) FIG. 10. Profiles of natural stress variables with θ along the streamline ψ = −10−12. As in Sec. III, we consider these equations along streamlines given by (3.1). Using (3.5) and (3.6), we have vθ r dλ dθ + 2µ∇ · w + κgλ + λ = 1 |v|2 , vθ r dµ dθ + ν∇ · w + κgµ + µ = 0, vθ r dν dθ + κgν + ν = |v|2, (5.3) where ∇ · w = 1 |v|4 ( (v2 θ − v 2 r ) ( ∂vθ ∂r + 1 r ∂vr ∂θ − vθ r ) + 4vrvθ ∂vr ∂r ) 121604-17 Evans, Palhares Junior, and Oishi Phys. Fluids 29, 121604 (2017) and |v|2 = v2 r + v2 θ . Central to the asymptotic solution is that the natural stress variables are constant along streamlines in the core region. Here we verify the behaviors FIG. 11. Verification of the exponents (5.5) for the variation of the natural stress variables across streamlines (5.4). The natural stress variables are plotted along the line θ = π 2 . The radial dependencies of (5.4) using the stream function (3.1) are 3/2 times the values in (5.5), which give the slopes plotted in the figure. λ = C1 C2 0 ( ψ C0 )n1 , µ = C2 ( ψ C0 )n2 , ν = C2 0C3 ( ψ C0 )n3 , (5.4) with the exponents for the three models being Oldroyd-B: n1 = − 6 5 , n2 = − 1 5 , n3 = 4 5 , PTT: n1 = − 10 11 , n2 = − 1 11 , n3 = 8 11 , Giesekus: n1 = − 7 8 , n2 = 0, n3 = 7 8 . (5.5) Figures 9 and 10 illustrate the numerical solution for the three natural stress variables along the streamlines ψ = �10�6 and ψ = �10�12 with κ = 0.1 for PTT and Giesekus. The pro- files are relatively constant in the main part of the domain and clearly deviate as expected near the stick and slip surface regions. There is a large variation between the values that the three natural stress variables take in all three models, with λ very large, ν very small, and µ between these (being more order 1). Where necessary, the values of the Oldroyd-B model have been scaled with the values shown in the legends to fit on the same graph. The ν variation for the PTT model is less convincingly constant than the other two models but improves with smaller and more realistic κ values such as 0.01 (not shown). Figure 11 confirms the exponents (5.5) for (5.4) along the line θ = π/2, where the radial exponents in (5.4) will be the respective powers multiplied by 3/2, on using the stream func- tion (3.1). The excellent agreement between the asymptotics and numerical results is evident. VI. DISCUSSION The asymptotic results in Refs. 30 and 31 for the stick- slip stress singularity of the PTT and Giesekus models have been validated numerically in a given Newtonian velocity field, along with their behaviors at the stick and slip bound- aries. However, it still remains to confirm these results through full numerical simulation, particularly the assumed Newto- nian velocity behavior near the singularity. The results of Sec. III B indicate that the stress singularity is caught for all models on radial distances smaller than 10�2 when Wi = 1. The scaling (2.15) suggests that this becomes 10−2Wi2 for other Weissenberg values. Thus for small Weissenberg numbers, the length scale of the singularity decreases and its effect is confined to progressively smaller distances as the Weissenberg number decreases. For large Weissenberg numbers, the opposite occurs with its length scale increas- ing to O(1) when the Weissenberg number is only of order 10. Another feature apparent in the numerical results of Sec. III B is that the PTT and Giesekus polymer stresses are far smaller (by two or more orders of magnitude) than those of Oldroyd-B, even away from the slip surface (where Oldroyd-B grows unboundedly). This undoubtedly is an important aspect which has allowed successful simulation of PTT in Refs. 26 and 27. 121604-18 Evans, Palhares Junior, and Oishi Phys. Fluids 29, 121604 (2017) The corresponding asymptotic results for the Oldroyd-B model have been derived in the regime of fixed Newtonian kinematics. Unlike the PTT and Giesekus models, this regime is not the flow field that occurs for the Oldroyd-B model near the singularity. This follows since the polymer stress is larger than the solvent stress and consequently the momentum equa- tion does not reduce to Stokes flow. Nevertheless, examining the behavior of the polymer stress equations has shed signifi- cant light on how the Oldroyd-B model responds in this flow situation. The presence of unbounded stress on the slip sur- face indicates that the model is not well behaved for these types of flows and parallels the infinite stresses that can be obtained in elongational flows. The analysis and numerics of Secs. III and IV for the Oldroyd-B model show that it appears to have no mechanism within its equations to arrest the growth of stress at the slip surface. This would explain the severe stresses seen in the numerical schemes of Refs. 23 and 26. This is in complete contrast to the PTT and Giesekus models, where it is the quadratic stress terms which stop the continual stress growth associated with the convected and affine defor- mation terms of the upper convected derivative of stress. It is of course still unknown if the stress response of the Oldroyd-B model behaves any better in viscoelastic velocity kinematics. However, if the model has a singularity at the stick-slip join, then as in the Newtonian velocity field, it is difficult to see a mechanism in the Oldroyd-B equations that can arrest the growth of its polymer stress along the slip surface. The numer- ical results in the literature are mixed, with Ref. 27 suggesting convergence can be obtained for reasonably large Weissenberg numbers up to 25, whereas the careful numerics in Ref. 24 suggest that this is not the case for the Oldroyd-B model. ACKNOWLEDGMENTS The authors would like to thank the financial support given by SPRINT/FAPESP Grant No. 2015/50094-7 and The Royal Society Newton International Exchanges Grant No. 2015/NI150225. C. M. Oishi acknowledges CNPq (Con- selho Nacional de Desenvolvimento Cientı́fico e Tecnologico), Grant No. 307459/2016-0. I. L. Palhares Junior acknowledges the financial support of FAPESP (Fundação de Amparo à Pesquisa do Estado de São Paulo) Grant No. 2014/17348-2 and BEPE No. 2016/20389-8. 1S. Richardson, “A ‘stick-slip’ problem related to the motion of a free jet at low Reynolds numbers,” Math. Proc. Cambridge Philos. Soc. 67, 477–489 (1970). 2R. I. Tanner, Engineering Rheology, 2nd ed. (Oxford University Press, 2000). 3O. Hassager, “Working group on numerical techniques in Proceedings of the Vth Workshop on Numerical Methods in Non-Newtonian Flow, Lake Arrowhead, USA,” J. Non-Newtonian Fluid Mech. 29, 2–5 (1988). 4N. P. Thien and R. I. Tanner, “A new constitutive equation derived from network theory,” J. Non-Newtonian Fluid Mech. 2(4), 353–365 (1977). 5N. Phan-Thien, “A nonlinear network viscoelastic model,” J. Rheol. 22(3), 259–283 (1978). 6H. Giesekus, “A simple constitutive equation for polymer fluids based on the concept of deformation-dependent tensorial mobility,” J. Non-Newtonian Fluid Mech. 11(1-2), 69–109 (1982). 7H. Giesekus, “A unified approach to a variety of constitutive models for polymer fluids based on the concept of configuration-dependent molecular mobility,” Rheol. Acta 21, 366–375 (1982). 8J. G. Oldroyd, “On the formulation of rheological equations of state,” Proc. R. Soc. A 200, 523–541 (1950). 9R. B. Bird, R. C. Armstrong, and O. Hassager, Dynamics of Polymeric Liquids, Volume 1: Fluid Mechanics (John Wiley & Sons, New York, 1977). 10D. V. Boger, “A highly elastic constant-viscosity fluid,” J. Non-Newtonian Fluid Mech. 3(1), 87–91 (1977). 11M. E. Mackay and D. V. Boger, “An explanation of the rheological prop- erties of Boger fluids,” J. Non-Newtonian Fluid Mech. 22(2), 235–243 (1987). 12M. Renardy, Mathematical Analysis of Viscoelastic Flows (SIAM, 2000). 13S. A. Trogdon and D. D. Joseph, “The stick-slip problem for a round jet I. Large surface tension,” Rheol. Acta 19(4), 404–420 (1980). 14C. J. Coleman, “A contour integral formulation of plane creeping Newtonian flow,” Q. J. Mech. Appl. Math. 34(4), 453–464 (1981). 15G. C. Georgiou, L. G. Olson, W. W. Schultz, and S. Sagan, “A singular finite element for Stokes flow: The stick–slip problem,” Int. J. Numer. Methods Fluids 9(11), 1353–1367 (1989). 16T. R. Salamon, D. E. Bornside, R. C. Armstrong, and R. A. Brown, “The role of surface tension in the dominant balance in the die swell singularity,” Phys. Fluids 7(10), 2328–2344 (1995). 17D. H. Michael, “The separation of a viscous liquid at a straight edge,” Mathematika 5(1), 82–84 (1958). 18H. K. Moffatt, “Viscous and resistive eddies near a sharp corner,” J. Fluid Mech. 18(1), 1–18 (1964). 19M. Elliotis, G. Georgiou, and C. Xenophontos, “Solution of the pla- nar Newtonian stick–slip problem with the singular function bound- ary integral method,” Int. J. Numer. Methods Fluids 48(9), 1001–1021 (2005). 20C. J. Coleman, “On the use of boundary integral methods in the analysis of non-Newtonian fluid flow,” J. Non-Newtonian Fluid Mech. 16(3), 347–355 (1984). 21J. M. Marchal and M. J. Crochet, “A new mixed finite element for cal- culating viscoelastic flow,” J. Non-Newtonian Fluid Mech. 26(1), 77–114 (1987). 22M. R. Apelian, R. C. Armstrong, and R. A. Brown, “Impact of the con- stitutive equation and singularity on the calculation of stick-slip flow: The modified upper-convected Maxwell model (MUCM),” J. Non-Newtonian Fluid Mech. 27(3), 299–321 (1988). 23J. Rosenberg and R. Keunings, “Numerical integration of differential viscoelastic models,” J. Non-Newtonian Fluid Mech. 39(3), 269–290 (1991). 24R. G. Owens and T. N. Phillips, “A spectral domain decomposition method for the planar non-Newtonian stick-slip problem,” J. Non-Newtonian Fluid Mech. 41(1-2), 43–79 (1991). 25T. R. Salamon, D. E. Bornside, R. C. Armstrong, and R. A. Brown, “Local similarity solutions for the stress field of an Oldroyd-B fluid in the partial- slip/slip flow,” Phys. Fluids 9(8), 2191–2209 (1997). 26A. Fortin, A. Zine, and J.-F. Agassant, “Computing viscoelastic fluid flow problems at low cost,” J. Non-Newtonian Fluid Mech. 45(2), 209–229 (1992). 27F. P. T. Baaijens, “Application of low-order discontinuous Galerkin methods to the analysis of viscoelastic flows,” J. Non-Newtonian Fluid Mech. 52(1), 37–57 (1994). 28V. Ngamaramvaranggul and M. F. Webster, “Viscoelastic simulations of stick-slip and die-swell flows,” Int. J. Numer. Methods Fluids 36(5), 539–595 (2001). 29G. Karapetsas and J. Tsamopoulos, “On the stick-slip flow from slit and cylindrical dies of a Phan-Thien and Tanner fluid model. I. Steady state,” Phys. Fluids 21(12), 123101 (2009). 30J. D. Evans, “Stick-slip and slip-stick singularities of the Phan-Thien– Tanner fluid,” J. Non-Newtonian Fluid Mech. 199, 12–19 (2013). 31J. D. Evans, “Stick-slip singularity of the Giesekus fluid,” J. Non-Newtonian Fluid Mech. 222, 24–33 (2015). 32R. I. Tanner and X. Huang, “Stress singularities in non-Newtonian stick- slip and edge flows,” J. Non-Newtonian Fluid Mech. 50(2-3), 135–160 (1993). 33M. Renardy, “The stresses of an upper convected Maxwell fluid in a Newto- nian velocity field near a re-entrant corner,” J. Non-Newtonian Fluid Mech. 50(2-3), 127–134 (1993). 34E. J. Hinch, “The flow of an Oldroyd fluid around a sharp corner,” J. Non- Newtonian Fluid Mech. 50(2-3), 161–171 (1993). 35J. M. Rallison and E. J. Hinch, “The flow of an Oldroyd fluid past a reentrant corner: The downstream boundary layer,” J. Non-Newtonian Fluid Mech. 116(2), 141–162 (2004). https://doi.org/10.1017/s0305004100045758 https://doi.org/10.1016/0377-0257(77)80021-9 https://doi.org/10.1122/1.549481 https://doi.org/10.1016/0377-0257(82)85016-7 https://doi.org/10.1016/0377-0257(82)85016-7 https://doi.org/10.1007/bf01534296 https://doi.org/10.1098/rspa.1950.0035 https://doi.org/10.1098/rspa.1950.0035 https://doi.org/10.1016/0377-0257(77)80014-1 https://doi.org/10.1016/0377-0257(77)80014-1 https://doi.org/10.1016/0377-0257(87)80038-1 https://doi.org/10.1007/bf01524013 https://doi.org/10.1093/qjmam/34.4.453 https://doi.org/10.1002/fld.1650091105 https://doi.org/10.1002/fld.1650091105 https://doi.org/10.1063/1.868746 https://doi.org/10.1112/s0025579300001388 https://doi.org/10.1017/s0022112064000015 https://doi.org/10.1017/s0022112064000015 https://doi.org/10.1002/fld.973 https://doi.org/10.1016/0377-0257(84)85019-3 https://doi.org/10.1016/0377-0257(87)85048-6 https://doi.org/10.1016/0377-0257(88)85002-x https://doi.org/10.1016/0377-0257(88)85002-x https://doi.org/10.1016/0377-0257(91)80018-f https://doi.org/10.1016/0377-0257(91)87035-v https://doi.org/10.1016/0377-0257(91)87035-v https://doi.org/10.1063/1.869342 https://doi.org/10.1016/0377-0257(92)85004-g https://doi.org/10.1016/0377-0257(94)85057-7 https://doi.org/10.1002/fld.145.abs https://doi.org/10.1063/1.3271495 https://doi.org/10.1016/j.jnnfm.2013.06.001 https://doi.org/10.1016/j.jnnfm.2014.08.012 https://doi.org/10.1016/j.jnnfm.2014.08.012 https://doi.org/10.1016/0377-0257(93)80028-a https://doi.org/10.1016/0377-0257(93)80027-9 https://doi.org/10.1016/0377-0257(93)80029-b https://doi.org/10.1016/0377-0257(93)80029-b https://doi.org/10.1016/j.jnnfm.2003.10.001 121604-19 Evans, Palhares Junior, and Oishi Phys. Fluids 29, 121604 (2017) 36J. D. Evans, “Re-entrant corner flows of Oldroyd-B fluids,” Proc. R. Soc. A 461, 2573–2603 (2005). 37G. Schleiniger and R. J. Weinacht, “A remark on the Giesekus viscoelastic fluid,” J. Rheol. 35(6), 1157–1170 (1991). 38G. Schleiniger and R. J. Weinacht, “Steady Poiseuille flows for a Giesekus fluid,” J. Non-Newtonian Fluid Mech. 40(1), 79–102 (1991). 39M. Renardy, “How to integrate the upper convected Maxwell (UCM) stresses near a singularity (and maybe elsewhere, too),” J. Non-Newtonian Fluid Mech. 52(1), 91–95 (1994). 40M. Renardy, “High Weissenberg number boundary layers for the upper convected Maxwell fluid,” J. Non-Newtonian Fluid Mech. 68(1), 125–132 (1997). https://doi.org/10.1098/rspa.2004.1410 https://doi.org/10.1122/1.550169 https://doi.org/10.1016/0377-0257(91)87027-u https://doi.org/10.1016/0377-0257(94)85060-7 https://doi.org/10.1016/0377-0257(94)85060-7 https://doi.org/10.1016/s0377-0257(96)01491-7