Numerical study of viscoelastic upstream instability

Abstract

In this work, we report numerical results on the flow instability and bifurcation of a viscoelastic fluid in the upstream region of a cylinder in a confined narrow channel. Two-dimensional direct numerical simulations based on the FENE-P model (the finite-extensible nonlinear elastic model with the Peterlin closure) are conducted with numerical stabilization techniques. Our results show that the macroscopic viscoelastic constitutive relation can capture the viscoelastic upstream instability reported in previous experiments for low-Reynolds-number flows. The numerical simulations reveal that the non-dimensional recirculation length (L_D) is affected by the cylinder blockage ratio (BR), the Weissenberg number (Wi), the viscosity ratio (β) and the maximum polymer extension (L). Close to the onset of upstream recirculation, L_D with Wi satisfy Landau-type quartic potential under certain parameter space. The bifurcation may exhibit subcritical behaviour depending on the values of L² and β. The parameters β and L² have nonlinear influence on the upstream recirculation length. This work contributes to our theoretical understanding of this new instability mechanism in viscoelastic wake flows.

1. Introduction

A mixture of a Newtonian fluid and high-molecular-weight polymers (even at extremely low concentration) exhibits viscoelasticity. In the low-Reynolds-number range, when the elasticity becomes the dominant source of nonlinearity, this flow can manifest many interesting and novel flow phenomena that have not been studied extensively, especially flow instabilities, such as symmetry breaking (Arratia et al. 2006; Poole, Alves & Oliveira 2007; Haward, Toda-Peters & Shen 2018; Haward, Hopkins & Shen 2020), secondary flow (Yue, Dooley & Feng 2008; Davoodi et al. 2018), time dependency and even elastic instability and turbulence (Shaqfeh 1996; Groisman & Steinberg 2000, 2001; Larson 2000; Schiamberg et al. 2006; Grilli, Vázquez-Quesada & Ellero 2013; Varshney & Steinberg 2019; Steinberg 2021). These instability phenomena widely occur in processes ranging from plastic manufacturing (Denn 2001; Varchanis et al. 2021) to human blood plasma (Brust et al. 2013; Thiébaud et al. 2014) and porous media flows (Walkama, Waisbord & Guasto 2020; Hopkins, Haward & Shen 2021). However, detailed understandings of the underlying mechanisms in many instances remain vague. In these flows, because the extra elastic stress plays a key role, it is essential to know its value and how it influences the flow field. Experimental measurement and numerical simulation are the two main approaches to obtain the elastic stress. Compared to a large number of reported experimental studies with regards to elastic instability phenomena, the numerical method is applied relatively less (Poole 2019) to this problem and deserves to be employed more frequently because of its advantages such as more feasible parametric studies and direct analysis of the elastic stress field. However, a notoriously difficult problem in numerical simulation of viscoelastic flow has challenged researchers for decades, namely the numerical instability when simulating flow at high Weissenberg number (Wi), which is called the high-Weissenberg-number problem (HWNP). The reader is referred to the recent review paper by Alves, Oliveira & Pinho (2021) for a more complete account of these issues. In the following, we focus on explaining an idealized flow model, that is, viscoelastic flow around a circular cylinder, and the recent studies of its upstream instability.

Viscoelastic flows past bluff bodies are frequently encountered in industrial applications (e.g. filtration processes and oil extraction) and natural phenomena (e.g. flow in soil and blood flow in cardiovascular valves and brain tissues) (Larson 1999; Iliff et al. 2012; Marsden 2014). In addition, porous media are frequently modelled by ordered and disordered arrays of microfluidic circular cylinders (Walkama et al. 2020; Hopkins et al. 2021). Placing one or more cylinders in channel or pipe inlet section is often used as a disturbance source to study the channel viscoelastic instability in experiments (Varshney & Steinberg 2018; Qin et al. 2019b). Viscoelastic flow around a circular cylinder has also been regarded as a benchmark case for numerical and experimental studies (Ultman & Denn 1971; Dhahir &Walters 1989). In Newtonian wake flows, a downstream recirculation zone and vortex shedding can be observed because of the global instability in these flows (Williamson & Rosirko 1988; Barkley 2006; Sipp & Lebedev 2007; Tang et al. 2020). However, in viscoelastic wake flows, interestingly, an upstream recirculation can form in cylinder wake flows confined between two plane plates.

The recirculation in front of a circular cylinder was firstly reported by Kenney et al. (2013). Later, Shi et al. (2015), Zhao, Shen & Haward (2016), Qin et al. (2019a), Haward et al. (2021) and Hopkins et al. (2022a,b) successively studied this phenomenon experimentally. These experiments were performed using Boger solutions, which are elastic without shear thinning, such as dilute polyethylene oxide or polyacrylamide solutions, or using wormlike fluids with a strong shear-thinning behaviour, such as cationic surfactant cetyltrimethylammonium bromide and stable hydrotropic salt 3-hydroxynaphthalene-2-carboxylate solutions. Regarding the geometry of the experimental set-up, the blockage ratio (BR), i.e. the ratio of cylinder diameter to channel width, in these experiments is not less than 50 % and the depth-to-diameter ratio (α) widely ranges from 0.5 to 5. Zhao et al. (2016) summarized the flow patterns in the Wi–BR space with five states: Newtonian-like state, bending streamlines, vortex growth upstream, unsteady downstream and three-dimensional time-dependent chaotic upstream, which are shown in figure 1. Qin et al. (2019a) found that this flow is inherently three-dimensional and observed symmetry breaking as well as strong upstream propagation effects via elastic waves. At low BR (~10 %), there also exists a mild upstream instability, with only streamline bending, which was reported by Ribeiro et al. (2014), Nolan et al. (2016) and Haward et al. (2018). Therefore, it seems that high BR is one of the conditions for (strong) upstream instability. Besides, the experiments by Shi & Christopher (2016) and Varshney & Steinberg (2017, 2018) also showed that when the cylinders are arranged in tandem at low BR, the recirculation zone is likely to appear in the upstream region of the rear cylinders. Few studies quantitatively investigated the length of upstream recirculation. A schematic diagram of the upstream recirculation length (l) is shown in figure 2, which is the distance from the foremost end of the upstream recirculation to the foremost end of the cylinder. Using the cylinder diameter D to normalize l as L_D = l/D, Zhao et al. (2016) obtained a Landau-type behaviour of the dimensionless upstream recirculation length L_D (L_D increasing with Wi). In the experiment of Qin et al. (2019a), L_D is almost linear with the Weissenberg number based on the first normal stress difference $W{i_N} = {N_1}/2\dot{\gamma }\eta (\dot{\gamma })$, where N ₁ is the first normal stress difference, $\dot{\gamma }$ is the characteristic shear rate defined by the mean centreline velocity and $\eta (\dot{\gamma })$ is the shear viscosity at this $\dot{\gamma }$. It is worth pointing out that the definitions of the Weissenberg number in these two studies are different.

Figure 1. Summary for the Wi–BR state space of flow patterns. Panels (a, b, d), (c) and (e) are adapted from Zhao et al. (2016), Hopkins et al. (2021) and Qin et al. (2019a), respectively.

Figure 2. Schematic diagram of upstream recirculation length.

In the simulations of viscoelastic fluid flow, an additional elastic stress divergence term is often linearly introduced to the right-hand side of the Navier–Stokes equations. The elastic stress is closed by the constitutive relation with conformation tensor, such as upper-convected Maxwell model (Olsson & Yström 1993), Oldroyd-B model (Oldroyd 1950), finite extensible nonlinear elastic (FENE) model series (Herrchen & Öttinger 1997), Phan-Thien and Tanner model (Thien & Tanner 1977), etc. The calculated Wi is often limited to a low value due to the aforementioned numerical instability. Varchanis et al. (2020) simulated the asymmetric viscoelastic flow around a cylinder, which agreed well with their experiment. Their numerical simulation was still limited to a low BR.

Flow around a cylinder in a channel is a typical mixed flow, including shear flow and tensile flow. Using full-field time-resolved flow-induced birefringence imaging, Zhao et al. (2016a) found that upstream instabilities are associated with high stress in the fluid that accelerates in the narrow gap between the cylinder surface and the channel wall when BR is large. In the narrowest area between the cylinder surface and the channel wall, the fluid parcels are strongly stretched, which results in high elastic stress there. This sharp growth in elastic stress may cause loss of convergence and trigger the well-known HWNP in numerical simulation (Hulsen, Fattal & Kupferman 2005). It is a challenging task to balance numerical accuracy and stability. Up to now, there have been no systematic numerical simulations to investigate this upstream flow phenomenon (Poole 2019). In this paper, we present numerical simulations of the viscoelastic upstream instability recently observed experimentally in front of a cylinder in a narrow channel. The FENE model with the Peterlin closure (FENE-P) is adopted to describe the rheological constitutive behaviour of a dilute polymer solution. The square root reconstruction method (Balci et al. 2011) is selected as the stabilization technique for the numerical simulation to handle HWNPs.

The rest of the paper is organized as follows. Section 2 introduces our specific problem, governing equations, solution method, mesh generation and grid convergence test. In § 3, we discuss our numerical results on the upstream instability. We conclude the present work in § 4. In Appendix A, we also provide more results on the validation, some additional results on the FENE-CR model, the effect of Péclet number and the effect of different boundary conditions.

2. Problem formulation and numerical method

2.1. Problem description

Figure 3(a) shows the computational domain, which consists of a pressure-driven flow from left to right past a circular cylinder in a confined channel. It is true that Qin et al. (2019a) revealed the inherent three-dimensional nature of the upstream recirculation zone. However, we study the corresponding flow instability in this two-dimensional setting as the first attempt. The width of the channel is H. The diameter of the cylinder is D. Three different blockage ratios (BR = D/H), i.e. 50 %, 62.5 % and 75 %, are considered. The distance between the cylinder centre and the inlet is 37.5H (upstream region l ₁) and the distance between the cylinder centre and the exit is 25H (downstream region l ₂). A parabolic flow profile $u(y) = 1.5\bar{U}(1 - 4{y^2}/{H^2})$ is adopted as the inlet velocity, where $\bar{U}$ is the average velocity. The zero-pressure and zero-velocity gradient conditions are adopted at the outlet. On the upper wall, the lower wall and the cylindrical wall, the no-slip boundary condition is imposed. The treatment of boundary conditions for the conformation tensor is discussed below and also in Appendix A.4.

Figure 3. (a) Schematic diagram of computational domain and (b) the mesh topology near the cylinder.

2.2. Governing equations

The viscoelastic fluids in the experiments of Qin et al. (2019a,b) and Pan et al. (2013) show both the shear thickening of elongational viscosity and the shear thinning of shear viscosity. Thus, the FENE-P model is adopted in the present study, which takes into account the finite elongation of polymer molecules and the bounded stress. The governing equations of the flow combined with the FENE-P constitutive model read as follows (Bird, Dotso & Johnson 1980; Bird, Armstrong & Hassager 1987):

(2.1)\begin{equation}\left. {\begin{array}{*{20}{c}} {\boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{u} = 0,}\\ {\rho \dfrac{{\partial \boldsymbol{u}}}{{\partial t}} + \rho (\boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{\nabla })\boldsymbol{u} ={-} \boldsymbol{\nabla }p + {\eta_s}\Delta \boldsymbol{u} + \boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{\tau },}\\ {\dfrac{{\partial \boldsymbol{c}}}{{\partial t}} + (\boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{\nabla })\boldsymbol{c} - (\boldsymbol{\nabla }\boldsymbol{u})\boldsymbol{\cdot }\boldsymbol{c} - \boldsymbol{c}\boldsymbol{\cdot }{{(\boldsymbol{\nabla }\boldsymbol{u})}^T} ={-} \dfrac{\boldsymbol{\tau }}{{{\eta_p}}},} \end{array}} \right\}\end{equation}

with

(2.2a–c)\begin{equation}\boldsymbol{\tau } = \frac{{{\eta _p}}}{\lambda }[f(\boldsymbol{c})\boldsymbol{c} - \theta \boldsymbol{I}],\quad f(\boldsymbol{c}) = \frac{1}{{1 - tr(\boldsymbol{c})/{L^2}}}\quad \textrm{and}\quad \theta = \frac{{{L^2}}}{{{L^2} - 3}},\end{equation}

where ρ, u, p, t, λ, c, τ and I are fluid density, flow velocity, pressure, time, statistical relaxation time of polymers, conformation tensor, stress tensor of polymers and identity tensor, respectively, $\eta _p$ and $\eta _s$ are the polymeric contribution to zero-shear-rate viscosity and solvent viscosity, respectively, tr denotes trace operator of tensor and the extensibility parameter L measures the maximum stretching of the polymer chains (Bird et al. 1980, 1987; Purnode & Crochet 1998). When L is set as ∞ and θ tends to 1, the FENE-P model returns to the Oldroyd-B model (Oldroyd 1950).

Note that another common FENE-type model is the FENE-CR model, which has been widely applied to predict the flow behaviour of viscoelastic fluids with a constant shear viscosity and a bounded elongational viscosity. We have also performed numerical simulations based on the FENE-CR model. More information about the FENE-CR model and a comparison of the numerical results obtained using the FENE-P model and the FENE-CR model are provided in Appendix A.3.

The main dimensionless parameters for this simulation are the Reynolds number Re and the Weissenberg number Wi, which are defined as $Re = \rho \bar{U}H/2{\eta _0}$ and $Wi = 2\lambda \bar{U}/H$, respectively. It is emphasized that the present study uses the half-channel width as the reference length. Thus, a factor of 2 should be applied to the corresponding Wi defined using the whole channel width as the reference length when comparing with published data. Here $\eta$₀ = $\eta_{p}$ + $\eta_{s}$ is the total viscosity of the solution at zero shear rate and β = $\eta_{s}$/$\eta_{0}$ is the solvent viscosity ratio, a measurement of polymer concentration and molecular characteristics of polymers. In this study, β is set to 0.59 in most of the cases. Besides, β = 0.9, 0.75, 0.45, 0.3 and 0.15 are also considered in some cases. Four values of L ² = 400, 2500, 10 000 and 40 000 are considered. It is suggested from the definitions of Re and Wi that the fluid viscoelasticity is more prominent than the fluid inertia for microscale flows. Fluid viscoelasticity can be characterized by Wi. In this work, we mostly only change Wi to probe how the elasticity influences the upstream instability.

The drag coefficient acting on the cylinder is computed as

(2.3)\begin{equation}{C_d} = \frac{1}{{{\eta _0}\bar{U}H}}\int_S {[ - p\boldsymbol{I} + \boldsymbol{\tau } + {\eta _s}\boldsymbol{\nabla }\boldsymbol{u}]\boldsymbol{\cdot }\boldsymbol{i}\boldsymbol{\cdot }\textrm{d}\boldsymbol{S}} ,\end{equation}

where S is a vector normal to each face element of the cylinder boundary, whose magnitude is equal to the area of face element, and i is a unit vector aligned with the streamwise direction.

The fluctuation frequency (f_f) of flow field and elastic stress wave in the x direction can be obtained by fast Fourier transform (FFT) of the time series of C_d. The Strouhal number St is defined as

(2.4)\begin{equation}St = {f_f}H/(2\bar{U}).\end{equation}

In the experiment of Qin et al. (2019a) and Pan et al. (2013), the Weissenberg number is defined using the first normal stress difference. In order to compare our numerical results with their experimental results, we calculate the first normal stress difference of the FENE-P model under simple shear flow as follows (Purnode & Crochet 1998):

(2.5) \begin{align}{N_1} = \frac{{(\textrm{1} - \beta ){\eta _0}{L^4}}}{{\lambda ({L^2} - 3)}}{\left[ \vphantom{{\Delta _1^{2/3} + \Delta _2^{2/3} - 2{{\left( {\frac{{2{L^6}}}{{27{\lambda^2}{{\dot{\gamma }}^2}{{({L^2} - 3)}^2}}}} \right)}^{1/3}}}} {\frac{{2{L^6}}}{{{\lambda^2}{{\dot{\gamma }}^2}{{({L^2} - 3)}^2}}}} \vphantom{{\Delta _1^{2/3} + \Delta _2^{2/3} - 2{{\left( {\frac{{2{L^6}}}{{27{\lambda^2}{{\dot{\gamma }}^2}{{({L^2} - 3)}^2}}}} \right)}^{1/3}}}}\right]^{ - 1/3}}\left[ {\Delta _1^{2/3} + \Delta _2^{2/3} - 2{{\left( {\frac{{2{L^6}}}{{27{\lambda^2}{{\dot{\gamma }}^2}{{({L^2} - 3)}^2}}}} \right)}^{1/3}}} \right],\end{align}

where

(2.6a,b)\begin{equation}{\Delta _1}= 1 + \sqrt {1 + \frac{{2{L^6}}}{{27{\lambda ^2}{{\dot{\gamma }}^2}{{({L^2} - 3)}^2}}}} ,\quad {\Delta _2} = 1 - \sqrt {1 + \frac{{2{L^6}}}{{27{\lambda ^2}{{\dot{\gamma }}^2}{{({L^2} - 3)}^2}}}} .\end{equation}

The shear viscosity for the special strain rate $(\dot{\gamma })$ reads

(2.7)\begin{equation}\eta = \beta {\eta _0} + (\textrm{1} - \beta ){\eta _0}{\left[ {\frac{{{L^6}}}{{4{\lambda^2}{{\dot{\gamma }}^2}{{({L^2} - 3)}^2}}}} \right]^{1/3}}[\Delta _1^{1/3} + \Delta _2^{1/3}].\end{equation}

In experiment, the material properties of a viscoelastic fluid may be complicated and not easy to fit the parameters in the existing viscoelastic constitutive models, including the existence of branched chains of polymer molecules, the spectrum distribution of relaxation time and the certain distribution of chain lengths of polymer molecules. Herein, we use a normal stress similarity to re-calibrate the rheological parameters. The Weissenberg number based on the first normal stress difference $W{i^{\prime}_N}$ can then be defined as the ratio of the first normal stress difference to twice the Newtonian shear stress (Qin et al. 2019a,b), i.e.

(2.8)\begin{equation}W{i^{\prime}_N} = {N_1}/[2\dot{\gamma }(\eta - {\eta _s})],\end{equation}

which is different from Wi for the FENE-P model. The $W{i^{\prime}_N}- Wi$ relationship for β = 0.59 is plotted in figure 4(a). For the Oldroyd-B model, L tends to infinity and we have

(2.9)\begin{equation}W{i^{\prime}_N} = Wi.\end{equation}

Figure 4. (a) The $W{i^{\prime}_N}- Wi$ relationship and (b) the Wi_N – Wi relationship for the FENE-P model with β = 0.59.

It is worth mentioning that the above definitions for Wi and $W{i^{\prime}_N}$ do not include the effect of β. Thus, we adopt another definition here proposed in Pan et al. (2013):

(2.10)\begin{equation}W{i_N} = {N_1}/(2\dot{\gamma }\eta ).\end{equation}

In particular, for the Oldroyd-B model with $L \to \infty $, one has

(2.11)\begin{equation}W{i_N} = (1 - \beta )Wi.\end{equation}

This definition reflects the influence of β. The value of Wi_N becomes smaller for larger β. The corresponding Wi_N–Wi relationship for β = 0.59 is shown in figure 4(b). Note that for β = 0.59 and L ² = 2500, Wi_N is equal to 9.391 when Wi = 100, which is close to the maximum Wi_N in the experiments of Qin et al. (2019a,b) and Pan et al. (2013).

2.3. Numerical method

The governing equations are solved by the open-source CFD platform OpenFOAM (Weller et al. 1998) and the rheoTool toolbox (Pimenta & Alves 2018). Because no numerical studies have reported the upstream instability in viscoelastic wake flows, we detail our numerical method below and provide extensive validation of the numerical code. In order to ensure the boundedness of tr(c) = c_kk in the FENE-P model, an implicit algorithm is used for pre-calculation before each time step (Richter, Iaccarino & Shaqfeh 2010). Tracing the transport equation of the conformation tensor yields

(2.12)\begin{equation}\frac{{\partial {c_{kk}}}}{{\partial t}} + (\boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{\nabla }){c_{kk}} = tr[(\boldsymbol{\nabla }\boldsymbol{u})\boldsymbol{\cdot }\boldsymbol{c} + \boldsymbol{c}\boldsymbol{\cdot }{(\boldsymbol{\nabla }\boldsymbol{u})^T}] - \frac{{({c_{kk}} - 3){L^4}}}{{\lambda ({L^2} - {c_{kk}})({L^2} - 3)}}.\end{equation}

By defining

(2.13)\begin{equation}\varphi ={-} \ln \left( {1 - \frac{{{c_{kk}}}}{{{L^2}}}} \right),\end{equation}

equation (2.12) can be rewritten as

(2.14)\begin{equation}\frac{{\partial \varphi }}{{\partial t}} + {u_j}\frac{{\partial \varphi }}{{\partial {x_j}}} = \frac{{{e^\varphi }}}{{{L^2}}}\left( {{c_{kj}}\frac{{\partial {u_k}}}{{\partial {x_j}}} + {c_{jk}}\frac{{\partial {u_j}}}{{\partial {x_k}}}} \right) + \frac{{{e^\varphi }}}{{\lambda ({L^2} - 3)}}(3{e^\varphi } + {L^2} - {L^2}{e^\varphi }).\end{equation}

The scalar φ is solved at each time step, and then saved for the calculation of the next time step. Once φ is obtained at a given time step, the conformation tensor c can be calculated by (2.13).

Because of the high singularity of the constitutive governing equations, a small global artificial dissipation term $\kappa \varDelta \boldsymbol{c}$ is added to the right-hand side of the transport equation of conformation tensor c in order to avoid divergence. In our simulations, we set Re = 0.0001. Molecular dissipation of polymer should be considered at this Re. The presence of this additional diffusive term can be justified by the diffusivity of polymer in solvent, which was estimated to be over the range of 10⁻⁵ to 10⁻⁷ cm² s⁻¹ (Haggerty, Sugarman & Prud'homme 1988). The Schmidt number Sc = $\eta_{0}$/ρκ, which is defined as the ratio between the zero-shear-rate viscosity and the polymer molecular diffusivity, is introduced to quantify the artificial dissipation. For dilute solution of polymer dissolved in water, the density ρ is around 1000 kg m⁻³. The viscosity may vary over a wide range. For example, $\eta_{0}$ is equal to 0.3 Pa s in the experiment of Qin et al. (2019a). Thus, Sc is estimated over the range of 10⁵ to 10⁸. The normalized transport equation of conformation tensor features two dimensionless numbers, Wi and the Péclet number (Pe = Re × Sc). The value of Pe, instead of Sc, will be specified in our discussions. In the present study, Re is set to be 0.0001. The corresponding range of Pe is from 10 to 1000. The effect of Pe on the flow behaviour is discussed in Appendix A.2. By comparing the numerical results and experimental results of Qin et al. (2019a), Pe is fixed at 40 for the final simulation.

In order to improve the stability of numerical calculation, a stabilization technique must be adopted. More information on common stabilization techniques can be found in Appendix A.1. We perform a preliminary study to test two stabilization techniques, i.e. the logarithmic reconstruction method and the square root reconstruction method. In the numerical results based on the logarithmic reconstruction method, we do not observe the upstream recirculation found in the experiments of Qin et al. (2019a). Moreover, the logarithmic reconstruction method was also adopted by Mokhtari et al. (2022) and Kumar & Ardekani (2022) and the upstream recirculation was not reported in their simulations. However, the square root reconstruction method can successfully predict the upstream recirculation, which is, therefore, adopted in the present simulation. A new symmetric tensor b is introduced satisfying $\boldsymbol{c} = \boldsymbol{b}\boldsymbol{\cdot }{\boldsymbol{b}^T}$. The transport equation of conformation tensor c can then be rewritten as

(2.15)\begin{equation}\frac{{\partial \boldsymbol{b}}}{{\partial t}} + (\boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{\nabla })\boldsymbol{b} = \boldsymbol{b}\boldsymbol{\nabla }\boldsymbol{u} + \boldsymbol{ab} + \frac{1}{{2\lambda }}(\theta \boldsymbol{\cdot }{({\boldsymbol{b}^T})^{ - 1}} - \boldsymbol{b}{e^\varphi }) + \frac{\kappa }{2}\Delta \boldsymbol{b} + \kappa \boldsymbol{h},\end{equation}

where

(2.16)\begin{equation}\boldsymbol{h} = {\boldsymbol{b}^{ - 1}}[{\textstyle{1 \over 2}}(\Delta \boldsymbol{b})\boldsymbol{b} + {({\partial _x}\boldsymbol{b})^2} + {({\partial _y}\boldsymbol{b})^2} + {({\partial _z}\boldsymbol{b})^2}].\end{equation}

As advised by Balci et al. (2011), the term $\kappa \boldsymbol{h}$ in (2.15) is ignored in our numerical study. Note that a in (2.15) is an antisymmetric tensor, which can be written in the form of components as

(2.17)\begin{equation}\boldsymbol{a}= \left( {\begin{array}{*{20}{@{}c@{}}} 0&{{a_{12}}}&{{a_{13}}}\\ { - {a_{12}}}&0&{{a_{23}}}\\ { - {a_{13}}}&{ - {a_{23}}}&0 \end{array}} \right).\end{equation}

The components of a can be calculated by solving the following equations:

(2.18)\begin{equation}\left. {\begin{array}{*{20}{c@{}}} {({b_{11}} + {b_{22}}){a_{12}} + {b_{23}}{a_{13}} - {b_{31}}{a_{23}} = {w_1},}\\ {{b_{23}}{a_{12}} + ({b_{11}} + {b_{33}}){a_{13}} + {b_{12}}{a_{23}} = {w_2},}\\ { - {b_{13}}{a_{12}} + {b_{12}}{a_{13}} + ({b_{22}} + {b_{33}}){a_{23}} = {w_3},} \end{array}} \right\}\end{equation}

where

(2.19)\begin{equation}\left. {\begin{array}{*{20}{c@{}}} {{w_1} = ({b_{12}}{u_{1,1}} - {b_{11}}{u_{2,1}}) + ({b_{22}}{u_{1,2}} - {b_{12}}{u_{2,2}}) + ({b_{23}}{u_{1,3}} - {b_{13}}{u_{2,3}}),}\\ {{w_2}= ({b_{13}}{u_{1,1}} - {b_{11}}{u_{3,1}}) + ({b_{33}}{u_{1,3}} - {b_{13}}{u_{3,3}}) + ({b_{23}}{u_{1,2}} - {b_{12}}{u_{3,2}}),}\\ {{w_3} = ({b_{13}}{u_{2,1}} - {b_{12}}{u_{3,1}}) + ({b_{23}}{u_{2,2}} - {b_{22}}{u_{3,2}}) + ({b_{33}}{u_{2,3}} - {b_{23}}{u_{3,3}}).} \end{array}} \right\}\end{equation}

Here u_{i ,j} are the components of $\boldsymbol{\nabla }\boldsymbol{u}$. For a detailed description of this method, the reader is referred to Balci et al. (2011).

The second-order backward scheme is used for time discretization. The time step Δt in our simulation is set small enough to ensure numerical stability. The MINMOD scheme is adopted for discretization of $\boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{\nabla }\boldsymbol{b}$. The second-order upwind scheme is adopted for discretization of $\boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{\nabla }\boldsymbol{u}$. For the tensor b, the linear extrapolation boundary with second-order accuracy is imposed on the upper and lower walls, and no-flux boundary condition with first-order accuracy on the cylindrical wall, to ensure the stability of numerical calculation. More discussion on the boundary condition for the tensor b can be found in Appendix A.4. The pUcoupled algorithm is adopted for the pressure–velocity coupling (Jareteg 2012; Pimenta & Alves 2019). Our simulation is performed using the rheoFoam solver module of rheoTool in OpenFOAM extend 4.0 (Pimenta & Alves 2018). The validation of the present numerical method is provided in Appendix A.1.

The subcritical bifurcation can be distinguished by checking whether the simulation result is path-dependent or not, as shown in figure 5. In the present study, we systematically vary the controlling parameter Wi and the initial condition to examine their effects on the final simulation result. In the increasing Wi process, we gradually and slowly increase Wi and use the numerical result of the current state to initialize the simulation for the next simulation. In the decreasing Wi process, we gradually and slowly decrease Wi and use the numerical result of the current state to initialize the simulation for the next state. We explicitly point out which process is applied to perform the simulation only if the results at the same Wi obtained for the increasing and decreasing Wi processes are not identical. For an instability indicator A around the linear critical condition (Wi_c, as shown in figure 5), finite-amplitude solutions can exist even if Wi is less than Wi_c, which results in a different increasing Wi path from a decreasing Wi path (the hysteresis phenomenon). This method for identifying subcritical bifurcation has been widely used in previous studies (Becherer, Morozov & van Saarloo 2009; Pan et al. 2013; Burshtein et al. 2017). In this paper, L_D, the root-mean-square upstream recirculation length RMS(L_D), the time-averaged drag coefficient $\overline {{C_d}} $, the root-mean-square drag coefficient C_drms and an asymmetry parameter I (defined in (3.8)) are selected as the instability indicators.

Figure 5. Method to identify subcritical bifurcation. The solid and dashed lines denote the stable and unstable solutions, respectively. Parameter Wi_c denotes the critical condition; IW and DW denote the increasing and decreasing Wi processes, respectively.

2.4. Mesh generation and grid convergence test

A block-structured mesh is generated for the computational domain using the commercial software ANSYS ICEM. The surrounding region of the cylinder is discretized by an O-type mesh, as shown in figure 3(b). The remainder of the computational domain is discretized using several blocks of quadrilateral meshes. In the y direction, Ny = 121 (for BR = 50 %, 201 for BR = 65 % and 75 %) grid points are unevenly distributed. In the x direction, Nx = Ny. The dense mesh is applied near the cylinder. The O-type mesh comprises Ns = 2(Nx + Ny − 2) = 800 (BR = 65 % and 75 %) or 440 (BR = 50 %) grid points uniformly distributed along the cylinder perimeter and Nr = 200 (BR = 65 % and 75 %) or 120 (BR = 50 %) grid points stretched over an exponential progression along the radial direction to ensure a fine mesh near the cylinder surface. We set the first cell side-by-side to the cylinder surface in the radial direction to 0.00125H. In the x direction, Nl ₁ = 701 grid points are set in the upstream region, and Nl ₂ = 301 grid points are unevenly arranged in the downstream region. In the y direction, Ny = 121 (201 for BR = 65 % and 75 %) grid points are unevenly distributed. The total number of cells in the computational domain is approximately 280 000 (BR = 65 % and 75 %) or 150 800 (BR = 50 %). The detailed mesh generation description can be found in our previous publication (Peng et al. 2021).

The parameter set (BR, Wi, L ², β) = (50 %, 30, 2500, 0.59) is selected to check the mesh dependency and the corresponding simulation is performed on three sets of meshes, which are listed in table 1. The effect of the mesh size on L_D is investigated. In our numerical simulation, L_D is measured as the horizontal distance from the most remote upstream stagnation point (u = 0 and v = 0) to the front end of the cylinder. The simulation results show that L_D for Mesh1 is time dependent and its variation range is recorded in table 1. Length L_D is steady for Mesh2 and Mesh3 and the difference between the results for Mesh2 and Mesh3 is 3 %. We thus consider that our numerical results for Mesh2 are accurate enough and the final simulations are performed on this mesh.

Table 1. Parameter set (BR, Re, Wi, L ², β) = (50 %, 0.0001, 30, 2500, 0.59). Here $N{l_1}$ and $N{l_2}$ denote the grid numbers of the upstream and downstream regions along the x direction, respectively; Nx and Ny denote the grid numbers in the x and y directions, respectively; Nr denotes the grid numbers radiating from the cylinder surface; and $\Delta r$ is the grid size of the innermost grid near the cylinder wall.

3. Results and discussion

In this section, we report the numerical results of the viscoelastic wake flow in a confined channel. The parameters considered are (BR, β) = (50 %, 0.59) and L ² = 400, 2500, 10 000 or 40 000 in § 3.1; β = 0.59, L ² = 400, 2500, 10 000 or 40 000 and BR = 10 % to 75 % to study the effect of BR in § 3.2; and BR = 50 % or 75 %, L ² = 2500 or 10 000 and β = 0.1–0.9 to study the effect of β in § 3.3. Table 2 summarizes the parameter space of the present study. In a certain parameter range, the flow is unsteady. We use the data in five statistical cycles to calculate the time-averaged power spectral density (PSD) and u_rms quantities after the flow reaches the statistically steady state.

Table 2. Parameter space of the present study.

Before we perform the flow simulations, the effect of the length of downstream region l ₂ on the upstream recirculation is examined. Here we consider the parameter set (BR, β, L ², Wi) = (50 %, 0.59, 2500, 60). Three different l ₂ of 25H, 50H and 100H are investigated. The upstream recirculation lengths at different l ₂ are summarized in table 3, which indicates that l ₂ has a very small effect on the upstream recirculation. To save computational time, l ₂ = 25H is selected for the final simulation.

Table 3. Effect of l ₂ on the upstream recirculation length.

3.1. High-Weissenberg-number simulations for the FENE-P model at BR = 50 %

Simulation is first performed for the case with the parameter set (BR, β, L ²) = (50 %, 0.59, 2500). The set (BR, β) = (50 %, 0.59) is often used as a benchmark case for numerical simulation. Due to the numerical instability and heavy calculation burden, Wi was often set to less than 0.5 in previous studies (e.g. Fan, Tanner & Phan-Thien 1999; Alves, Pinho & Oliveira 2001; Hulsen et al. 2005). In our numerical simulation, relatively high-Wi flows (up to 100) are calculated by considering polymer molecular diffusion and using a small time step. Lee, Hwang & Cho (2021) indicated that the effect of molecular dissipation of polymer can be negligible only when Pe is large (Pe > ~10⁵) and Wi is small (Wi < ~1). Since the maximum Wi considered is 100 and Pe is equal to 40 in the present study, molecular dissipation of polymer should be taken into account.

The streamline distributions and the u velocity profiles are shown in figure 6. In a Newtonian fluid, the flow is almost symmetric with respect to the cylinder as shown in figure 6(a). As Wi increases, the velocity decreases gradually along the directions approaching the cylinder both upstream and downstream, and the downstream velocity decreases faster. Figure 6(f) indicates that the downstream u velocity almost linearly increases with x except for the region near the cylinder wall for different Wi, which is consistent with the experimental results of Haward et al. (2018). Note that the strain rate at the rear stagnation point of the cylinder is equal to zero (Haward et al. 2019). At Wi = 15, a recirculation zone starts to form upstream of the cylinder. Correspondingly, a region of negative u velocity can be observed in the upstream recirculation region as shown in figure 6(f). With a further increase in Wi, the upstream recirculation bubble becomes larger.

Figure 6. Streamlines for Newtonian (a) and viscoelastic (b–e) fluids at BR = 50 %, β = 0.59 and L ² = 2500 for different Wi and the corresponding u velocity profiles along y = 0 (f). The contours in (a–e) are coloured by $\sqrt {{u^2} + {v^2}} /\bar{U}$. All the streamlines are almost symmetric about the horizontal centreline of the cylinder.

To better understand the upstream recirculation, we follow Lee et al. (2007) and discuss the flow-type parameter defined as

(3.1)\begin{equation}\xi (x,y) = \frac{{|\dot{\boldsymbol{\gamma }}|- |\boldsymbol{\varOmega }|}}{{|\dot{\boldsymbol{\gamma }}|+ |\boldsymbol{\varOmega }|}},\end{equation}

where $|\dot{\boldsymbol{\gamma }}|= \sqrt {{\textstyle{1 \over 2}}\boldsymbol{D}:\boldsymbol{D}}$ and $|\boldsymbol{\varOmega }|= \sqrt {{\textstyle{1 \over 2}}\boldsymbol{\varOmega }:\boldsymbol{\varOmega }}$ are the magnitudes of the deformation rate tensor $\boldsymbol{D} = {\textstyle{1 \over 2}}(\boldsymbol{\nabla }\boldsymbol{u} + \boldsymbol{\nabla }{\boldsymbol{u}^T})$ and the vorticity tensor $\boldsymbol{\varOmega } = {\textstyle{1 \over 2}}(\boldsymbol{\nabla }\boldsymbol{u} - \boldsymbol{\nabla }{\boldsymbol{u}^T})$, respectively, which can be locally evaluated using the components of the velocity vectors. Parameter ξ is all-coordinate invariant (Lee et al. 2007). Here, ξ = −1 indicates solid body rotation, ξ = 0 simple shear and ξ = 1 pure extension. The flow strength in the extensional regions is quantified by the principal strain rate or the eigenvector of D (Astarita 1979):

(3.2)\begin{equation}{\lambda _1} = {\textstyle{1 \over 2}}\sqrt {{{({D_{11}} - {D_{22}})}^2} + 4D_{12}^2} .\end{equation}

In order to better describe the tensile strength in the extension area, a dimensionless parameter λ ₂ is defined by multiplying λ ₁ and ξ:

(3.3)\begin{equation}{\lambda _2} = \xi {\lambda _1}.\end{equation}

A local stretch Weissenberg number is defined as

(3.4)\begin{equation}Wi_{local}^{elastic} = \frac{\lambda }{{{u^2} + {v^2}}}\left[ {{u^2}\frac{{\partial u}}{{\partial x}} + uv\frac{{\partial v}}{{\partial x}} + uv\frac{{\partial u}}{{\partial y}} + {v^2}\frac{{\partial v}}{{\partial y}}} \right],\end{equation}

which describes the local effect on the Weissenberg number caused by stretching (positive $Wi_{local}^{elastic}$) or compression (negative $Wi_{local}^{elastic}$) of fluid parcels along the flow streamline direction.

The ξ and λ ₂ distributions are shown in figures 7(a) and 7(b), respectively. In a Newtonian fluid, both ξ and λ ₂ distribute symmetrically against the x and y axes with respect to the cylinder. The fluid parcel around the cylinder exhibits high extension. Just upstream of the cylinder, the flow velocity is reduced, resulting in compression. In the narrow gap between the cylinder surface and the channel wall with x < 0, flow velocity increases, resulting in stretching. Figure 7 shows that the upstream compression or stretching regions further expand upstream along the middle line (y = 0) and the gap middle line (y = 3/8H) as Wi increases. The λ ₂ distribution demonstrates that the strongest tensile strength region appears in the gap between the cylinder surface and the channel wall.

Figure 7. (a) Flow-type parameter distributions and (b) dimensionless λ ₂ distributions.

The local elastic stretch Weissenberg numbers $(Wi_{local}^{elastic})$ at y = 0 and y = 3/8H are extracted and shown in figure 8. As defined above, positive and negative $Wi_{local}^{elastic}$ denote fluid stretching and compression, respectively. At y = 0, the initial decrease in the flow velocity results in a negative $Wi_{local}^{elastic}$ for a long distance. Parameter $Wi_{local}^{elastic}$ initially decreases and then increases along the x direction. Generally, the location of the valley point moves upstream and the corresponding minimum $Wi_{local}^{elastic}$ decreases with increasing Wi. For Wi = 60, the minimum $Wi_{local}^{elastic}$ is about −19.8. At y = 3/8H, the increase in the flow velocity results in a positive $Wi_{local}^{elastic}$ for a long distance. Before approaching 2x/H = 0, Wi rapidly increases to the peak value and then experiences a sudden drop. The maximum $Wi_{local}^{elastic}$ is about 68.1 for Wi = 60, which occurs in the gap between the cylinder surface and the channel wall. The absolute $Wi_{local}^{elastic}$ in the gap (= 68.1) is much larger than that in front of the cylinder (= 19.8). The large $Wi_{local}^{elastic}$ in the gap reflects that the high elastic stress is concentrated there, which extends upstream near the wall region as shown in figure 9. A similar observation has been reported by Zhao et al. (2016) in wormlike fluids.

Figure 8. Variation of the local elastic stretch Weissenberg number along the x direction at (a) y = 0 and (b) y = 3/8H.

Figure 9. Contours of $1/(1 - {c_{kk}}/{L^2})$ for different Wi; $1/(1 - {c_{kk}}/{L^2})$ is positively correlated with the elastic stress.

A well-known dimensionless parameter which rationalizes these streamline instabilities is the M parameter introduced by McKinley, Pakdel & Öztekin (1996):

(3.5)\begin{equation}M = {\left[ {\frac{{\lambda \bar{U}}}{\Re }\frac{{{\tau_{11}}}}{{({\eta_0}|\dot{\boldsymbol{\gamma }}|)}}} \right]^{1/2}},\end{equation}

where $\lambda \bar{U} = l$ is the characteristic length over which perturbations to the base stress and velocity fields relax, $\Re$ is the streamline radius of curvature and ${\tau _{11}}$ is the streamwise tensile stress. We convert this definition of the criterion to a form amenable to local evaluation in flows by the substitution of characteristic values with local values. To consider local effective relaxation time $({\lambda _{eff}} = \lambda /f)$ and ignore solid-like rotational flows, Cruz et al. (2016) proposed a modified Pakdel–McKinley criterion:

(3.6)\begin{equation}{M^\ast } = {\left[ {\frac{{{\lambda_{eff}}|\boldsymbol{u}|}}{r}\frac{{{\tau_{11}}}}{{||\boldsymbol{\tau }|{|_F}}}} \right]^{1/2}},\end{equation}

where $r = |\boldsymbol{u}{|^3}/|\boldsymbol{u} \times \boldsymbol{\dot{u}}||$. Here we note that $\dot{\boldsymbol{u}}$ is the material derivative of the velocity vector, which is equivalent to $\boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{\nabla }\boldsymbol{u}$ for steady-state flow. The Frobenius norm $(||\cdot |{|_F})$ is used, so that the resulting ${\tau _{11}}/||\boldsymbol{\tau }||_F$ will vary between zero, when the normal stress is weak, and one, when the tensile normal stress dominates, in highly elastic shear or extensional flows. Therefore in strongly extensional flows, M* is approximately given by ${\lambda _{eff}}|\boldsymbol{u}|/r$.

Contours of M* for different Wi are shown in figure 10. The largest M* region occurs in the gap between the cylinder and the outer wall. In the upstream recirculation regions, we can also find non-negligible, though smaller, values of M*. From the definition of M*, high M* happens in highly elastic shear or extensional flows. Although the flow curvature in the recirculation zone is large, we only see a relatively small value of M* in the recirculation region because the extension in this region is not larger than that in the gap regions (figure 8). The distribution of M* suggests that the instability is most likely related to the high-shear region near the cylinder, which could be classified as the standard curve-streamline shear flow instability. This conclusion would be also consistent with that reported in Davoodi, Dominques & Poole (2019) where the effect of strain rate on the start point of purely elastic instability in elongational-dominated flows was investigated.

Figure 10. Contours of M* in the viscoelastic cylindrical wake flow with upstream instability for different Wi.

The above results indicate that $Wi_{local}^{elastic}$ and the stress dominate in the gap, compared with those in front of the cylinder. In this sense, the flow across the gap can be regarded as a main flow, while the upstream recirculation is a secondary flow. A primary–secondary flow model shown in figure 11 can be applied to characterize this upstream recirculation, i.e. a high-speed stretching gap flow and a relatively low-speed upstream recirculation.

Figure 11. The primary–secondary flow model for the present flow configuration.

The time evolutions of the drag coefficient for Wi = 70, 80, 90 and 100 are plotted in figure 12. To save computing time, the numerical result of Wi = 60 (at time of $2t\bar{U}/H = 525$) is set as the initial value for these higher-Wi flows. The drag acting on the cylinder increases with Wi when Wi > ~0.5 (refer to data in table 5 for low-Wi range and data here in high-Wi range), no matter whether the flow is steady (Mokhtari et al. 2022) or unsteady. The time-averaged drag coefficient $(\overline {{C_d}} )$ is ~532 at Wi = 80 but ~799 at Wi = 82.5. A sharp increase in $\overline {{C_d}} $ is observed when Wi increases from 80 to 82.5 as shown in figure 12(a), which results from the enhancement in the additional extensional viscosity due to flow fluctuations (Browne & Datta 2021) as shown in figure 12(b). A similar sudden increase of drag was also reported in the numerical simulation of Grilli et al. (2013). Therefore, the increase in $\overline {{C_d}} $ with increasing Wi is sharper in unsteady flow than in steady flow.

Figure 12. The solid black lines denote the temporal histories of the drag coefficient at (BR, β, L ²) = (50 %, 0.59, 2500) for (a) Wi = 70, (b) Wi = 80, (c) Wi = 90 and (d) Wi = 100. The dashed brown line in (d) denotes the temporal history of L_D.

When Wi < ~80, our calculation indicates that the flow is steady after developing for a long time. However, the flow becomes unsteady when Wi ≥ 80. At Wi = 80, the fluctuation of drag coefficient exhibits a fully developed periodic state with a single frequency demonstrated by the FFT analysis shown in figure 13(a). For Wi = 100 in figure 12(d), the variation of drag coefficient roughly maintains a periodic state that features more discrete frequencies demonstrated by the FFT analysis shown in figure 13(b). Interestingly, the quantitative analysis demonstrates St ₄ ≈ 4St ₁, St ₃ ≈ 3St ₁ and St ₂ ≈ 2St ₁, indicating the nonlinear effect of the period-doubling. Figure 12(d) also shows the corresponding time evolution of L_D for one period, which exhibits a trend similar to that of the C_d curve. The maximum drag of the cylinder corresponds to the longest upstream recirculation. The corresponding streamlines at four time instants within one period are shown in figure 14.

Figure 13. The FFT of the drag coefficients for (BR, β, L ²) = (50 %, 0.59, 2500) at (a) Wi = 80 and (b) Wi = 100.

Figure 14. Streamlines at time instants (a) t = t ₁, (b) t = t ₂, (c) t = t ₃ and (d) t = t ₄ marked in figure 12(d) for (BR, β, L ², Wi) = (50 %, 0.59, 2500, 100).

A sharp increase of drag between Wi = 80 and 82.5 under the increasing Wi process (i.e. the discontinuity on the $\overline {{C_d}} - Wi$ curve near Wi = 82.5) shown in figure 15(a) implies a subcritical transition. We use the result of Wi = 82.5 as the initial field to study this transition behaviour and gradually reduce Wi (i.e. the decreasing Wi process). It is found that the $\overline {{C_d}} - Wi$ curve maintains continuity until Wi = 67.5 and then experiences a sudden drop between Wi = 67.5 and 65 as shown in figure 15(a). Both figures 15 and 16 show that the unsteady flow becomes steady when Wi is decreased from Wi = 67.5 to Wi = 65. The hysteresis phenomenon shown in figure 15(b) implies that the transient transition is a subcritical Hopf bifurcation.

Figure 15. Variations of (a) the time-averaged drag coefficient $(\overline {{C_d}} )$ and (b) the root-mean-square drag coefficient (C_drms) with Wi for (BR, β, L ²) = (50 %, 0.59, 2500). Here IW and DW denote the increasing and decreasing Wi processes, respectively.

Figure 16. The temporal histories of the drag coefficient at (BR, β, L ²) = (50 %, 0.59, 2500) for (a) Wi = 67.5 and (b) Wi = 65 under the decreasing Wi process. The numerical result of Wi = 82.5 is selected as the initial field. Here DW denotes the decreasing Wi process.

For L ² = 10 000 and 40 000, we also observe the upstream recirculation. The non-dimensional recirculation length (L_D) as a function of Wi for L ² = 2500, 10 000 and 40 000 is summarized in figure 17(a). At high Wi, the flow becomes unsteady for L ² = 2500, as discussed above. For L ² = 40 000, the flow also becomes unsteady when Wi ≥ 70. However, for L ² = 10 000, the unsteady behaviour is only observed when 20 ≤ Wi ≤ 50. When the flow is unsteady and L_D varies with time, the maximum and minimum values of L_D at the corresponding Wi are provided, as shown in figure 17(a). Note that we only present the result up to Wi = 60 for L ² = 40 000 here. A further increase in Wi causes simulation divergence due to numerical instability, even if we use a very small time step. We observe a sudden increase of L_D before simulation divergence.

Figure 17. (a) The L_D–Wi relationship and (b) the L_D–Wi_N relationship at (BR, β) = (50 %, 0.59). When the flow becomes unsteady and L_D varies with time at high Wi, the time-averaged L_D is provided while the error bar indicates the corresponding maximum and minimum values. Here IW and DW denote the increasing and decreasing Wi processes, respectively.

Length L_D for L ² = 40 000 is shorter than those for L ² = 2500 or L ² = 10 000. However, Qin et al. (2019a) implies the upstream recirculation length is positively correlated with extensional viscosity. Variation of L_D with Wi_N is plotted in figure 17(b). As mentioned above, the definition of Wi_N is similar to that in Pan et al. (2013). The time-averaged L_D measured by Qin et al. (2019a) is also plotted in figure 17(b). Behaviour of L_D at L ² = 2500 is mostly consistent with the experimental results of Qin et al. (2019a). It is worth pointing out that in their experiment, the upstream recirculation length varies strongly with time when Wi_N is beyond 4.35. However, the fluctuation of L_D is not very obvious in our simulation when Wi_N > 13.5 for L ² = 10 000. Note that the present simulation is two-dimensional and the present rheological model and parameters may be different from those for the material used in the experiment. Besides, Re may be different between our simulation and their experiment. Thus, the temporal behaviour of L_D deviates from that reported by Qin et al. (2019a).

Figure 17 clearly shows that there is a critical Weissenberg number Wi_c for the onset of the upstream recirculation at each L ². Close to the transition, L_D increases smoothly from zero to a non-zero value with increasing Wi, and the relationship between L_D and Wi can be well described by a simple Landau-type quartic potential minimized as

(3.7)\begin{equation}Wi = W{i_c}(gL_D^2 + hL_D^{ - 1} + 1),\end{equation}

where g denotes the growth rate coefficient and h quantifies system imperfection that biases a transition to a favoured branch. Figure 18 shows the relationship between L_D and Wi for both the numerical results and the fitted data around the onset of the upstream recirculation region. The corresponding values of g and h are listed in table 4. The fitted Wi_c and g increase with L ² while h is always close to zero. When L ² = 10 000, the fitted curve only collapses well with the numerical results when Wi is less than 40. This implies that at BR = 50 %, variation in Wi may cause inherent change in flow state when L ² = 10 000. It may result from the velocity increase in the gap between the cylinder surface and the channel wall. The time-averaged streamlines and contours of the dimensionless u velocity for Wi = 40 and 50 are shown in figures 19(a) and 19(b). The streamlines indicate that the upstream recirculation elongates when Wi increases from 40 to 50. Meanwhile, the contours of the dimensionless u velocity show that the maximum velocity in the gap between the cylinder surface and the channel wall increases from 3.31 at Wi = 40 to 3.65 at Wi = 50. The maximum time-averaged dimensionless x-direction velocity in the flow field $\textrm{max(}\bar{u}/\bar{U})$ is extracted and plotted in figure 19(c). Obviously, a significant increase in velocity occurs between Wi = 40 and 50. At high Wi, the elastic stress is concentrated near the channel and cylinder walls, which acts to narrow the gap and accelerate the flow velocity there. For example, the maximum velocity in the gap is 3.65 for (Wi, L ²) =(50, 10 000); however, it is about 3 for a Newtonian fluid. The effect of BR on the upstream recirculation is discussed in the next section.

Figure 18. Variation of L_D with Wi at BR = 50 % for (a) L ² = 400, (b) L ² = 2500, (c) L ² = 10 000 and (d) L ² = 40 000. The solid line denotes the data fitted with Landau-type quartic potential while the symbols $\bullet$, ♦, ★ and $\ast$ denote the numerical results. The fitted Wi_c for each L ² is summarized in table 4. Here IW and DW denote the increasing and decreasing Wi processes, respectively.

Figure 19. Time-averaged flow streamlines and contours of the dimensionless u velocity at (BR, β, L ²) = (50 %, 0.59, 10 000) for (a) Wi = 40 and (b) Wi = 50. (c) Variation of the maximum time-averaged dimensionless u velocity $\textrm{max}(\bar{u}/\bar{U})$ with Wi.

Table 4. The fitting parameters for BR = 50 % in (3.7).

3.2. Effect of BR on viscoelastic upstream instability

Previous experiments indicated that the upstream recirculation occurs when BR is not less than 50 %. When BR is low, such as 10 %, only upstream streamline bending occurs (Haward et al. 2018). Thus, there is a critical value of BR which is around 50 % and determined by Wi and L ², below which the upstream recirculation will not occur. When BR is 50 % and Wi is low in the simulations discussed in § 3.1, the velocity in the midline (y = 0) changes from 1.5 at the inlet to 0 at the cylinder surface, which results in a compression $Wi_{local}^{elastic}$. However, the velocity along y = 3/8H changes from about 1.5 at the inlet to about 3.0 in the gap, which leads to a stretching $Wi_{local}^{elastic}$. The maximum compression $Wi_{local}^{elastic}$ is almost equal to the maximum stretching $Wi_{local}^{elastic}$ in this situation. When BR exceeds 50 %, the gradient of flow velocity in the gap is steeper than that upstream of the cylinder. That is, the stretching $Wi_{local}^{elastic}$ upstream is larger than the compression $Wi_{local}^{elastic}$ in the gap. The compression $Wi_{local}^{elastic}$ in the gap plays a leading role, which results in a primary–secondary flow as shown in figure 11. Moreover, with the increase of Wi, the flow velocity inside the gap may slightly increase (discussed in last paragraph in § 3.1). This indicates that the blockage effect is more severe and the apparent BR is increased (the increase of the flow rate in the centreline of the gap can be equivalent to the increase of BR). Therefore, the upstream recirculation also occurs when BR = 50 %.

It is reasonable to speculate that a certain boundary exists in the space of (BR, Wi) to distinguish the regimes with and without the upstream recirculation. Thus, we studied the effect of BR on Wi_c for (β, L ²) = (0.59, 2500), as shown in figure 20. Note that Wi does not exceed 100 in this test. Here Wi_c is not obtained by fitting (3.7). Instead, at a given BR, multiple simulations are performed near Wi_c. The bisection method is adopted to determine Wi_c by checking whether the upstream recirculation occurs. Our results imply that the upstream recirculation only occurs when BR is more than 16.7 %, i.e. the Wi _c–BR curve asymptotically approaches BR ~ 16.7 %. Hopkins et al. (2022a) found that for a viscoelastic wormlike micellar solution, the upstream recirculation is observed only when BR is greater than certain thresholds and Wi_c is almost inversely proportional to BR. Our numerical simulation results qualitatively agree with their experimental results. We also select five typical cases with (BR, Wi) = P ₁ (25 %, 20), P ₂ (25 %, 30), P ₃ (50 %, 20), P ₄ (62.5 %, 20) and P ₅ (75 %, 20) to demonstrate the instantaneous streamlines. Here P ₁ is under the Wi_c–BR curve and no upstream recirculation exists. In contrast, P ₂ to P ₅ are all located above the Wi_c–BR curve, and the upstream recirculation does occur. For a fixed Wi, as BR increases (P ₃, P ₄ and P ₅), the upstream recirculation becomes longer and behaves asymmetrically. This asymmetric behaviour is explained later in this section.

Figure 20. Variation of Wi_c with BR for (β, L ²) = (0.59, 2500). The instantaneous streamlines for five typical cases with (BR, Wi) = P ₁ (25 %, 20), P ₂ (25 %, 30), P ₃ (50 %, 20), P ₄ (62.5 %, 20) and P ₅ (75 %, 20) are plotted.

In the rest of this subsection we mainly consider high BR of 62.5 % and 75 % and discuss its effect on the upstream recirculation. Figures 21(a) and 21(b) plot L_D as a function of Wi for BR = 62.5 % and 75 %, respectively. When (BR, L ²) = (62.5 %, 2500), (62.5 %, 10 000), (75 %, 400) and (75 %, 2500), the growth of upstream recirculation with Wi close to the transition also satisfies (3.7), as shown in figures 22(a,b) and 23(a,b). However, when (BR, L ²) = (62.5 %, 40 000), (75 %, 10 000) and (75 %, 40 000), the L_D–Wi relationship for the increasing Wi process does not satisfy (3.7), as shown in figures 21(a,b) and 23(d). In fact, there is a sudden increase of L_D when Wi is larger than a certain value, indicating that the L_D–Wi relationship is discontinuous there. For example, when BR = 62.5 % and L ² = 10 000, L_D equals zero at Wi = 14.5. However, L_D suddenly changes to 2.186–2.5697 when Wi > 14.5. Thus, the bifurcation for these parameter sets is subcritical. For decreasing Wi, the results for (BR, L ²) = (75 %, 10 000) and (75 %, 40 000) still satisfy (3.7), as shown in figure 23(b,c). The flow instability is affected by external disturbance. In previous experiments (Kenney et al. 2013; Shi et al. 2015; Zhao et al. 2016; Qin et al. 2019a; Haward et al. 2021), no sudden increase in the L_D–Wi relationship was observed, which may be caused by strong disturbance in those experiments. The appearance of upstream recirculation is a more stable form. These scenarios are consistent with our results obtained in the decreasing Wi process, where strong disturbance is also introduced by the initial fields. It is worth pointing out that when (BR, L ²) = (62.5 %, 40 000), L_D jumps suddenly with increasing Wi around Wi = 15–16, as shown in figure 21(a). This implies that hysteresis may exist. However, we are not able to carry out further investigation due to serious numerical instability.

Figure 21. The L_D–Wi relationship for (a) BR = 62.5 % and (b) BR = 75 % (for the decreasing Wi process, the numerical results of Wi = 10 are selected as the initial fields for L ² = 10 000, while the numerical results of Wi = 9 are selected as the initial fields for L ² = 40 000). Variation in the root-mean-square L_D with Wi at L ² = 40 000 for (c) BR = 62.5 % and (d) BR = 75 %. Here IW and DW denote the increasing and decreasing Wi processes, respectively.

Figure 22. The growth of upstream recirculation with Wi at BR = 62.5 % for (a) L ² = 2500 and (b) L ² =10 000. The solid line denotes the data fitted with Landau-type quartic potential while the symbols ♦ and ★ denote the numerical results.

Figure 23. The growth of upstream recirculation with Wi at BR = 75 % for (a) L ² = 400, (b) L ² = 2500, (c) L ² = 10 000 and (d) L ² = 40 000. The solid line denotes the data fitted with Landau-type quartic potential while the symbols $\bullet$, ♦, ■ and $\ast$ denote the numerical results. Here DW denotes the decreasing Wi process.

Figures 21(c) and 21(d) plot the root-mean-square L_D as a function of Wi at L ² =40 000 for BR = 62.5 % and 75 %, respectively. Comparison between figures 21(a) and 21(c) indicates that the critical Wi for the onset of the upstream recirculation and the non-zero fluctuation in the upstream recirculation region almost coincide at Wi ~ 15. The variations of the root-mean-square L_D with Wi for both the increasing and decreasing Wi processes in figure 21(d) also indicate that bifurcation is subcritical.

Parameter Wi_c for different BR can be obtained from figures 18 and 21. With an increase in BR, Wi_c decreases. With an increase in L ², Wi_c becomes larger (the increasing Wi process) for most cases. Yamani & McKinley (2023) defined an important parameter Wi/L for flow instability. This parameter implies that L ² is positively related to Wi_c, which is consistent with the present results. Therefore, the onset of the upstream recirculation is promoted and Wi_c is smaller when L ² is smaller. Moreover, Wi_c is influenced by how the flow field is initialized, i.e. the increasing or decreasing Wi process mentioned in § 2.3. The value of Wi_c for the decreasing Wi process may be much lower than that for the increasing Wi process. For example, Wi_c at (BR, L ²) = (75 %, 10 000) is about 2.75 for the decreasing Wi process and about 6.75 for the increasing Wi process.

At high Wi, such as (BR, Wi) = (62.5 %, 20) or (75 %, 10), L_D only slightly changes with the variation in Wi. Length L_D may reach the saturation value in this situation. For BR = 62.5 % or 75 %, the saturation value is large when L ² is large, which is different from the trend for BR = 50 %. A large L ², i.e. a longer molecular chain length, means higher potential of extensional viscosity. Qin et al. (2019a) imply the upstream recirculation length is positively correlated with extensional viscosity, which is qualitatively consistent with our results at BR = 62.5 % or 75 %.

When BR = 62.5 % or 75 %, the upstream recirculation exhibits strong time-dependence behaviour at high Wi when L ² = 40 000. The u velocity contours and streamlines at time instant $2t\bar{U}/H = 1710$ for (BR, L ², Wi) = (75 %, 40 000, 9) are shown in figure 24(a). Under this parameter set, the flow is unsteady and the length and shape of the upstream recirculation vary with time. However, the fluctuation amplitude in the upstream recirculation region is not the maximum in the flow field. Instead, the strongest temporal fluctuation occurs in the gap region, as shown by the u_rms distribution in figure 24(b). Downstream of the cylinder, the high-velocity region is concentrated on both sides of the centreline, as shown in figure 24(a). Correspondingly, the flow fluctuation also is high there.

Figure 24. (a) The u velocity contours and streamlines at time instant $2t\bar{U}/H = 1710$ and (b) the corresponding u_rms distribution. The parameter set used in this simulation is (BR, L ², Wi) = (75 %, 40 000, 9).

We carefully examine the symmetry of the flow field and find that the symmetry against the horizontal centreline is obviously broken for BR = 62.5 % and 75 % when L ² is large, such as L ² = 40 000. Figure 25 plots the variations of the averaged velocity (space average) in the upper and lower gaps with time at (BR, L ², Wi) = (75 %, 40 000, 9). Both averaged velocities show strongly oscillating behaviour. The corresponding oscillating amplitudes for upper and lower averaged velocities are opposite, in order to satisfy the continuity equation. The upper averaged velocity fluctuates around a mean value of 4.168, which is larger than the 3.898 for the lower averaged velocity.

Figure 25. Time histories of the averaged velocity (space average) in the upper and lower gaps for (BR, L ², Wi) = (75 %, 40 000, 9).

A parameter, following Hopkins et al. (2021), is defined to evaluate the flow asymmetry:

(3.8)\begin{equation}I = \frac{{|\overline {{Q_u}} - \overline {{Q_l}} |}}{{\overline {{Q_u}} + \overline {{Q_l}} }},\end{equation}

where $\overline {{Q_u}} $ and $\overline {{Q_l}} $ are the volumetric flow rate through the upper and lower gaps, respectively. Variation of I with Wi for (BR, L ²) = (75 %, 40 000) is presented in figure 26. The flow symmetry is only maintained at Wi = 0, i.e. I = 0 for Newtonian flow. When Wi > 0, the viscoelastic flow exhibits complex asymmetric behaviour for high BR and high L ². For the increasing Wi process, I initially increases with Wi until Wi = 5.5. Then I decreases with further increase in Wi until Wi = 6.75. Note that for (BR, L ²) = (75 %, 40 000), the onset of the upstream recirculation occurs at Wi_c ~ 6.75. When Wi > 6.75, I increases again with Wi. For the decreasing Wi process, I decreases with decreasing Wi from Wi = 9 to 5.5. However, I increases with a further decrease in Wi until Wi = 4.75, which corresponds to Wi_c below which the upstream recirculation disappears. When Wi < 4.75, I decreases again with decreasing Wi. The increasing and decreasing Wi paths do not coincide for 4.75 < Wi < 6.75, which indicates a typical hysteresis phenomenon.

Figure 26. Variation of flow asymmetry parameter I with Wi for (BR, L ²) = (75 %, 40 000). Here IW and DW denote the increasing and decreasing Wi processes, respectively.

3.3. Effect of β and L ² on viscoelastic upstream instability

Figure 18(b) shows that the L_D–Wi relationship for (BR, L ²) = (50 %, 2500) approximately satisfies a Landau-type quartic potential near Wi_c. However, figure 21(b) shows that the variation of L_D with Wi for (BR, L ²) = (75 %, 10 000) under the increasing Wi process does not approximately fit a Landau-type quartic potential near Wi_c, but shows the hysteresis phenomenon. Thus, these two parameter sets are specially selected to investigate the effect of β. For (BR, L ²) = (50 %, 2500), β ranges from 0.1 to 0.9 and Wi is no more than 30. For (BR, L ²) = (75 %, 10 000), β ranges from 0.15 to 0.9 and Wi is no more than 20. Within these parameter spaces, the flow fluctuations are small and L_D only varies weakly with time. Thus, the effect of β on time-dependent stability is not discussed here.

In experiments, β decreases with increasing polymer concentration. The FENE-P model can describe the shear-thinning behaviour of a viscoelastic fluid, usually for the case when L ² is small. The smaller the value of β, the more obvious the shear-thinning effect in the flow. In our study, L ² is large and the influence of β on shear thinning can be ignored (Tamano et al. 2020). However, a variation in β may affect the distribution of elastic stress as discussed below.

We first consider (BR, L ²) = (50 %, 2500). The variations of L_D with Wi for different β are shown in figure 27. Near the upstream recirculation transition point, the L_D–Wi curves all approximately satisfies the Landau-type quartic potential. A smaller β corresponds to a smaller Wi_c. For example, Wi_c is between 23 and 24 at β = 0.9, while Wi_c is between 8 and 9 at β = 0.2. For the same Wi, a lower β corresponds to a larger L_D. For example, L_D is 0.3242 at β = 0.9 while L_D is 2.1010 at β = 0.2, for a fixed Wi = 30.

Figure 27. The L_D–Wi relationship at (BR, L ²) = (50 %, 2500) for different β.

For (BR, L ²) = (75 %, 10 000), L_D for different (Wi, β) is evaluated and shown in figure 28. Parameter β has a nonlinear effect on L_D. For both the increasing and decreasing Wi processes, before the onset of upstream recirculation, the L_D–Wi relationships are identical and satisfy (3.7) at β = 0.9 or 0.75. However, the L_D–Wi relationships for the increasing and decreasing Wi processes at β = 0.59, 0.45, 0.3, or 0.15 are different. For the increasing Wi process, L_D experiences a sudden jump with increasing Wi, and the L_D–Wi relationship shows a discontinuity. Interestingly, the discontinuity point (Wi_c) is located between 6.5 and 7.0 for all β. These features indicate that the transition is subcritical bifurcation at β = 0.59, 0.45, 0.3 or 0.15. Thus, we can conclude that β affects the bifurcation type, which due to the amplifying effect on elastic stress with decreasing β. Note that β also shows a similar effect in other viscoelastic flows. For example, in elasto-inertial pipe flow, recent experimental and theoretical studies indicated supercritical bifurcation for large β and subcritical bifurcation for small β (Samanta et al. 2013, Chandra, Shankar & Das 2020, Choueiri et al. 2021, Wan, Sun & Zhang 2021). For the decreasing Wi process, all the L_D–Wi relationships are still uninterrupted to satisfy (3.7).

Figure 28. (a) The L_D–Wi relationship and (b) the L_D–β relationship at (BR, L ²) = (75 %, 10 000). In (b), the larger L_D is selected if L_D obtained in the increasing and decreasing Wi processes (denoted by IW and DW, respectively) at the same parameters is not identical.

The effect of L ² on the upstream flow behaviour has been discussed at various points in the above subsections. We thus provide a brief summary in this subsection. Parameter L ² has a noticeable effect on Wi_c for the onset of the upstream recirculation. Generally, a larger L ² corresponds to a higher Wi_c as shown in figures 17(a) and 21(a,b). Also, L ² affects subcritical behaviour. For large BR with high L ² (such as L ² = 40 000 for BR = 62.5 %, and L ² = 10 000 or 40 000 for BR = 75 %), the upstream instability is subcritical. The effect of L ² on subcritical behaviour is similar to that of 1 − β discussed in § 3.3. For larger BR (e.g. 62.5 % and 75 %), the unsteady flow is more likely to appear at a larger L ².

4. Concluding remarks

In this paper, we confirm that the existing macroscopic viscoelastic constitutive relationship (FENE-P model) is still qualitatively applicable for predicting the viscoelastic upstream instability through numerical simulation. By applying the square root reconstruction method and considering the molecular dissipation effect to stabilize the numerical simulations, we studied the viscoelastic flow over a circular cylinder in a narrow channel at a very low Re and over a wide range of Wi.

Upstream recirculation is observed when Wi is beyond a certain critical value Wi_c and BR is larger than around 16.7 %. The occurrence of upstream recirculation is related to high local stretch Weissenberg number $Wi_{local}^{elastic}$ and the high elastic stress in the narrow gap between the cylinder surface and the channel wall, which is consistent with that reported by Zhao et al. (2016). This flow instability around the narrow gap with a high local stretch can also be interpreted by the curve-streamline shear flow instability (Davoodi et al. 2019). For BR = 50 %, our simulation results on the upstream recirculation are basically consistent with the experimental results of Qin et al. (2019a) for the time-averaged behaviour.

Higher BR could precipitate the onset of the upstream recirculation, that is, Wi_c becomes lower. The viscosity ratio β has less effect on Wi_c (the discontinuous point, the increasing Wi process), but has a nonlinear influence on L_D at high BR (75 %) and relatively high L ² (10 000). For BR = 50 % and L ² = 2500, a lower β corresponds to a lower Wi_c. Close to the onset of upstream recirculation, L_D with Wi satisfy the Landau-type quartic potential when BR, 1 − β and L ² are not very large, e.g. (BR, β, L ²) = (50 %, 0.59, 2500). Close to the unsteady transition, the occurrence of the hysteresis phenomenon (signalling a subcritical bifurcation) depends on L ² and β (Burshtein et al. 2017). For the cases with large L ², small β and large BR, the flow often exhibits subcritical behaviour. In this situation, strong disturbance should be added in the initial flow field so that the corresponding simulation results can be consistent with the experimental results.

When Wi exceeds a critical Wi ~ 0.5 at BR = 50 % (the critical Wi depends on BR), the drag on the cylinder increases with Wi, although the flow remains steady. When Wi continues to increase beyond a certain threshold, the upstream recirculation becomes unsteady. The drag on the cylinder increases more sharply in this range.

Our results advance the understanding of the underlying mechanism of this novel flow instability, which may also be relevant to a number of interesting phenomena observed in viscoelastic flows, e.g. stationary dead zones associated with complex geometries such as multiple tandem circular cylinders (Shi & Christopher 2016) and porous media (Kawale et al. 2017).

A cylinder in a channel is often used as a prototype configuration to experimentally study the elastic instability. The upstream instability may extend downstream when Wi is larger. It is interesting to study the elastic wave in the channel through numerical simulations. However, only the elastic instability in a channel without a cylinder was studied (Morozov 2022). Our simulation could provide a platform for flows in a channel with cylinder(s).

The upstream recirculation phenomena are governed by many parameters and the present study does not explore all the parameter space of (BR, Re, Wi, L ², β). These parameters need to be investigated in more detail. In experiments of Shi et al. (2015), upstream instability eventually occurs when the elastic Mach number $Ma = \sqrt {Re \times Wi}$ is beyond ~10 at Re ~ 1, which is affected by both channel geometry and fluid properties. We believe that Re is an important dimensionless number for such instability. However, the effect of Re is not discussed in this paper and will be studied in our future work.