Computational Predictions for Boger Fluids and Circular Contraction Flow under Various Aspect Ratios ^†

Abstract

This work puts forward a modeling study contrasted against experimental, with focus on abrupt circular contraction flow of two highly-elastic constant shear-viscosity Boger fluids, i.e., a polyacrylamide dissolved in corn-syrup PAA/CS (Fluid-1) and a polyisobutylene dissolved in polybutene PIB/PB (Fluid-2), in various contraction-ratio geometries. Moreover, this work goes hand-in-hand with the counterpart matching of experimental pressure-drops observed in such 4:1 and 8:1 aspect-ratio contraction flows, as described experimentally in the literature. In this study, the experimental findings, for Boger fluids with severe strain-hardening features, reveal significant vortex-evolution characteristics, correlated with enhanced pressure-drop phasing and normal-stress response in the corner region. It is shown how such behavior may be replicated through simulation and the rheological dependencies that are necessary to bring this about. Predictive solutions with an advanced hybrid finite-element/volume (fe/fv) algorithm are able to elucidate the rheological properties (extensional viscosity and normal-stress response) that rule such vortex-enhancement evolution. This is accomplished by employing the novel swanINNFM(q) family of fluids, through the swIM model-variant, with its strong and efficient control on elongational properties.

1. Background and Introduction

Quantitative comparison between numerical predictions, experimental observations, and complex flow, occurring in contraction and contraction–expansion flows, has occupied the attention of the rheological scientific community over decades; see, for instance, Boger [1] and Boger et al. [2], López-Aguilar et al. [3], Tamaddon-Jahromi et al. [4], Nigen and Walters [5], Binding et al. [6], Pérez-Camacho et al. [7], Tamaddon-Jahromi et al. [8], López-Aguilar et al. [9], López-Aguilar et al. [10], Webster et al. [11], and reviews of Walters and Webster [12], White et al. [13], and Owens and Phillips [14]. This may be recorded in terms of pressure-drops, and vortex-activity in the recess-zones nearby salient and re-entrant corners of these geometries. In Boger [1] and Boger et al. [2], attention was given to two experimental studies with highly-elastic constant-viscosity ‘Boger’ fluids and circular contractions. In the first of these two studies [2], comparison was made between the flow of two Boger fluids, with basically the same principal characteristic relaxation-times, in three contraction-ratio geometries (α_aspect) of 2:1, 4:1, and 16:1. Findings with increasing shear-rate disclosed two distinct kinematical patterns. The first fluid under α_aspect = {4, 16}, a polyacrylamide in corn-syrup PAA/CS solution (Fluid-1, Figure 1), showed continual salient-corner vortex-growth, with separation-line adjustment in shape from concave-to-convex. In contrast, the second test polyisobutylene in polybutene PIB/PB solution (Fluid-2, Figure 2) displayed a sequential combination of salient-corner and lip vortices. Then, as deformation-rate increased, the lip kinematical structure completely engulfed the shrinking salient-corner vortex, giving way to a single large recirculating entity of convex shape (elastic-corner vortex). Note, in both experiments and at high shear-rates, the vortex extended in coverage up to the re-entrant corner. Under such different scenarios, Boger et al. [2] concluded that measurement of steady and dynamic shear properties alone were insufficient to characterize the response of such elastic liquids in circular contraction flow.

In the second and further study, Boger and Binnington [15] studied two Boger fluids to produce streak-like photographic observations for α_aspect = 4 circular contractions, of sharp and rounded-corner configurations. The first fluid was the organically-based international test fluid Ml, (polyisobutylene dissolved in polybutene (PIB/PB)). The second fluid was referred to as fluid P1 (0.03% polyacrylamide dissolved in corn-syrup (PAA/CS)). Both fluids exhibited significant elasticity, while at the same time, a constant shear viscosity. There, these two fluids exhibited distinctly different vortex-enhancement paths for a given aspect-ratio choice. For example, in the rounded-corner geometry, the salient-corner vortex appeared almost constant in size with the M1-fluid, whilst lip-vortex formation was observed under the P1 fluid case. As for the abrupt contraction, the M1-fluid displayed marginal vortex-growth, whilst vortex-enhancement was more active for the P1-fluid. Once more, these major differences between the responses in complex-flow of solutions with similar shear properties in fixed geometries render their extensional features as the subjacent explanation for such diversity of trends in vortex-activity.

In addition, Rothstein and McKinley [16] explored experimentally the creeping flow of a dilute 0.025 wt% polystyrene/polystyrene (PS/PS) Boger fluid. These authors, covering a large range of Deborah numbers, devoted attention to circular contraction–expansion flow-settings of various aspect-ratios (α_aspect = {2, 4, 8}) and re-entrant corner curvature degrees. For a relatively low aspect-ratio of α_aspect = 2 of sharp-corners, a steady lip-vortex was observed at the re-entrant corner. For aspect-ratios of between α_aspect = 4 and 8, lip-vortex formation was absent, but a salient-corner vortex was recorded, which grew with the flow-rate increase. Rounding the re-entrant corner shifted such landmarks and trends to higher values of flow-rates, but did not change qualitatively the structure and evolution of the overall flow-field.

Sato and Richardson [17] performed simulations for planar α_aspect = 4 contraction flow. These were based on a hybrid finite volume/element method, embedded in a time-stepping procedure within a pressure-correction scheme. These authors reported lip-vortex formation as a pseudo-transient phenomenon, appearing at Reynolds-number levels of Re = 0.01, and being triggered by an instantaneous increase in Deborah numbers (De) from six to twelve. Subsequently, such transient lip-vortex faded through the time-stepping process, as a steady-state solution was approached at the limiting value of elasticity of De = 12. Similarly, for the same α_aspect = 4 planar contraction flow, Olsson [18] also observed the transient presence of a lip-vortex, but using the Giesekus rheological equation-of-state; whilst employing a method-of-lines technique for time-integration and a discretization based on finite-differences.

In keeping with the above developments, the present study considers counterpart predictive solutions generated with a hybrid-subcell finite-element/volume algorithm (fe/fv) [19,20,21], incorporating some novel advanced stabilization techniques [22,23]. Attention is directed towards contraction-ratios of α_aspect = {2, 4, 8}, covering in particular the correlation of pressure-drop enhancement, vortex-dynamics (lip-vortex formation), and flow-structure (normal-stress response); the context is one of Boger fluids and creeping flow conditions. This range of contraction-ratios was held sufficient for present comparison purposes, as gathered from our prior work on contraction-expansion ratio comparison in López-Aguilar et al. [9]. An appeal is also made to our companion study in Tamaddon-Jahromi et al. [4], where the focus of attention there was solely on the α_aspect = 8 contraction-ratio problem.

2. Governing Flow Equations, Material Functions, Problem Specification, and Numerical Algorithm

Following the principles of conservation of momentum and mass, the non-dimensional equations that govern the flow response of viscoelastic fluids under creeping incompressible and isothermal conditions, are:

$$\nabla \cdot \mathit{u}=0$$

(1)

$$Re\frac{\partial \mathit{u}}{\partial t}=\nabla \cdot \mathit{T}-Re\mathit{u}\cdot \nabla \mathit{u}-\nabla p.$$

(2)

As such, a domain bounded in space and time (x, t) is considered, over which spatial–temporal differential operators apply. Then, field variables $\mathit{u},p$, and T represent fluid velocity, hydrodynamic pressure, and stress-tensor, respectively. The stress-tensor may be expressed as:

$$\mathit{T}=\mathit{\tau }+2\beta \mathit{d}.$$

(3)

The stress-tensor is decomposed into two parts by means of the Elastico-Viscous Stress Splitting (EVSS) assumption, where the total stress T is composed by two contributions, one for the polymer, to which the viscoelastic nature is addressed through τ, and another of Newtonian-like response of the form $2\beta \mathit{d}$. In this, $\mathit{d}=(\nabla \mathit{u}+\nabla {\mathit{u}}^{†})/2$ represents the rate-of-deformation tensor, where tensor transpose is denoted with the superscript $†$. In addition, the non-dimensional group Reynolds number is defined as $Re=\frac{\rho U_{char}{L}_{char}}{{\mu }_0}$, through characteristic scales of $U_{char}$ on fluid velocity (mean flow-velocity over the characteristic-length), and, for length, $L_{char}$, as the constriction radius. The material density is represented with $\rho$ and the characteristic viscosity taken as a zero shear-rate viscosity (${\mu }_0={\mu }_p+{\mu }_s$). Here, ${\mu }_p$ and ${\mu }_s$ are the polymeric viscosity and the solvent viscosity components, respectively, so that the solvent-fraction parameter can be defined as $\beta =\frac{{\mu }_s}{{\mu }_0}$. Creeping flow conditions are presumed throughout, so that Reynolds numbers are typically O(10⁻²) or smaller.

2.1. Constitutive Modeling

To complete the equation set, one needs a state law on stress, which is provided by the swanINNFM(q) model formalism (see Tamaddon-Jahromi et al. [4,8] and López-Aguilar et al. [9]). This model is soundly-based, being derived from two well-respected models, the Finitely Extensible Non-linear Elastic dumbbell Chilcott-Rallison FENE-CR model (Chilcott and Rallison [24]) and a White–Metzner model (White and Metzner [25]).

The base FENE-CR model may be written in a configuration-tensor $\mathit{A}$ form (bold-face symbols denote tensorial quantities), as:

$$Wi\stackrel{\nabla }{\mathit{A}}+f\left[tr\left(\mathit{A}\right)\right](\mathit{A}-\mathit{I})=0.$$

(4)

Here, $\stackrel{\nabla }{\mathit{A}}$ stands for the upper-convected material-derivative of the configuration-tensor ($\mathit{A}$), defined as:

$$\stackrel{\nabla }{\mathit{A}}=\frac{\partial \mathit{A}}{\partial t}+\mathit{u}\cdot \nabla \mathit{A}-{(\nabla \mathit{u})}^{†}\cdot \mathit{A}-\mathit{A}\cdot (\nabla \mathit{u}).$$

(5)

The FENE-CR structural-functional $f\left[tr\left(\mathit{A}\right)\right]$ is:

$$f\left[tr\left(\mathit{A}\right)\right]=\frac{1}{1-\frac{tr\left(\mathit{A}\right)}{L^2}}.$$

(6)

Then, Kramer’s rule relates configuration and extra-stress tensors as follows:

$$\mathit{\tau }=\frac{1-\beta }{Wi}f\left[tr\left(\mathit{A}\right)\right]\left(\mathit{A}-\mathit{I}\right).$$

(7)

In the above, $L$ is the extensibility parameter for the FENE-CR model, related to the dumbbell chain-length, and $\mathit{I}$ is the identity tensor. In addition, the non-dimensional Weissenberg group-number is defined as $Wi={\lambda }_1\frac{U_{char}}{L_{char}}$, where dependency upon the fluid relaxation-time (${\lambda }_1$) and a characteristic rate ($\frac{U_{char}}{L_{char}}$) is observed. Then, rise in Wi may be generated through deformation-rate increase, fixing the fluid elastic character through ${\lambda }_1$, whilst the flow-rate Q-dynamics are increased, i.e., $Wi={\lambda }_1\frac{Q}{\pi L_{char}^3}$, considering that $Q=A{U}_{char}=\pi L_{char}^2{U}_{char}$.

To arrive at the swanINNFM(q) model, one needs to consider a rate-dependent viscosity in the above developments, as under the generalized White–Metzner model, taking this to be extension-rate dependent alone. Then, using Equations (3) and (7) above, the new resulting family of swanINNFM(q) models, in its single relaxation-time swIM model-variant [3,4,8,9,10,11], may be articulated through the amended total-stress tensor, as:

$$\mathit{T}=\frac{1-\beta }{Wi}f\left[tr\left(\mathit{A}\right)\right]\left(\mathit{A}-\mathit{I}\right)\varphi \left(\dot{\epsilon }\right)+2\beta \varphi \left(\dot{\epsilon }\right)\mathit{d},$$

(8)

where the dissipative extensional-function $\varphi \left(\dot{\epsilon }\right)$ is taken as a quadratic-form from the truncated Taylor-series approximation of the cosh-exponential expression available. This dissipative extensional-function is defined as $\varphi \left(\dot{\epsilon }\right)=1+{\left({\lambda }_D\dot{\epsilon }\right)}^2$, with parameterization on a dissipative material time-scale parameter ${\lambda }_D$, and functionality on a generalized strain-rate invariant $\dot{\epsilon }=\frac{II{I}_d}{I{I}_d}$. Here, $II{I}_d$ and $I{I}_d$ represent the third and the second invariants of $\mathit{d}$, respectively. Fuller details on the development of this swanINNFM(q)-family of fluids are supplied in Debbaut and Crochet [26], Debbaut et al. [27], Tamaddon-Jahromi et al. [4,8], López-Aguilar et al. [3,9,10], and Webster et al. [11].

2.2. Material Functions

The relevant swIM-model rheometrical functions are provided in Figure 3, noting a constant shear-viscosity. These are the extensional viscosity ${\eta }_e$ and the first normal-stress difference in shear $N_{1Shear}$, where variation over model parameters ($\beta$, $L$, ${\lambda }_D$) is presented. Their functional forms are, respectively:

$${\eta }_e=3\varphi (\dot{\epsilon })\beta +3\varphi (\dot{\epsilon })\left(1-\beta \right)\left[\frac{f^2}{f^2-fWi\dot{\epsilon }-2W{i}^2{\dot{\epsilon }}^2}\right],$$

(9)

$$N_{1Shear}=\frac{2(1-\beta )Wi{\dot{\gamma }}^2}{f}.$$

(10)

In Figure 3a, a solvent-content $\beta$-variation extensional viscosity ${\eta }_e$-response is exposed for the swIM model, under ${\lambda }_D$ = 0.075, $L$ = 5, and $\beta$ = {0.9, 0.8, 0.5, 1/9}. Firstly, the swIM model-response appears bounded by the two extremes of its behavior; under ${\lambda }_D$ = 0, the dissipative extensional influence disappears, and the swIM model extensional hardening finds a plateau for moderate-to-high extension-rates; on the other extreme, at $L$ → ∞, an Oldroyd-B-like response is recovered, with infinite extensional viscosity predictions. Within these two bounds, swIM $\beta$-decrease, which may be interpreted as an increase in solute-content, renders a rise in the plateau-level observed at intermediate shear-rates in the range of 1 ≤ ${\lambda }_1\dot{\epsilon }$ ≤ 70 units. Beyond such plateaued-stage, a steep rise is witnessed as a result of the influence of the extensionally-driven dissipative mechanism promoted by $\varphi \left(\dot{\epsilon }\right)=1+{\left({\lambda }_D\dot{\epsilon }\right)}^2$; here, $\beta$-decrease shifts this ${\eta }_e$-rise to higher extension-rates. In Figure 3b, the first normal-stress in shear $N_{1Shear}$-response is plotted for both Oldroyd-B and swIM under $\beta$-decrease. Here, in contrast to the stiff quadratic rising Oldroyd-B $N_{1Shear}$-trend, swIM provides a softer trend for shear-rates beyond ${\lambda }_1\dot{\gamma }$~10 units; in addition, $\beta$-decrease shifts $N_{1Shear}$-rise to lower shear-rates. Such predictive capabilities of the swIM model are contrasted against $N_{1Shear}$ data reported by Boger et al. [2] for both Boger fluids formed by diluted solutions of PAA/CS and PIB/PB. This data is presented in dimensionless form, taking as characteristics time and viscosity scales, the characteristic time and the constant shear-viscosity reported experimentally by Boger et al. [2] as {λ₁, η₀} = {0.380 s, 97.5 P} for the PAA/CS solutions and {λ₁, η₀} = {0.149 s, 251 P} for the PIB/PB case. Stark matching is recorded between the experimental rheometrical $N_{1Shear}$ and the predictions achieved using the swIM model. Particular to Figure 3b, a window of experimental-data capture is defined in the ranges of solvent-fraction 0.5 ≤ β ≤ 0.9, extensibility-parameter 5 ≤ L ≤ 12, and dissipative-parameter 0 ≤ λ_D ≤ 0.1, which is used in subsequent sections for the simulation of contraction complex flow of those Boger fluids.

In Figure 3c, swIM extensional viscosity response with $L$-variation is reported. Particularly, this parametric study is performed under $\beta$ = {0.9, 0.7}, ${\lambda }_D$ = {0, 0.075}, and $L$ = {5, 8, 10, 12}. swIM extensional viscosity under $L$-increase is analogous to that observed under solvent-fraction $\beta$-decrease, displaying an intermediate plateaued region, and followed by a sharp increase. Interestingly, the cumulative response of increasing both solvent-fraction $\beta$, the extensibility-parameter $L$ and the extensional-dissipative time-scale ${\lambda }_D$ of $\beta$ = 0.7, $L$ = 12 and ${\lambda }_D$ = 0.3, exposes the strength of this model to boost hardening in extensional viscosity. The effects of such variations in first normal-stress in shear $N_{1Shear}$ are, under $L$-increase, to enhance elasticity beyond ${\lambda }_1\dot{\gamma }$ ~ 10 units (Figure 3d). Notably, the coincidence between swIM predictions and $N_{1Shear}$ experimental data-trends reported by Boger et al. [2] holds.

In Figure 3e, the influence of ${\lambda }_D$-variation over extensional viscosity is provided under $\beta$ = 0.9 and $L$ = 5. Here, the departure from the FENE-CR trend at intermediate extensional-rates appears sooner as the level of ${\lambda }_D$ is larger, even vanishing the plateaued section for ${\lambda }_D$ = 0.4 in comparison with smaller ${\lambda }_D$-cases. In terms of first normal-stress in shear $N_{1Shear}$ (Figure 3f), as devised for this swIM model, ${\lambda }_D$-increase does not affect response in shear deformations.

2.3. Problem Specification and Numerical Scheme

The meshes used to discretize the problem, described in a 2D-domain, under three aspect-ratios α_aspect = {2, 4, 8}, are displayed in Figure 4. Their characteristics in terms of node number and degrees-of-freedom are provided in

Table 1. Mesh characteristics, different contraction aspect-ratios.

Mesh αaspect	Elements	Nodes	$\mathbf{Degrees}\mathbf{of}\mathbf{Freedom}\left(\mathit{u},\mathit{p},\mathit{\tau }\right)$
2	2762	5787	36,235
4	2987	6220	38,937
8 (refined)	2016	4191	26,234
8 (medium)	1707	3634	22,768
8 (coarse)	868	1897	11,897

. On mesh-refinement and solution-convergence, one may refer to the counterpart study in Tamaddon-Jahromi et al. [4], where this topic is well-covered. One may note the additional fine meshing used in the α_aspect = 4 instance, around the re-entrant corner zone, which is necessary to pursue a stringent lip-vortex search.

Boundary conditions. On flow boundary-conditions, the flow specification is as follows. No-slip is taken on boundary walls, and shear-free symmetry is imposed on the flow-centerline (see also below). At flow-inlet, velocity and stress are specified, according to the flow-rate setting, akin to full-developed shear-flow. There, under vanishing convective terms, the partial differential equations PDE for stress-configuration collapses to a temporal ordinary differential equation ODE system, providing evolution to the algebraic stress equivalent forms. At outlet, only the arbitrary level of pressure is set. Then, through a steady-state solution-continuation procedure, initial conditions from a prior flow-rate solution may be accessed, and under more severe-flow parameter selection, feed-forward exit-procedures may be used on velocity-gradients and stress components to accelerate convergence (see López-Aguilar et al. [9,22]). Such feed-forward procedure overwrites fully-developed polymeric-stress $\mathit{\tau }$ and velocity-gradient $\nabla \mathit{u}$ components from inter-field regions towards the outlet-edge neighborhood prior solution approximation in each time-step. This helps to reduce noise proliferation originated at the outlet, which reflects back towards the internal field, and is particularly useful under moderate-to-high flow-rates (López-Aguilar et al. [9,22]). Generally, a flow-rate increase mode is adopted through a series of steady-state solutions, as appropriate and as prescribed elsewhere (Tamaddon-Jahromi et al. [4,8], López-Aguilar et al. [3,9,10]).

On numerical-to-experimental scaling. In the α_aspect = 4.08 contraction flow of Boger [1] and Boger et al. [2], a Weissenberg number definition was introduced as ${Wi}^{exp}={\lambda }^{exp}\dot{\gamma }$; where ${\lambda }^{exp}$ is a relaxation-time and $\dot{\gamma }$ is a characteristic shear-rate (downstream wall shear-rate in tubular entry flow). Accordingly, respective relaxation-times were identified of ${\lambda }^{exp}$ = 0.149 s for the polyisobutylene (PIB) in polybutene (PB) Boger fluid, and ${\lambda }^{exp}$ = 0.308 s for the polyacrylamide (PAA) in water and corn syrup (CS) Boger solution. A functional relationship may be derived between these two experimental and computational Wi-definitions; experimentally of (${Wi}^{exp}={\lambda }^{exp}\dot{\gamma }$) and computationally of ($W{i}^{comp}={\lambda }_1\frac{U_{char}}{L_{char}}$). For one-to-one comparison purposes, this establishes appropriate scaling factors of ($\frac{1}{{\lambda }^{exp}}$), between experimental findings and computational predictions per fluid-instance; yielding: $W{i}_{}^{comp}=6.71W{i}_{}^{exp}$ for fluid PIB/PB and $W{i}_{}^{comp}=3.25W{i}_{}^{exp}$ for fluid PAA/CS (see López-Aguilar et al. [3] and Tamaddon-Jahromi et al. [4] for more details on such scaling). Based on these scaling factors, $W{i}_{}^{exp}$ as in Boger [1] compares as 10.5 ≤ ${Wi}^{comp}$ ≤ 16.6 under fluid PIB/PB, and 2.05 ≤ ${Wi}^{comp}$ ≤ 5.3 under fluid PAA/CS. In practice, one notes below for the α_aspect = 4, that computationally, a slight lip-vortex appears for 2.5 ≤ ${Wi}^{comp}$ ≤ 5.5 (lower rate range); whilst experimentally, a lip-vortex is only observed for fluid PIB/PB (in the higher ${Wi}^{comp}$-range) and not recorded with fluid PAA/CS.

Hybrid subcell finite-element/finite-volume scheme. The numerical method used in this work is based on a hybridized finite-element (fe) and finite-volume (fv) spatial-discretization scheme. Such a formulation comprises both time-stepping and fractional-staged (three) equation-structure. On the momentum-mass conservation equation doublet, finite-element (fe) discretization, grafted upon a Taylor–Petrov–Galerkin structure, is selected following incremental pressure-correction strategy. On the constitutive stress-equation, finite-volume (fv) discretization is employed. Such a space-time discretization agrees with equation-type specification. Accordingly, Galerkin-type (fe) approximation is applied on parent triangular tesselations; whilst a subtended subcell/cell-vertex finite-volume (fv) discretization is used for the rheological equation-of-state on stress-tensor components. An element-by-element iterative solution-procedure, space-efficient in its implementation, is utilized for discretized equations. The pressure-equation is solved with a direct Choleski-reduction method. Then, the (fv)-component on stress is treated in a direct single-iteration implementation. The conservation-form for the stress equation is non-linear, and contains inhomogeneous source terms. This demands both median-dual-cell treatment for source terms and fluctuation-distribution for fluxes (upwinding). Additionally, quadratic interpolation is chosen for velocity, whilst linear interpolation is specified for pressure on the parent fe triangular-cell grid. For the finite-volume implementation, four fv-subcells per parent fe-cell are obtained, being the fv-sub-cell constructed via the interconnection of the parent fe-cells mid-side nodes. In such a structured arrangement, stress variables are located at the vertices of fv-sub-cells, and solution projection between is unnecessary. On the child subcell-level, this provides for a subcell-vertex fv-method for which trial-solutions are interpolated linearly. The resulting formulation is consistent in time and holds an accuracy of second-order. Further details on the numerical scheme and its detailed implementation characteristics, can be found in Wapperom and Webster [19], Webster et al. [20], and Aboubacar and Webster [21].

Stabilization techniques. Additional and latest aspects of improved stabilization techniques for viscoelastic flow employed are summarized as follows. The set of such stabilization techniques comprises the velocity-gradient VGR-correction [22,23], the use of velocity-gradient recovery, a discrete continuity correction over the flow-domain, and additional compatibilizing conditions on the flow-centerline, with pure-extension shear-free inhomogeneous extensional deformation at the centerline symmetry. Additionally, the absolute-value ABS-f correction [22,23] regularizes the problem through absolute-value imposition on the structure-network functional (f) in the rheological equation-of-state and within the Kramer’s rule transformation in the momentum equation. Positive-definiteness of the problem is promoted through the use of configuration-tensor form in the constitutive equation [3,4,8,9,11,22,23].

3. Results—Computational Predictions and Flow-Structure versus Pressure-Drop Correlation

Comparison of predictive solutions is presented across the three geometric aspect-ratios, in turn, of α_aspect = {2, 4, 8}. In this, it is informative to consider ramping-up through flow-rate (low, medium, high), where, due to the variation in dynamics per geometry, these ranges themselves will vary per geometry. Specific insight is drawn through comparison against the counterpart experimental patterns of Boger [1] and Boger et al. [2], and particularly, when focusing on Boger fluids of two different solvent-fractions of $\beta$ = 0.87 and $\beta$ = 0.73, as extracted above in $N_{1Shear}$-match with swIM model in Section 2.2. The Results section is organized in two main subsections. Firstly, 3.1 8:1 contraction flow includes: flow-structure (vortices and first normal-stress difference N₁) and pressure drops, renders the main findings of this work, and evidences the matching of experimental pressure-drops using the swIM model. Here, conspicuously, vortex-development phasing is correlated with pressure-drop enhancement, as described by Binding and Walters [28], providing theoretical explanation to experimental features of this benchmark circular contraction problem. Moreover, viscoelastic response in the recess-zones, observed through first normal-stress N₁-fields, appears directly linked with the vortex-formation and evolution, with salient-corner, lip, and elastic corner vortex capture, where the tracking of the shape and size of vortex-structures is recorded. This lies as a major finding on relating pressure-drop enhancement, vortex-evolution (flow kinematics), and flow structure [9]. The second section, Section 3.2 Predictive capabilities of the swIM—Vortex-dynamics across α_aspect = {2, 4, 8} circular contraction flow, provides insight into the influential swIM model parameters that permit the prediction of the elusive lip-vortex. Particular attention is paid to swIM solvent-fraction β, extensibility-parameter L and dissipative parameter λ_D variation, for which exploration of their vortex-dynamics and lip-vortex-formation are explored, and where intermediate λ_D-values at relatively high extensibility L-features appear as a proper combination for lip-vortex capture. Such parametrical-specification correlates with the precise control of extensional properties provided by the swanINNFM(q) model-family, embodied here through its swIM variant.

3.1. 8:1 Contraction Flow: Flow-Structure (Vortices and First Normal-Stress Difference N₁) and Pressure Drops

Notably and overall, all geometric aspect-ratios explored display elastic-corner vortices (ecv) at large flow-rates, as it is apparent in Figure 1 and Figure 2. Hence, to describe the evolution of such kinematic structures and its relationship with pressure-drops with flow-rate rise, one may begin with the α_aspect = 8 aspect-ratio case (see Tamaddon-Jahromi et al. [4]), as this instance provides the strongest dynamics and the sharpest distinction in flow-pattern features and viscoelastic pressure-drop rise arising experimentally. Such an α_aspect = 8 geometry is reflective of higher ratios and reveals the evolving streamline patterns and exaggerated pressure-drop trends of Figure 5. First, recorded in relatively low flow-rate range of Q/Q₀ ≤ 0.035 units, salient-corner vortices (scv) at low-rates arise, accompanied by pressure-drops that concur with Newtonian equivalents (here, Q₀ and Δp₀ are, respectively, characteristic flow-rate and pressure-drop taken from experiments [1]). Then, with flow-rate increase, in the low-to-mid flow-rate range of 0.12 ≤ Q/Q₀ ≤ 0.2 units, such relatively simple behavior gives way to co-existent salient-corner/lip vortices (lv), instance for which viscoelastic pressure-drop initiates its departure from Newtonian response. Further flow-rate increase drives coalescence of the co-existent scv and lv, and marks the entry to the onset of large elastic-corner vortices (ecv) in the mid-to-high flow-rate regime of Q/Q₀ ≥ 0.35 units, in which stark departure in pressure-drop between simple Newtonian and viscoelastic Boger fluids is apparent. Such vortex-evolution, from salient-corner vortex to elastic-corner vortex, has been proposed by Binding and Walters [28] as a cause for the pressure-drop enhancement observed experimentally.

Correlation of complex flow features and Boger-fluid normal-stress difference in complex flow. It was established in López-Aguilar et al. [9] that the various vortex-structures and flow-stages (scv, lv, and ecv) tie in closely with the corner-patterns sustained in N₁-fields from complex flow, whilst vortex-evolution with Q-increase across such scv, lv, and ecv stages are closely driven the extensional viscosity response in ideal extensional deformation, as predicted with the swIM-model. The present study would concur with this, as evidenced in Figure 6, where sample streamlines are contrasted against N₁-fields for α_aspect = 8 case. Note that in the circular contraction complex flow at hand, the first normal-stress difference is defined as N₁ = τ_zz−τ_rr, where τ_zz and τ_rr represent the normal-stresses in z and r directions, respectively. Here, the phases of scv growth, co-existence of scv-lv, and ecv-domination are mirrored in the N₁ field-data of Figure 6. Under this evidence, one may conclude that elasticity and, hence, non-linearity, manifest through first normal-stress difference in shear and extension (extensional viscosity), is observed to strongly influence the formation of such vortex-structures and its counterpart energetic pressure-drop effects in departure from Newtonian equivalent levels.

3.2. Lip-Vortex Predictive Capabilities of the SwIM—Vortex-Dynamics across α_aspect = {2, 4, 8} Circular Contraction Flow

Once that prediction of experimental pressure-drops and their corresponding vortex-structure from Boger [1] and Boger et al. [2], have been achieved using the swIM model for an α_aspect = 8, a parametric study is worthy to explore the predictive capabilities of this swanINNFM(q) family of models, in terms of the elusive task of predicting the appearance and persistence of lip-vortex structures in less demanding α_aspect = {2, 4} circular contraction geometries. Here, one may note that the versatile swIM model, with its control of normal-stress response, both in elongational (extensional-viscosity) and shear first normal-stress in shear deformations, is able to answer such a question. Specifically, different flow-structures and their evolution are observed by varying parameters of solvent-fraction $\beta$, extensibility $L$-parameter, and extensional dissipation-parameter ${\lambda }_D$. Each of these material parameters may be associated with corresponding variations in rheological properties (see Figure 3), observed through extensional viscosity η_e and normal-stress difference in shear N_1Shear. Under alternative parameter-variation, one notes the selection of implied base-values of: $L$ = 5 for $\beta$–variation, and $\beta$ = 0.9 for $L$-variation. Then, the dissipation-parameter may be set either at ${\lambda }_D$ = 0 or ${\lambda }_D$ = 0.075. From this perspective, counterpart findings may then be explored comparatively for flows in lower and alternative geometric contraction-ratios. Here, parameters influencing normal-stress response in both shear and extensional deformations appear with a key-role in the development of varied flow-structure formation in the contraction flows analyzed.

3.2.1. 8:1 Contraction Flow: Flow-Rate and Solvent-Fraction Adjustment (1/9 ≤ β ≤ 0.9)

This analysis starts with the α_aspect = 8 circular contraction case, with focus on three flow-rate regimes that highlight essential vortex-structure features.

Low flow-rates. Wi = [1, 2] Under swIM [$L$ = 5, ${\lambda }_D$ = 0.075], only salient-corner vortex (scv) activity is recorded at Wi = 1 (Figure 7a), with no apparent lip-vortex (lv) activity. Then whilst remaining at Wi = 1, vortex-size and intensity strengthen considerably with solvent-fraction $\beta$-reduction. With increase in flow-rate level, at Wi = 2 (Figure 7b), the first appearance for lv-formation is found at $\beta$ = 0.7, attendant with scv-presence. At lower solvent-fractions, such lip-vortex presence considerably enhances, to even take over the salient-corner vortex-intensity and produce a single strong elastic-corner vortex with a prominently convex separation-line; see Boger [1] and Boger et al. [2] for analogous experimental trends. Consistently, one observes ecv Ψ_min-intensity nearly doubling from $\beta$ = 0.5 to $\beta$ = 0.3; and tripling from $\beta$ = 0.5 to $\beta$ = 1/9, as recorded in Figure 7.

Medium Wi = [3, 5] and high [Wi = 16] flow-rates. Addressing the medium flow-rate regime, and comparing against that at low flow-rates, lip-vortex formation appears earlier with solvent-fraction $\beta$-decline. So, for example, at Wi = 3, first lip-vortex detection is noted with more dilute systems at $\beta$ = 0.9, see Figure 8a. Flow-patterns and trends outlined at lower flow-rates are then repeated under solvent-fraction $\beta$-decrease. This precedes pronouncedly concave ecv-formation by $\beta$ ≤ 0.6. This theme is continued into the second yet higher-rate level in the medium flow-rate regime of Wi = 5, shown in Figure 8b. In the high flow-rate regime, still larger ecv Ψ_min-intensities abound in Figure 9, through some 8-times increase in vortex-size from $\beta$ = 0.9 to $\beta$ = 1/9 and dramatic vortex separation-line adjustment.

Extensional dissipative-parameter${\lambda }_D$-variation. In passing, it is worth reflecting on the influence of dissipation extensional-influence and ${\lambda }_D$-rise on the results above; specifically in the medium flow-rate range Wi = 3, adopting other parameters of $L$ = 5, $\beta$ = 0.9. In Figure 10, this is illustrated through streamline-patterns across the range 0.1 ≤ ${\lambda }_D$ ≤ 1.0, in a regime of strong lip-vortex activity, both in intensity and spatial occupation. This is useful for insight upon the α_aspect = 4 ratio case reviewed below. In this data, one can detect a clear amplification of the lip-vortex in the range 0.1 ≤ ${\lambda }_D$ ≤ 0.4, prior to coalescence with the salient-corner vortex (0.5 ≤ ${\lambda }_D$), the lip-vortex being the dominant feature both prior to and post-coalescence. As a consequence, the eye of the vortex-center, subsequent to coalescence, is driven towards the re-entrant corner (response observed experimentally with Boger Fluids 1 and 2; Boger [1]), announcing the onset of a phase of strong elastic corner-vortex domination, and further, unsteady vortex-oscillation at still larger flow-rates [8].

3.2.2. Comparison across Geometric Aspect-Ratios α_aspect

Predictive solutions with swIM [$L$ = 5, ${\lambda }_D$ = 0.1, $\beta$ = 0.9] provide direct comparison on major findings across the various contraction-ratios, taken against increasing flow-rate (Q) and charted through rise in Wi, according to low, moderate, and high flow-rate regimes. Note that to truly amplify the detail, scaling is applied in the cross-stream direction. Accordingly, streamline field-patterns are reported in Figure 11 and the line-graph of Figure 12. This data conveys the corresponding trend in movement on vortex-intensity. Notably, lip-vortex appearance is only recorded in the larger contraction-ratio α_aspect = 8.

Low flow-rate (Wi = 1). In the first row of Figure 11, comparable vortex-structures are discerned across all three aspect ratios, gathering common concave-shaped separation-lines (referenced to the salient-corner recess), and with practically identical scv Ψ_min-intensity of O(10⁻⁴).

Intermediate flow-rate regime [Wi = 3, Wi = 5]. Here, more interesting distinction can be drawn. At (Wi = 3), the second row of Figure 11, both α_aspect = {2, 4} solutions retain scv Ψ_min-intensity of O(4 × 10⁻³), yet proving one order-of-magnitude larger rotational intensity than at corresponding Wi = 1. With (α_aspect = 8, Wi = 3), there is a relatively marginal decline in scv-intensity to O(3 × 10⁻³) noted, whilst also supporting some energy transfer into the onset of a lip-vortex, of intensity one order-of-magnitude lower, i.e., O(5 × 10⁻⁴). As above, here vortex separation-lines retain their concave-shaped form. At the more dynamic level of Wi = 5, corresponding to the third row of Figure 11, for α_aspect = 8, there is a sudden burst of activity, with an increase of two orders-of-magnitude in lip-vortex rotational intensity, from Wi = 3 (row-two) to Wi = 5 (row-three). Simultaneously, the standing scv-intensity triples to O(9 × 10⁻³). Now for the first time, conspicuously, the vortex separation-line begins to adjust in shape around the lip-vortex zone, depicting somewhat of a more convex-to-concave delineation. In contrast, neither of the lower ratio α_aspect = {2, 4} solutions pick up any sign of lip-vortex activity, whilst their scv-intensities reflect levels comparable to those of {lip-vortex, α_aspect = 8}, being slightly larger at O(2 × 10⁻²).

High flow-rate regime of Wi = [10, 20]. Rows four and five of Figure 11 now display vortex separation-lines of convex shapes for all three ratio-solutions. Nevertheless, the α_aspect = 8 Wi = 10-solution (row-four of Figure 11) is disparate, in that its elastic-corner vortex has an evolution history that passes through coexistent lv-scv structures. In contrast, both α_aspect = {2, 4} solutions provide an ecv delivered from a growing salient-corner vortex directly, without any intermediate transition. The relative position of vortex-centre loci across contraction-ratio, clearly displays dependency upon their evolution history through flow-rate rise, with α_aspect = 8 locating the ecv closer to the re-entrant corner. These vortex-evolution patterns concur well with Boger Fluid-1 (PAA/CS) findings of Boger [1]. Particularly notable is the trend observed for aspect-ratios (α_aspect ≥ 4), and the movement of the vortex-eye gravitating towards the re-entrant corner with increasing flow-rate (Figure 11).

Informative overall trends can be gathered from vortex-intensity (−Ψ_min) data listed in

Table 2. Vortex-intensity (−Ψ_min); α_aspect = {2, 4, 8}; Wi = [1, 20]; swIM [$L$ = 5, $\beta$ ≤ 0.9]; ${\lambda }_D$ = 0.1.

Wi/αaspect	2	4	8
1	0.0006	0.0009	0.0008
3	0.0043	0.0048	0.0005 0.0029(lv)
5	0.0244	0.0277	0.0093 0.0271(lv)
10	0.0841	0.2874	0.3028
20	0.2654	1.132	1.589

and its graphical representation in Figure 12. From a united trend at low flow-rates, there is a pronounced rise in vortex-intensity at (Wi ≥ 5), in instances α_aspect = {4, 8} above α_aspect = 2. By Wi = 20, separation is clearly apparent between all three instances, with the largest ecv-intensity attracted by α_aspect = 8 (with lv-formation marked); and trends in α_aspect = 4 solutions follow closely those under α_aspect = 8. Notably, trends in α_aspect = 2 ratio solutions are the least dynamic in adjustment, providing smooth and continuous rise in scv-intensity.

3.2.3. α_aspect = 4 and α_aspect = 2 Ratios: Lip-Vortices, Rise in Wi, Extensibility-Parameter $L$, and Solvent-Fraction $\beta$ Switch

In this section, a parametric study on swIM extensibility parameter $L$ and solvent fraction $\beta$ is carried out, to discern the possibility of lip-vortex formation in less stringent contraction-ratios of α_aspect = {2, 4}. Recall that Boger [1] and Boger et al. [2] only observed a lip-vortex experimentally for the PIB/PB-based fluid at {${Wi}^{exp}={\lambda }^{exp}\dot{\gamma }=2.3$; ${Wi}_{{\lambda }_1}^{comp}=6.71{Wi}_{\lambda }^{exp}=15.4$} under the α_aspect = 4.08 ratio flow. This is performed through manipulation of extensional properties modulated by $L$-variation and solute-content promoted by $\beta$-decrease. Particularly, one turns to predictions for fluids with more pronounced extensional features and slightly higher solvent-fractions, i.e., $L$ = 12, ${\lambda }_D$ = 0.1, and $\beta$ = 0.7 (as noted in Boger Fluids 1 and 2 of Boger [1]). One should mention that relatively more diluted fluids ($\beta$ = 0.8; Figure 13) give signs of lip-vortex formation, but such kinematical structures remain difficult to track under such conditions.

Lip-vortex-capture under α_aspect = 4, $\beta$ = 0.7, $L$ = 12 and λ_D = 0.1 settings. Figure 14 illustrates a successful lip-vortex prediction under α_aspect = 4. Firstly, an order-of-magnitude increase is detected in lv-intensity from Wi = 4 (ψ_min = −0.0005) to Wi = 5.7 (ψ_min = −0.0053). This proves to be the largest lip-vortex observed for the various different α_aspect = 4 trial-settings. As such, Wi = 5.7 represents a critical level, beyond which steady-state solution-tractability fails. What is apparent is the fine balance in rheology here, between both solvent-fraction settings and hardening-levels. This interplay clearly has strong impact upon such localized issues as lip-vortex appearance (or not, as the case may be). Clearly experimentally, both such outcomes were observed in Boger [1]; with Boger Fluid-2 substantiating a bulb-like lip-vortex.

Seeking lip-vortices through${\lambda }_D$-rise under$\beta$= {0.8, 0.7} and$L$= 12. Following on from the parameter adjustment and findings on vortex activity in Tamaddon-Jahromi et al. [4] under α_aspect = 8 flow, one may be peaked to further investigate the distinct influence of the extensional dissipation-parameter ${\lambda }_D$. In particular, seeking the segregated impact of extensional viscosity alone on α_aspect = 4 lip-vortex response, as identified above. Figure 15 provides further evidence supplied with ${\lambda }_D$-rise, at the associated two solvent-fraction levels of $\beta$ = {0.8, 0.7}. From this data, it is clear that early dissipation-factor ${\lambda }_D$-rise (with its η_e-boosting control) does strengthen lip-vortex activity, prior to this being subsumed by the more dominant salient-corner vortex, as the latter fingers its way towards the re-entrant corner. This response is simply intensified with reduction in solvent-fraction, with $\beta$ = 0.7 subsuming lv-activity at ${\lambda }_D$ = 0.4, whilst $\beta$ = 0.8 performs likewise by ${\lambda }_D$ = 0.3.

On α_aspect = 2 flow. Finally, the flow dynamics for the α_aspect = 2 contraction-ratio are so mild, that one needs to locate large polymeric-composition as high as $\beta$ = 1/9 to detect any indication of lv-activity. Even then, such activity is limited to the temporal evolution phase alone, and vanishes prior to a steady-state being established (also noted by others for planar α_aspect = 4 contraction flow, see Sato and Richardson [17], Olsson [18]). In this instance, the influence of extensibility structure $L$-parameter adjustment on vortex-structures is barely noticeable.

4. Conclusions

Findings in this study go some way to replicating the experimental observations of Boger [1] and Boger et al. [2], on flow structures generated in three circular contraction flows of different aspect-ratios of α_aspect = {2, 4, 8}, and with two different Boger fluids, of essentially the same rheological response. Predictions reveal the delicate rheological balances at play, to match such response, taking into account the separate and combined influences of extensional viscosity and first normal-stress difference. Accordingly, elusive lip-vortices have been captured, as well as transition phases between salient and lip-vortices, their coalescence and the domination of elastic-corner vortices, and their relationship pressure-drops and first normal-stress difference in complex combined shear-to-extensional flow. Predictions permit fine interrogation of such phenomena. The flow history through flow-rate Q-rise provides the key to recognizing the dominant rheological characteristics here, where the swIM-model’s precise tuning of normal-stress response on both shear and extension permits the capture of lip-vortices in the three geometrical configurations explored. Moreover, a link between flow-structure (through vortex-morphology and first normal-stress difference in the vicinities of the constriction) and pressure-drop enhancement staging is reported. In this respect, vortex-shape concurs well with N₁-field structure in the recess zones, for which connection between kinematics and rheological response is exposed. The pressure-drop enhancement mechanism, associated with the non-linear viscoelastic features of the Boger solutions, correlates well with vortex-dynamics phasing, where transition is recorded from early concave salient-corner vortices, coexistent lip and salient-corner vortices, and evolution to convex elastic-corner vortices, as reported experimentally [28] and numerically [9].

Of particular note across contraction-ratio α_aspect = {2, 4, 8} solutions is the gradually shifting upwards of the vortex-centre loci. It is attracted to the re-entrant corner, yet the relative positions of vortex-centre loci clearly display dependencies upon their salient-corner to elastic-corner vortex evolution history through flow-rate rise. These vortex evolution patterns concur well with the Boger Fluid-1 (PAA/CS) and Fluid-2 (PIB/PB) findings of Boger [1], prominent for aspect-ratios (α_aspect ≥ 4), and the movement of the vortex-eye gravitating towards the re-entrant corner with increasing flow-rate, finally leading to elastic-corner vortex domination and viscoelastic pressure-drop enhancement.