We verify numerically the theoretical stress singularities for two viscoelastic models that occur at sharp corners. The models considered are the Giesekus and Phan-Thien–Tanner (PTT), both of which are shear thinning and are able to capture realistic polymer behaviors. The theoretical asymptotic behavior of these two models at sharp corners has previously been found to involve an integrable solvent and polymer elastic stress singularity, along with narrow elastic stress boundary layers at the walls of the corner. We demonstrate here the validity of these theoretical results through numerical simulation of the classical contraction flow and analyzing the $270\xb0$ corner. Numerical results are presented, verifying both the solvent and polymer stress singularities, as well as the dominant terms in the constitutive equations supporting the elastic boundary layer structures. For comparison at Weissenberg order one, we consider both the Cartesian stress formulation and the alternative natural stress formulation of the viscoelastic constitutive equations. Numerically, it is shown that the natural stress formulation gives increased accuracy and convergence behavior at the stress singularity and, moreover, encounters no upper Weissenberg number limitation in the global flow simulation for sufficiently large solvent viscosity fraction. The numerical simulations with the Cartesian stress formulation cannot reach such high Weissenberg numbers and run into convergence failure associated with the so-called high Weissenberg number problem.

## I. INTRODUCTION

The construction of approximate solutions through asymptotic techniques is widely adopted in fluid mechanics, for example, van Dyke (1964), Lagerstrom (1988), Ockendon and Ockendon (1995), and Schlichting and Gersten (2000). It is particularly useful in identifying and quantifying singular behavior, such as boundary layers and stress singularities. In a mathematical context, asymptotic analysis is a powerful technique, often providing accurate approximate solutions where numerical schemes can encounter difficulty or even fail. Here, we consider such a situation that arises for viscoelastic fluids in a contraction domain that possesses sharp corners, where both stress singularities and boundary layers are encountered.

A variety of fluid flow problems include sharp corners at the boundaries of the flow domain. Examples of such cases include flows around forward- or backward-facing steps (Demuren and Wilson, 1994; Wang, 1994), cavity flows (both over and inside) (Gatski and Grosch, 1985), and flows in expanding or contracting channels (Fearn *et al.*, 1990; Shapira *et al.*, 1990; Battaglia *et al.*, 1997; Alleborn *et al.*, 1997; and Drikakis, 1997). These situations have received much attention for Newtonian fluids, where the sharp change in the slope of the boundary at the corner causes singularities in both the vorticity and stress. The identification of the singularity goes back to the work of Dean and Montagnon (1949) and Moffatt (1964). For Newtonian fluids, these corner singularities not only dominate the flow behavior near the corner, but also can extend to the entire domain, where reversed flow regions and flow separation can occur. Such sharp corner singularities can also create serious numerical difficulties in accurately resolving the flow field around the corner. Convergence of numerical schemes is often suppressed, unless specialized techniques, such as mesh refinement or singular elements, are used. Moreover, these difficulties may affect the accuracy of the overall numerical simulation in the whole computational domain. The approach of incorporating the stress singularity within an analytical expansion (Phillips, 1989) or numerical scheme (Holstein and Paddon, 1981; Floryan and Czechowski, 1995; Hawa *et al.*, 2002; Egger *et al.*, 2014; and Deliceoğlu *et al.*, 2019) for Newtonian fluids has been particularly successful.

The flow of viscoelastic fluids in such flow geometries though has posed far more problems. So much so that the contraction flow problem became a benchmark problem for the testing of numerical schemes (Hassager, 1988) back in 1987. A comprehensive discussion and reference for early work can be found in Owens and Phillips (2002), while the more recent literature can be found in Niethammer *et al.* (2018), Alves *et al.* (2021), Afonso *et al.* (2011), Fernandes *et al.* (2017), and Davoodi *et al.* (2022). The source of the difficulties can be apportioned to (1) a lack of knowledge of the stress singularity at sharp corners, and (2) the more challenging mathematical nature of the viscoelastic constitutive equations.

First, in regard to the stress singularity, it has only been relatively recently, since the work of Hinch (1993) and subsequent authors (Renardy, 1995, 1997; Rallison and Hinch, 2004; and Evans, 2005a; 2005b; 2008a; 2008b; 2010a; 2010b), that the asymptotic behavior has been determined for the UCM model, Oldroyd-B, PTT, and Giesekus. These models all share the common feature of narrow stress boundary layers at the walls, upstream and downstream of the singularity, which is a feature that Newtonian fluids do not possess. The Newtonian asymptotic stress behavior is uniformly valid around a sharp corner, a far simpler structure than the three-region asymptotic structure of the viscoelastic models. Second, the viscoelastic constitutive equations when coupled with the momentum equation give a mixed elliptic–hyperbolic system of equations (Joseph, 1990; Owens and Phillips, 2002; and Gerritsma and Phillips, 2001; 2008). This hyperbolic aspect of the viscoelastic equations ensures that the effect of any degradation in the numerical solution does not remain local, unlike Newtonian fluids (Renardy, 1989). This accounts for the immense difficulty that early approaches encountered. The problems of the stress singularity and narrow stress boundary layers are exacerbated as the relaxation time (elastic memory) in the viscoelastic models increases. This is notoriously referred to as the high Weissenberg number problem (HWNP) (Joseph *et al.*, 1985; Keunings, 1986; and Renardy, 2000a) and has been a major stumbling block in computational rheology for several decades (see, e.g., the reviews Keunings, 2000; Keunings and Euler, 2001; Walters and Webster, 2003; and Owens and Phillips, 2002). It refers to the inability of numerical schemes to resolve the stress singularities and boundary layers, causing computations to break down at relatively low Weissenberg numbers. The observed critical Weissenberg number at which convergence fails varies with both the numerical method, and geometry and rheological model. However, a significant step forward was made through the stabilization approach of Fattal and Kupferman (2004), who suggested that the HWNP was primarily a numerical issue of stiffness in capturing steep spatial gradients and their subsequent convection. The numerical consequence of inaccurately approximating transported stress profiles is the local multiplicative growth of the errors, resulting in a numerical blowup (Fattal and Kupferman, 2005). The approach to alleviate this numerical stiffness was to introduce the matrix logarithm of the stress (Fattal and Kupferman, 2004). Further discussion on this approach and subsequently other stabilization approaches can be found in Afonso *et al.* (2011; 2012) and Niethammer *et al.* (2018). However, this approach does not totally resolve the HWNP issue, and in order to make further progress, the stress singularity and boundary layers need to be adequately addressed by any numerical scheme.

In this current work, we consider two realistic viscoelastic models of Giesekus (1982) and Phan-Thien–Tanner (PTT) (1977), which capture the shear thinning behavior and non-zero normal stress differences of many common polymeric fluids (Renardy, 2000b). The stress singularity and boundary layer structures have been theoretically determined for these models in Renardy (1997), Hagen and Renardy (1997), and Evans (2010a; 2010b). We consider these two models together, since their constitutive equations differ by the presence of a term quadratic in the elastic stress [when the PTT model is considered in its “simplified” form of upper convected derivative and linear stress function (Poole *et al.*, 2019)]. However, the difference in this term leads to both different stress singularities and boundary layer structures for the two models. We also remark that these viscoelastic models are currently used in complex computational modeling (e.g., Minaeian *et al.*, 2022; Bayat *et al.*, 2022; Jeyasountharan *et al.*, 2022; and Maqbool *et al.*, 2022) evidencing the need of further study in the context of comparisons between theoretical and numerical results.

The goal of the work is to verify whether the theoretical stress singularity and boundary layer structures of the Giesekus and PTT models can be captured numerically in full numerical simulations. This would then confirm that the theoretical results are borne out in practical simulations and are the relevant structures to address in future work. We contrast two formulations of the constitutive equations. One in which the stress tensor is expressed in a fixed basis, which we term the Cartesian stress formulation (CSF), and another in a basis aligned with the velocity flow field, termed the natural stress formulation (NSF). We demonstrate that the NSF more accurately captures both the stress singularity and boundary layer structures, and moreover is not subject to a limiting Weissenberg number (for sufficiently large values of the solvent viscosity fraction). In addition, we have included a study concerning the dominance analysis of the constitutive equations around the reentrant corner. This investigation leverages further interpretations and constructions of data-driven methods (Callaham *et al.*, 2021) for viscoelastic models in complex flows.

We remark that because the stresses and deformation rates at a reentrant corner are infinite, the local Weissenberg number at the corner should also be viewed as infinite (Renardy, 2000a), regardless of the value of the global Weissenberg number. This consequently justifies performing simulations at order one Weissenberg numbers, where the main features of interest will still be manifest. The layout of the paper is as follows: in Secs. II and III, we summarize the models and asymptotic results, and then, in Sec. IV, we present the numerical simulations confirming the theoretical asymptotic behaviors.

## II. GOVERNING EQUATIONS

The dimensional governing equations in the present work are the mass and momentum equations given, respectively, by

where **v** is the fluid velocity, *p* is the pressure, *ρ* is the density, and $\tau $ is the extra-stress tensor. The extra-stress is written as the sum of Newtonian solvent and elastic (polymer) contributions in the form

where *η _{s}* is the solvent viscosity, $D=12(\u2207v+(\u2207v)T)$ is the rate-of-strain tensor, and $\tau p$ the elastic extra-stress. As constitutive equations for the elastic extra-stress, we adopt the PTT (the affine and linear stress function versions) and Giesekus models in the form

where *λ _{p}* is the relaxation time, $\tau p\u25bd$ is the upper-convected derivative of the elastic extra-stress,

*η*is the polymeric viscosity, and

_{p}*κ*is a model dependent parameter [often termed

*ϵ*for PTT and influences the elongational behavior (Phan-Thien and Tanner, 1977), while it is termed the mobility factor

*α*for Giesekus] controlling the influence of the quadratic stress terms

We non-dimensionalize using a characteristic length *L* and flow speed *U* as follows:

with $\eta 0=\eta s+\eta p$ being the total viscosity, which give (upon dropping ^{*}s) the dimensionless governing equations

with dimensionless parameters of the Reynolds number $Re=\rho UL/\eta 0$, the Weissenberg number $Wi=\lambda pU/L$ (the dimensionless relaxation time), and retardation parameter $\beta =\eta s/\eta 0\u2208[0,1]$ (a dimensionless retardation time or solvent viscosity fraction, so that $\eta p/\eta 0=1\u2212\beta $). Explicitly, the upper-convected stress derivative is defined as

The system of Eqs. (2.6)–(2.8) is considered here for planar flows, for which there are two notable representations of the constitutive equations (2.8), when expressed in component form (see, e.g., Evans and Oishi, 2017; Evans *et al.* 2019b; 2019a; 2020). The first and most common gives the Cartesian stress formulation (CSF), where the elastic extra-stress tensor is expressed with respect to the usual fixed Cartesian orthonormal basis, so that the Cartesian stress components for (2.8) satisfy

where

An alternative stress formulation introduces the conformation tensor **A** through the transformation

so that (2.8) becomes

The conformation tensor may then be decomposed using the velocity field and an orthogonal vector to express the elastic part of the extra-stress as

with

This is referred to as the natural stress formulation (NSF) and results in the components of (2.17) in (2.16) satisfying

where

and

It is notable that the coupling of the terms in the upper convective derivative part (square brackets) of the CSF equations (2.11)–(2.13) takes place through the spatial velocity gradients. In the equivalent terms in the NSF equations (2.18)–(2.20), the coupling takes place through the temporal velocity gradients. This emphasizes that the two formulations are significantly different and that the upper convected derivative of the natural stress variables achieves significant decoupling for steady flows.

## III. REVISITING THE ASYMPTOTIC RESULTS

In this section, we present a summary of the asymptotic results for the reentrant corner viscoelastic flows of the PTT and Giesekus models (Renardy, 1997; Evans, 2010a; 2010b). They are presented for the specific case of the $270\xb0$ corner angle, which arises in the contraction geometry shown in Fig. 1(a). The numerical simulations presented later in Sec. IV are obtained for the specific benchmark contraction case of 4:1.

The key feature of the asymptotic stress behaviors at a reentrant corner for both PTT and Giesekus models is that the Newtonian solvent stress dominates the elastic polymer stress. The velocity field is thus Newtonian in character and allows the classical works of Dean and Montagnon (1949) and Moffatt (1964) to be used to construct its behavior local to the reentrant corner. Consequently, away from the walls, but close to the corner, in a region we shall call the core region (Hinch, 1993), we can write the stream function and pressure, respectively, as

where *n* = 0.5445 and $(r,\theta )$ are the polar coordinates centered at the reentrant corner and

with

and *C*_{0} is a constant, which is determined by the properties of the flow away from the reentrant corner. In these expressions, we have used the numerical value of the first eigenvalue that controls the dominant behavior in the separable solution (3.1). This flow field gives the order of magnitude estimates determined in Evans (2010a; 2010b) and summarized in Table I for the main flow variables of velocity, velocity gradient (and pressure), elastic (polymer) stress, and natural stress variables.

Model . | v
. | $\u2207v$ (and p)
. | $Tp$ . | B.L. thickness . | λ
. | μ
. | ν
. |
---|---|---|---|---|---|---|---|

PTT | $r0.5445$ | $r\u22120.4555$ | $r\u22120.3286$ | $r1.1518$ | $r\u22121.4176$ | $r\u22120.1268$ | $r1.1638$ |

Giesekus | $r0.5445$ | $r\u22120.4555$ | $r\u22120.2796$ | $r1.2278$ | $r\u22121.3686$ | $r0.0$ | $r1.3686$ |

Model . | v
. | $\u2207v$ (and p)
. | $Tp$ . | B.L. thickness . | λ
. | μ
. | ν
. |
---|---|---|---|---|---|---|---|

PTT | $r0.5445$ | $r\u22120.4555$ | $r\u22120.3286$ | $r1.1518$ | $r\u22121.4176$ | $r\u22120.1268$ | $r1.1638$ |

Giesekus | $r0.5445$ | $r\u22120.4555$ | $r\u22120.2796$ | $r1.2278$ | $r\u22121.3686$ | $r0.0$ | $r1.3686$ |

The interesting features of the viscoelastic models are that the results in Table I for the elastic stresses do not hold uniformly around the corner. Rather, close to the walls viscometric elastic stresses dominate that are appropriate to shearing behavior. Consequently, stress boundary layers arise that reconcile the two behaviors. These stress boundary layers are the same that occur in the high Weissenberg limit discussed in Hagen and Renardy (1997), which appear to manifest themselves at stress singularities even in Weissenberg order one flows. The asymptotic estimates of the boundary layer (B.L.) widths are given in Table I. Figure 2 illustrates these boundary layer widths for the PTT and Giesekus models, along with three selected streamlines that pass close to the singularity. The PTT model is noted to have a slightly wider boundary layer than the Giesekus model, since its radial exponent for the B.L. thickness is smaller than that for Giesekus. A similar behavior in boundary layer thickness is seen in the stick-slip problem (Evans *et al.*, 2017).

The actual boundary layer equations for both models and formulations are now recorded below. These sets of equations assume a set of local Cartesian axes, orientated with *x* along the solid boundary and *y* normal into the fluid.

Boundary layer equations for PTT-CSF are as follows:

Boundary layer equations for PTT-NSF are as follows:

Boundary layer equations for Giesekus–CSF are as follows:

We remark that Eqs. (3.11) and (3.12) are more commonly written using $T\xaf22p=T22p+(1\u2212\beta )Wi$ and then take the more compact form

Boundary layer equations for Giesekus–NSF are as follows:

In both models, the elastic stresses of the CSF and NSF are linked in the boundary layers through the equations

A convenient summary of the terms that dominate in both the core and boundary layer regions, for both models and both formulations, is given in Table II. These are the leading-order terms that hold in the respective regions and it should be noted that not all terms in the quadratic stress expressions and the rate-of-strain term are present, as recorded in the explicit boundary layer Eqs. (3.4)–(3.17). The subscript on the upper convected derivative refers to the appropriate component and explicitly given in Eqs. (3.4)–(3.6) or (3.10)–(3.12). Moreover, it is noted for the Giesekus model that the constant terms in (3.12) arise from the leading order behavior $T22p\u223c\u2212(1\u2212\beta )/Wi$ in the linear and quadratic stress terms. Purely for convenience, we absorb the additional constant term from the linear stress component into the *g*_{22} quadratic stress term for the purposes of Table II. Such behavior of the *T*_{22} component in the boundary layer is peculiar to the Giesekus model and does not occur for PTT. These results hold for small or order one Reynolds number, order one (or large) Weissenberg numbers and $0<\beta <1$. We would expect these results to breakdown as the solvent viscosity vanishes. The case *β* = 0 thus needs separate consideration and for which currently no complete asymptotic results are known (Evans and Sibley, 2008; 2009).

Formulation . | Downstream and upstream BL . | Core . |
---|---|---|

CSF | $Wi[Tp11\u25bd+\kappa 1\u2212\beta g11]=0$ | |

$Wi[Tp12\u25bd+\kappa 1\u2212\beta g12]=2(1\u2212\beta )D12$ | $Tp\u25bd=0$ | |

$Wi[Tp22\u25bd+\kappa 1\u2212\beta g22]=2(1\u2212\beta )D22$ | ||

$g11\u223c(T22p)2,g12\u223cT11pT12p,g22\u223c{T11pT22p,PTT(T12p)2+(1\u2212\beta )2Wi2(\kappa \u22121)\kappa Giesekus$ | ||

NSF | $Wi[(v\xb7\u2207)\lambda +2\mu \u2207\xb7w]+\kappa g\lambda =0$ | $(v\xb7\u2207)\lambda =0$ |

$Wi[(v\xb7\u2207)\mu +\nu \u2207\xb7w]+\kappa g\mu =0$ | $(v\xb7\u2207)\mu =0$ | |

$Wi[(v\xb7\u2207)\nu ]+\kappa g\nu ={0,PTT|v|2,Giesekus$ | $(v\xb7\u2207)\nu =0$ | |

$g\lambda \u223c\lambda 2|v|2,g\mu \u223c\lambda \mu |v|2,g\nu \u223c{\lambda |v|2(\nu \u2212|v|2),PTT(\mu 2+\kappa \u22121\kappa )|v|2,Giesekus$ |

Formulation . | Downstream and upstream BL . | Core . |
---|---|---|

CSF | $Wi[Tp11\u25bd+\kappa 1\u2212\beta g11]=0$ | |

$Wi[Tp12\u25bd+\kappa 1\u2212\beta g12]=2(1\u2212\beta )D12$ | $Tp\u25bd=0$ | |

$Wi[Tp22\u25bd+\kappa 1\u2212\beta g22]=2(1\u2212\beta )D22$ | ||

$g11\u223c(T22p)2,g12\u223cT11pT12p,g22\u223c{T11pT22p,PTT(T12p)2+(1\u2212\beta )2Wi2(\kappa \u22121)\kappa Giesekus$ | ||

NSF | $Wi[(v\xb7\u2207)\lambda +2\mu \u2207\xb7w]+\kappa g\lambda =0$ | $(v\xb7\u2207)\lambda =0$ |

$Wi[(v\xb7\u2207)\mu +\nu \u2207\xb7w]+\kappa g\mu =0$ | $(v\xb7\u2207)\mu =0$ | |

$Wi[(v\xb7\u2207)\nu ]+\kappa g\nu ={0,PTT|v|2,Giesekus$ | $(v\xb7\u2207)\nu =0$ | |

$g\lambda \u223c\lambda 2|v|2,g\mu \u223c\lambda \mu |v|2,g\nu \u223c{\lambda |v|2(\nu \u2212|v|2),PTT(\mu 2+\kappa \u22121\kappa )|v|2,Giesekus$ |

## IV. UNSTEADY PLANAR 4:1 CONTRACTION FLOW

We now use the full numerical simulation of the mass and momentum conservation equations (2.6) and (2.7) in combination with the constitutive equations for CSF Eqs. (2.11)–(2.13) or NSF Eqs. (2.18)–(2.20). The numerical methodology for solving the set of equations is based on our previous work (Evans *et al.*, 2019b; 2019a; 2020). In summary, a Cartesian non-uniform mesh is employed to discretize the domain of the contraction flow geometry, with a more refined mesh close to the corners, as described in Fig. 3. The equations are discretized using finite-difference approximations, while a projection scheme is employed in order to uncouple the velocity and pressure fields in Eqs. (2.6) and (2.7). Explicit time discretizations are used for solving the constitutive equations of CSF and NSF, while the convective terms presented in these equations are approximated by an upwind methodology.

In the simulations, we have fixed some parameters as $\kappa =0.25$ and $Re=10\u22122$. The other parameters *β* and $Wi$ vary, as described in this section. Two different mesh sizes were used: M1 and M2 with the minimum size and channel lengths given in Table III.

### A. Numerical verification: Corner vortex analysis and Couette correction

Preliminary to investigating the asymptotic behaviors at the sharp corners, it is necessary to confirm first that we have complete attached flow at such corners. The asymptotic results presented in Sec. III are premised on the absence of a separating streamline at the reentrant corner, either due to the presence of a lip vortex or the extension of the salient corner vortex. This is confirmed for both models in Figs. 4 and 5, for a range of Weissenberg values with $\beta =1/2$ and 1/9, and fixed $\kappa =0.25$. No presence of lip vortices was detected (in keeping with the numerical results for 4:1 contractions in Alves *et al.* (2003; 2004), and it is clear that the larger value $\beta =1/2$ confines the salient corner vortex more than the traditional benchmark value of $\beta =1/9$. Further, it is noteworthy that our simulations have no upper Weissenberg limit in the NSF formulation for either *β* value. However, the restriction of the salient corner vortex and time required to reach steady-state makes the $\beta =1/2$ solvent viscosity fraction case optimum to investigate the reentrant corner behavior and which we exclusively use in Secs. IV B and IV C.

Our first simulation results are listed in Table IV and compare against the benchmark data of Alves *et al.* (2003) for linear PTT (no benchmark data being available for the Giesekus model). Recorded are the size of the corner vortex *X _{R}* and the Couette correction coefficient

*C*for the pressure, both of which are defined in Alves

*et al.*(2003). For consistency, we use the same minimum mesh spacing as in the benchmark case (Alves

*et al.*, 2003), this relatively coarse mesh being denoted as

*M*

_{1}, with the same model parameter values $\beta =1/9$ and $\kappa =0.25$. The CSF simulation results are in good agreement with the benchmark data for Weissenberg values up to 100, the simulations then ceasing to converge at higher Weissenberg numbers. This breakdown signifies the classical HWNP, which all numerical schemes using the CSF have historically encountered.

$Wi$ . | X (CSF)
. _{R} | X (NSF)
. _{R} | X (Alves _{R}et al., 2003)
. | C (CSF)
. | C (NSF)
. | C (Alves et al., 2003)
. |
---|---|---|---|---|---|---|

0 | 1.4922 | 1.4922 | 1.5002 | 0.3753 | 0.3753 | 0.3741 |

0.5 | 1.4922 | ⋯ | 1.5060 | 0.1696 | ⋯ | 0.1672 |

1 | 1.5309 | ⋯ | 1.5420 | 0.0980 | ⋯ | 0.0951 |

2 | 1.6112 | ⋯ | 1.6390 | 0.0284 | ⋯ | 0.0261 |

10 | 2.1799 | ⋯ | 2.1310 | −0.1543 | ⋯ | −0.1113 |

50 | 2.7240 | 1.4173 | 2.4930 | −0.1653 | −1.0352 | −0.1444 |

100 | 2.7240 | 1.4921 | 2.570 | −0.1168 | −0.7405 | −0.0917 |

1000 | ⋯ | 1.4543 | ⋯ | ⋯ | −0.1847 | ⋯ |

$Wi$ . | X (CSF)
. _{R} | X (NSF)
. _{R} | X (Alves _{R}et al., 2003)
. | C (CSF)
. | C (NSF)
. | C (Alves et al., 2003)
. |
---|---|---|---|---|---|---|

0 | 1.4922 | 1.4922 | 1.5002 | 0.3753 | 0.3753 | 0.3741 |

0.5 | 1.4922 | ⋯ | 1.5060 | 0.1696 | ⋯ | 0.1672 |

1 | 1.5309 | ⋯ | 1.5420 | 0.0980 | ⋯ | 0.0951 |

2 | 1.6112 | ⋯ | 1.6390 | 0.0284 | ⋯ | 0.0261 |

10 | 2.1799 | ⋯ | 2.1310 | −0.1543 | ⋯ | −0.1113 |

50 | 2.7240 | 1.4173 | 2.4930 | −0.1653 | −1.0352 | −0.1444 |

100 | 2.7240 | 1.4921 | 2.570 | −0.1168 | −0.7405 | −0.0917 |

1000 | ⋯ | 1.4543 | ⋯ | ⋯ | −0.1847 | ⋯ |

Focusing upon $\beta =1/2$, Table V gives the vortex size and Couette correction for both PTT and Giesekus. The CSF and NSF results are comparable for Weissenberg numbers for which both schemes converged. However, an upper Weissenberg limit was encountered in the CSF, which the NSF did not encounter for either model. In fact on the finest mesh *M*_{2}, relatively low Weissenberg numbers only were obtainable with the CSF, with the Giesekus model being more restrictive.

. | PTT . | Giesekus . | ||||||
---|---|---|---|---|---|---|---|---|

Wi . | CSF . | NSF . | CSF . | NSF . | ||||

. | X
. _{R} | C
. | X
. _{R} | C
. | X
. _{R} | C
. | X
. _{R} | C
. |

1 | 1.5059 | 0.2306 | 1.5059 | 0.2282 | 1.5059 | 0.2668 | 1.5059 | 0.2646 |

10 | 1.5682 | 0.2305 | 1.5682 | 0.2285 | ⋯ | ⋯ | 1.6330 | 0.3275 |

100 | ⋯ | ⋯ | 1.5059 | 0.2662 | ⋯ | ⋯ | 1.5059 | 0.3493 |

1000 | ⋯ | ⋯ | 1.5059 | 0.3599 | ⋯ | ⋯ | 1.5059 | 0.3641 |

10 000 | ⋯ | ⋯ | 1.5059 | 0.3640 | ⋯ | ⋯ | 1.5059 | 0.3648 |

. | PTT . | Giesekus . | ||||||
---|---|---|---|---|---|---|---|---|

Wi . | CSF . | NSF . | CSF . | NSF . | ||||

. | X
. _{R} | C
. | X
. _{R} | C
. | X
. _{R} | C
. | X
. _{R} | C
. |

1 | 1.5059 | 0.2306 | 1.5059 | 0.2282 | 1.5059 | 0.2668 | 1.5059 | 0.2646 |

10 | 1.5682 | 0.2305 | 1.5682 | 0.2285 | ⋯ | ⋯ | 1.6330 | 0.3275 |

100 | ⋯ | ⋯ | 1.5059 | 0.2662 | ⋯ | ⋯ | 1.5059 | 0.3493 |

1000 | ⋯ | ⋯ | 1.5059 | 0.3599 | ⋯ | ⋯ | 1.5059 | 0.3641 |

10 000 | ⋯ | ⋯ | 1.5059 | 0.3640 | ⋯ | ⋯ | 1.5059 | 0.3648 |

Figure 6 records the behavior of the velocity components, pressure, and Cartesian extra-stress components with radial distance from the reentrant corner along the ray $\theta =\pi /2$. In each plot, the two cases $\beta =1/9$ and $\beta =1/2$ are presented for comparison, and this is done for three selected order one Weissenberg numbers with $\kappa =0.25$ kept fixed. The theoretical value of the radial exponent for these variables is stated in each plot. The velocity components are a reasonable fit over the range of Weissenberg numbers, as are the Cartesian stress components (although these appear deteriorate as Weissenberg increases). However, the pressure behaviors do not capture the expected theoretical asymptotics.

### B. Numerical verification: Sharp corner asymptotics for $Wi=O(1)$

In this section, we present results for the parameter case $Wi=1,\beta =1/2,\kappa =0.25$ for both models. Figures 7–9 give the limiting behaviors of the velocity, pressure, Cartesian and natural stress components at the reentrant corner. Shown are the radial behaviors along selected fixed angles, uniformly spaced by $\pi /4$ around the $3\pi /2$ corner, together with the theoretical slopes of the variables from Table I. Figure 7 confirms the asymptotic Newtonian behavior of the velocity field for both CSF and NSF, with the NSF giving better convergence behavior. Less convincing is the pressure behavior, which we note is better for NSF, but still seems to require a far finer mesh for confirmation. Figure 8 gives the Cartesian stress components using the CSF, with the results indicating slower convergence and loss of resolution at the end points, than compared to the natural stress components in Fig. 9. Notable for the CSF in Fig. 8 is that the convergence rates to the limiting behaviors vary significantly with the angle of approach, with both models generally behaving similarly apart from the *T*_{22} component along the downstream ray $\theta =\pi /4$. It is also noteworthy that the convergence of the natural stress components is significantly less affected by the ray angle of approach and showing good convergence properties uniformly around the reentrant corner.

Next, we examine the dominant terms in the constitutive equations close to the reentrant corner. To enable this, we record in Table VI the radial sizes of the component terms occurring in both the CSF and NSF constitutive equations. These may be deduced directly from Table I and help order the sizes of terms in the constitutive equations in the asymptotic limit $r\u21920$ as the radial distance from the corner vanishes. First, we discuss the CSF for both models and present in Fig. 10, for each component equation in (2.11)–(2.13) the term whose modulus is a maximum. Table VI records the asymptotic sizes (with respect to the radial distance from the corner) of the terms in the CSF equations, which should hold in the core region near the corner but away from the walls. The plots do confirm the dominance of the upper convected derivative of stress for each component. However, Table VI indicates that the quadratic stress terms followed by the rate-of-strain terms are not that much smaller, and their appearance in Fig. 10 is not surprising. To see the relative ordering of the terms, we plot in Fig. 11 the modulus of the individual groups of terms in the constitutive equations along a given streamline that passes close to the reentrant corner (and shown in the plots in Fig. 10). The orientation of the polar coordinates is such that *θ* = 0 is the downstream wall and $\theta =3\pi /2\u22484.7$ is the upstream wall. The intention is to determine whether the theoretical asymptotic balances summarized in Table II hold for the boundary layers at the walls and the core region away from the walls. Figure 11 gives the sizes of the groups of terms in (2.11)–(2.13), plotted by component for the CSF. Away from the walls, it is clear that the upper convective derivative dominates in the 11- and 12-component equations, but that the quadratic stress terms are comparable in the 22-component equation. Near the upstream wall ($\theta \u22484.7$), the linear stress terms do appear smaller than the quadratic stress terms in the 11- and 12-component equations, and remain so along the rest of the streamline. The 22-component equation suggests near the upstream wall the linear stress term is comparable with the quadratic stress term, which is consistent with the balance for this component in the Giesekus equation [see Table II and Eq. (3.12)]. Figure 12 gives the behavior of the terms along the same streamline but now parameterized with the radial distance from the corner and starting from near the upstream wall. It confirms and reinforces the previous remarks.

CSF | $Tp\u25bd$ | $\kappa (1\u2212\beta )g(Tp)$ | $Tp/Wi$ | $2(1\u2212\beta )D/Wi$ |

PTT | $r\u22120.7841$ | $r\u22120.6572$ | $r\u22120.3286$ | $r\u22120.4555$ |

Giesekus | $r\u22120.7351$ | $r\u22120.5592$ | $r\u22120.2796$ | $r\u22120.4555$ |

NSF (2.18) | $v\xb7\u2207\lambda $ | $2\mu \u2207\xb7w$ | $\kappa Wig\lambda $ | $1Wi(\lambda \u22121|v|2)$ |

PTT | $r\u22121.8743$ | $r\u22121.6713$ | $r\u22121.7462$ | $r\u22121.4176$ |

Giesekus | $r\u22121.8241$ | $r\u22121.5445$ | $r\u22121.6482$ | $r\u22121.3686$ |

NSF (2.19) | $v\xb7\u2207\mu $ | $\nu \u2207\xb7w$ | $\kappa Wig\mu $ | $1Wi\mu $ |

PTT | $r\u22120.5823$ | $r\u22120.3807$ | $r\u22120.4554$ | $r\u22120.1268$ |

Giesekus | $r\u22120.4555$ | $r\u22120.1759$ | $r\u22120.2796$ | r^{0} |

NSF (2.20) | $v\xb7\u2207\nu $ | ⋯ | $\kappa Wig\nu $ | $1Wi(\nu \u2212|v|2)$ |

PTT | $r0.7083$ | $r0.7604$ | $r1.089$ | |

Giesekus | $r0.9131$ | $r1.089$ | $r1.089$ |

CSF | $Tp\u25bd$ | $\kappa (1\u2212\beta )g(Tp)$ | $Tp/Wi$ | $2(1\u2212\beta )D/Wi$ |

PTT | $r\u22120.7841$ | $r\u22120.6572$ | $r\u22120.3286$ | $r\u22120.4555$ |

Giesekus | $r\u22120.7351$ | $r\u22120.5592$ | $r\u22120.2796$ | $r\u22120.4555$ |

NSF (2.18) | $v\xb7\u2207\lambda $ | $2\mu \u2207\xb7w$ | $\kappa Wig\lambda $ | $1Wi(\lambda \u22121|v|2)$ |

PTT | $r\u22121.8743$ | $r\u22121.6713$ | $r\u22121.7462$ | $r\u22121.4176$ |

Giesekus | $r\u22121.8241$ | $r\u22121.5445$ | $r\u22121.6482$ | $r\u22121.3686$ |

NSF (2.19) | $v\xb7\u2207\mu $ | $\nu \u2207\xb7w$ | $\kappa Wig\mu $ | $1Wi\mu $ |

PTT | $r\u22120.5823$ | $r\u22120.3807$ | $r\u22120.4554$ | $r\u22120.1268$ |

Giesekus | $r\u22120.4555$ | $r\u22120.1759$ | $r\u22120.2796$ | r^{0} |

NSF (2.20) | $v\xb7\u2207\nu $ | ⋯ | $\kappa Wig\nu $ | $1Wi(\nu \u2212|v|2)$ |

PTT | $r0.7083$ | $r0.7604$ | $r1.089$ | |

Giesekus | $r0.9131$ | $r1.089$ | $r1.089$ |

Figures 13–16 repeat the previous CSF analysis but for NSF. The color map in Fig. 13 shows the dominant terms in the component equations (2.18)–(2.20). The convective term of the NS variable *λ* does dominate for both models in a region close to the corner. However, the quadratic stress terms appear to intrude very close to the corner, and this is confirmed when we follow the sizes of terms along the specified streamline close to the corner in Figs. 14 for PTT and 15 for Giesekus. In these two figures, results using the full numerical simulation code are denoted by “Full.” For comparison, we solve the NSF equations along the theoretical streamline given in (3.1) and present the same groups of terms in the plots denoted “Simplified” as a consistency check. Table VI gives the asymptotic sizes of the respective terms, and although the convective terms of the natural stresses dominate, it is seen that the quadratic stress terms are only slightly smaller. This would suggest that a much finer mesh is needed before the convective NS terms clearly dominate. It is clear that the streamline plots for PTT are supporting that the first three terms in the *λ* and *μ* NS equations and the first two terms in the *ν* NS equation in Table VI are balancing, with notably the terms $(\lambda \u22121/|v|2),\mu $ and $(\nu \u2212|v|2)$ being consistently smaller. We note that this is also the expected balance in the wall boundary layers as recorded in Table II. A similar behavior is seen in Fig. 15 for the Giesekus model, with notably the $(\nu \u2212|v|2)$ term now no longer being small and consistent with the predicted balance in Table VI (again also being the boundary layer balance in Table II). Nevertheless, despite the unclear dominance of all the convective NS variables at the reentrant corner, their behavior is still sufficiently close to the required theoretical asymptotic behaviors as evidenced by the NS variables in Fig. 9. Noticeably in Fig. 9, the *μ* variable has more difficulty converging for the Giesekus model, which may be understandable given the lack of dominance of its convective derivative in Figs. 13(d) and confirmed by Figs. 15(c) and 16(d). We note though that the behavior along the theoretical streamline in Fig. 15(d) would suggest the convected *μ* term should dominate strongly, yet its real behavior in Fig. 15(c) appears subdued. In a similar way of Fig. 12, we have described the behavior of the terms along the same streamline in Fig. 16 in order to further clear the above discussion.

### C. Numerical study at high Weissenberg numbers

In this section, we explore the effect of increasing the Weissenberg number for both models. The parameter values $\beta =1/2,\kappa =0.25$ are kept fixed, while the Weissenberg number is increased from unity to 10^{4}. The NSF was used exclusively for the simulations, due to the CSF breaking down for Weissenberg around 10^{2}. Behavior of the velocity components and pressure are given in Fig. 17 and confirm the Newtonian velocity field over the range of Weissenberg numbers. However, the pressure behavior is not confirmed. Convergence for the natural stress variables is confirmed in Fig. 18 over the range of Weissenberg numbers, and it is evident that convergence to the limiting behavior slows as the Weissenberg number increases.

Figures 19 and 20 show the behavior of the groups of terms in the NSF equations for large Weissenberg numbers. The balances in Table II again hold along the chosen streamline $\psi =3.6\xd710\u22125$ passing close to the reentrant corner. It is clear that the sizes of the terms grow extremely large at the walls, where there is now theoretically the occurrence of high Weissenberg boundary layers (Hagen and Renardy, 1997) that are present along all the walls of the domain. The thicknesses of these boundary layers are $Wi\u22121/3$ for PTT and slightly thinner $Wi\u22121/2$ for Giesekus. It is noticeable for both models that the behavior along the streamline is very different for $Wi=103$ and $Wi=104$ from $Wi=102$, where the terms involving $\u2207\xb7w$ now dominate the quadratic stress in the *λ* and *μ* component equations. This is a direct consequence of the Weissenberg number being high. The transition to this behavior occurs between 10^{2} and 10^{3}. Figure 21 illustrates that the dominance of the convected natural stress variables is more evident than the $Wi=1$ case, where small regions near the corner in the Giesekus case [panels (d) and (f) in Fig. 21] have started to emerge.

## V. DISCUSSION

The focus of this work was to verify the theoretical asymptotics for the PTT and Giesekus models at reentrant corners. This was investigated here for the $270\xb0$ corner that arises in the benchmark 4:1 contraction problem. The two models both contain quadratic stress terms that give them quite different behavior as compared to the Oldroyd-B model. While all three models share a similar asymptotic structure at a reentrant corner, their stress and boundary layer behaviors are very different. The key behavior for the PTT and Giesekus models is the dominance of the solvent stress over the polymer stress in contrast to Oldroyd-B, with thicker wall boundary layers. These features make the simulation of such models more attainable than Oldroyd-B.

The numerical simulation results of Sec. IV do verify the predicted theoretical behavior for the velocity and stress fields, which has been demonstrated over a large range of Weissenberg numbers. Not only has the theoretical stress and velocity gradient singularities of Table I been confirmed, the boundary layer balances of Table II are also seen to hold. Care was taken to ensure that a separating streamline was not present at the reentrant corner that may occur due to a lip vortex or extension of the salient corner vortex. This ensures that the reentrant corner remains sharp with the full corner angle producing a consistent stress singularity. This was one of the motivating reasons for choosing $\beta =1/2$ rather than the benchmark value of $\beta =1/9$. For the parameter values considered here of solvent viscosity fraction $\beta =1/2$ and the model-specific quadratic stress parameter $\kappa =0.25$, neither PTT nor Giesekus appears to form a lip vortex at the reentrant corner. This is noteworthy, as such a structure is present in the Oldroyd-B model and grows with the Weissenberg number (Alves *et al.*, 2003), which would make the corner singularity verification more difficult.

In the derivation of the theoretical asymptotic results in Evans (2010a; 2010b), an artificial small parameter *ϵ* based on the length scale was used as the asymptotic expansion parameter and order the sizes of terms. An estimate for its value can be deduced from Fig. 9 at around $10\u22123$ (as an upper bound) for $Wi=1$, which is the length scale on which the NSF asymptotics begin to appear. We remark that far smaller length scales are required for all components of the CSF to capture their asymptotics, as seen in Fig. 8, where convergence of some components at some angles has yet to attain their limiting behaviors. Further, the boundary layer thicknesses would be at most $3.5\xd710\u22124$ for PTT and $2\xd710\u22124$ for Giesekus using Table I, which are quite narrow to resolve, even for the finest mesh M2 used here. The color maps Figs. 10 and 13 seem to confirm this.

The simulations for PTT and Giesekus did not encounter an upper Weissenberg number value using the NSF equations. In contrast, the CSF equations broke down at Weissenberg around 10^{2}, as found also in Alves *et al.* (2003). It appears that the NSF is not only advantageous over the CSF near-stress singularities but also for global simulation of these models. Future work will examine the properties of the NSF and compare with the log conformation approach for the CSF.

One aspect of the numerical results that need improving is the pressure variable. While the velocity field behaviors are acceptable, accurate computations of the pressure variable seem more difficult to obtain as evidenced in Figs. 6, 7, and 17, and moreover suffer with increasing Weissenberg number. This may be due to the numerical treatment of the momentum equation as a diffusion equation. Certainly as the solvent fraction *β* reduces and/or the Weissenberg number increases, a different approach to the projection scheme for the velocity and pressure may be required. One feature we noted was that the pressure slopes improved dramatically for the NSF when truncated channel lengths ($8L\xd78L$) were used, as can be seen in Fig. 22. According to this figure, it is possible to see a better pressure convergence of the numerical scheme as compared with CSF for both fluid models. Therefore, high-order projection schemes for viscoelastic fluid flows could be considered for future work.

## ACKNOWLEDGMENTS

The authors would like to thank the financial support given by SPRINT/FAPESP Grant No. 2018/22242–0, The Royal Society Newton International Exchanges Grant No. 2015/NI150225, and the Núcleo de Processamento de Alto Desempenho of the Universidade Federal do Rio Grande do Norte—NPAD/UFRN to allow us to access their computer facilities. This research was developed with computational resources from Centro de Ciências Matemáticas Aplicadas à Indústria (CeMEAI) supported by FAPESP (Fundação de Amparo à Pesquisa do Estado de São Paulo) Grant No. 2013/07375–0. C. M. Oishi acknowledges CNPq (Conselho Nacional de Desenvolvimento Científico e Tecnólogico), Grant No. 307459/2016–0 and FAPESP Grant No. 2021/13833–7. F. Ruano Neto acknowledges the financial support of FAPESP Grant No. 2021/05727–2.

## AUTHOR DECLARATIONS

### Conflict of Interest

The authors have no conflicts to disclose.

### Author Contributions

**Jonathan Evans:** Conceptualization (equal); Data curation (equal); Formal analysis (equal); Investigation (equal); Methodology (equal); Project administration (equal); Supervision (equal); Validation (equal); Writing – original draft (equal); Writing – review & editing (equal). **Irineu Lopes Palhares Junior:** Data curation (equal); Formal analysis (equal); Software (equal); Validation (equal); Visualization (equal); Writing – review & editing (equal). **Cassio Oishi:** Conceptualization (equal); Formal analysis (equal); Investigation (equal); Methodology (equal); Project administration (equal); Supervision (equal); Writing – original draft (equal); Writing – review & editing (equal). **Fabiano Ruano Neto:** Data curation (equal); Investigation (equal); Software (equal); Visualization (equal); Writing – review & editing (equal).

## DATA AVAILABILITY

The data that support the findings of this study are available from the corresponding author upon reasonable request.