1. Introduction
Mixing is ubiquitous in both natural and industrial environments. Applications span a vast range of Reynolds numbers and length scales (see figure 1 in Ottino Reference Ottino1990 for illustration). From the coffee we drink, household cleaning products and oil extraction to the human digestive system and pharmaceutical production, various materials undergo mixing processes daily.
Despite its prevalence, mixing remains one of the more challenging paradigms in engineering to systematically define, frame and understand (Spencer & Wiley Reference Spencer and Wiley1951; Ottino Reference Ottino1990; Villermaux Reference Villermaux2019). In the simplest context, mixing entails the homogenisation of a passive tracer (level-1); however, it can be more intricately tied to the flow dynamics, as when rheology depends on tracer concentration (level-2), or chemical reactions occur during the process (level-3) (Dimotakis Reference Dimotakis2005).
A significant body of literature focuses on level-1 mixing, i.e. mixing of a passive dye in fluid. Even within this limited scope, the parameter space is extensive: mixing may be in line, active or passive, or take place in a stirred tank. Factors like domain geometry, impeller shape, size, position, stirring protocol and speed all influence mixing behaviour. Additionally, fluid rheology plays a critical role. Given the vast range of parameters and the complexity of the problem, most studies have focused on Newtonian fluids, and mixing remains an active research area (see Warhaft Reference Warhaft2000; Peltier & Caulfield Reference Peltier and Caulfield2003; Wunsch & Ferrari Reference Wunsch and Ferrari2004; Caulfield Reference Caulfield2021 and references therein).
Many fluids in polymer processing, food engineering, bioengineering, physiology and chemical engineering are non-Newtonian, with a subset exhibiting yield stress, such as polymeric gels, muds, paints and cosmetics. Yield-stress fluids are highly viscous materials that flow only when the applied shear stress exceeds a threshold known as the yield stress (Balmforth, Frigaard & Ovarlez Reference Balmforth, Frigaard and Ovarlez2014; Coussot Reference Coussot2014; Bonn et al. Reference Bonn, Denn, Berthier, Divoux and Manneville2017). It was recognised early on that turbulent mixing in these fluids is economically and technologically impractical. Spencer & Wiley (Reference Spencer and Wiley1951) proposed streamline mixing, achieved by continuously deforming the fluid to (a) increase the surface area of the interface and (b) distribute it throughout the material volume. With few exceptions (Derksen Reference Derksen2013; Daneshvar Garmroodi & Karimfazli Reference Daneshvar Garmroodi and Karimfazli2024), studies of mixing in yield-stress fluids focus on the mixing of passive dyes.
Initial efforts to understand mixing of non-Newtonian fluids were dedicated to establishing a relationship between impeller speed (in stirred tanks) and the fluid shear rate. According to Metzner & Otto (Reference Metzner and Otto1957), mixing flows of non-Newtonian fluids were qualitatively understood at best in the 1950s, a sentiment that remained accurate for decades. One of the first studies on mixing yield-stress fluids was conducted by Solomon et al. (Reference Solomon, Elson, Nienow and Pace1981), who experimentally identified the well-mixed regions in Xanthan gum and Carbopol solutions in stirred tanks. The coexistence of flowing and stagnant regions presented a challenge, as the latter remained unmixed. Whitcomb & Macosko (Reference Whitcomb and Macosko1978) introduced the term cavern to describe the well-mixed region where the fluid was yielded.
Although debate remains regarding the existence of a true yield stress and optimal measurement methods (Barnes Reference Barnes1999; Divoux, Barentin & Manneville Reference Divoux, Barentin and Manneville2011; Dinkgreve et al. Reference Dinkgreve, Paredes, Denn and Bonn2016), viscoplastic models are widely used to analyse and predict flows of yield-stress fluids (Mitsoulis & Tsamopoulos Reference Mitsoulis and Tsamopoulos2017). These models consider the fluid rigid below the yield stress and a flow with shear rate-dependent viscosity when the yield stress is exceeded.
Solomon et al. (Reference Solomon, Elson, Nienow and Pace1981) used viscoplastic models and a yield criterion to estimate cavern size by assuming it to be spherical. Subsequent studies explored various impellers and developed cavern size estimates based on different simplified geometries; see, for example, Galindo et al. (Reference Galindo, Argüello, Velasco, Albiter and Martínez1996), Tanguy et al. (Reference Tanguy, Bertrand, Labrie and Brito-De La Fuente1996), Amanullah, Hjorth & Nienow (Reference Amanullah, Hjorth and Nienow1998), Pakzad et al. (Reference Pakzad, Ein-Mozaffari, Upreti and Lohi2013), Sossa-Echeverria & Taghipour (Reference Sossa-Echeverria and Taghipour2015) and Ameur (Reference Ameur2020). Generally, a yield stress and shear-thinning viscosity reduce the mixing rate and extent and the cavern size. However, a mechanistic understanding of how the flow dynamics and mixing are interlinked, especially with respect to rheological parameters, remains elusive.
Seminal works by Aref (Reference Aref1984) and Ottino (Reference Ottino1989) demonstrated the role of chaotic flows in effective mixing, emphasising the necessity of three-dimensional or transient two-dimensional flows. A parallel branch of research has since modelled mixing by analysing dynamical systems, with a primary focus on chaotic mixing in two-dimensional time-periodic flows; for an overview of studies of Newtonian fluids, see Aref et al. (Reference Aref2017).
Niederkorn & Ottino (Reference Niederkorn and Ottino1994) numerically investigated chaotic advection in journal bearing flows for shear-thinning fluids, examining tracer advection, periodic points, stretching along unstable manifolds and stretching rate of fluid elements. They found that shear-thinning viscosity decreases the amount of stretching. Fan, Phan-Thien & Tanner (Reference Fan, Phan-Thien and Tanner2001) studied advective mixing of viscoplastic fluids between eccentric cylinders, showing that tracer coverage may depend on initial tracer position, and chaotic advection can be achieved through alternating cylinder rotations. Their results highlight qualitative transitions but do not predict when such transitions will occur.
Experimental studies by Wendell et al. (Reference Wendell, Pigeonneau, Gouillart and Jop2013) explored mixing in a rotating tank filled with a yield-stress fluid stirred by cylindrical rods rotating with constant angular velocity in an eggbeater configuration. Higher period ratios and lower yield stresses enhanced mixing, although mixing efficiency decreased in resonance conditions. The intermittent yielding and unyielding near the tank wall was hypothesised to explain the decreased mixing of yield-stress fluids. Further experiments by Boujlel et al. (Reference Boujlel, Pigeonneau, Gouillart and Jop2016) in the same set-up measured mixing rate using dye concentration variance, showing that mixing consists of rapid stretching and folding followed by slower diffusion-dominated mixing. They concluded that mixing rate is proportional to the volume of highly sheared fluid during each rod rotation.
In summary, while the qualitative impact of yield stress on mixing is understood – yield stress limits cavern size and filament stretching, thus decreasing mixing rate – a mechanistic description connecting the fluid dynamics to transitions in yield-stress fluid mixing remains absent. Decades after Niederkorn & Ottino (Reference Niederkorn and Ottino1994), design procedures still rely heavily on empiricism with a limited fundamental understanding of the fluid mechanics (see Paul et al. Reference Paul, Atiemo-Obeng and Kresta2004; Uhl Reference Uhl2012).
The primary objective of this manuscript is to identify and elucidate the mechanisms behind different mixing regimes and localisation in yield-stress fluids within a periodically stirred domain. We consider the simplest and most common viscoplastic model, the Bingham model, thus neglecting the thixotropy and elasticity typically associated with real yield-stress fluids. Potential influences of dye concentration on fluid properties and density are similarly neglected. We also adopt a minimalistic stirrer geometry and stirring strategy to avoid the added complexity of geometric effects and to isolate the roles of yield stress and localisation. We consider an infinite two-dimensional domain filled with a quiescent viscoplastic fluid stirred by a cylinder moving at constant speed along a circular path. By exploring a range of mixing speeds and yield stresses, we aim to characterise the flow and the mixing dynamics and establish a mechanistic link between them.
This problem is closely related to the classical flow past a circular cylinder. Over the past few decades, numerous studies have examined the flow of yield-stress fluids around a cylinder, either being drawn through the fluid or moving due to buoyancy (referred to as the resistance and mobility problems, respectively Putz & Frigaard Reference Putz and Frigaard2010). These studies typically focus on regimes where inertial effects are negligible and vortex shedding does not occur. Numerical, experimental and analytical results have been developed for predicting the terminal velocity, drag forces and the critical yield stress beyond which motion is arrested (see, e.g. Tokpavi et al. Reference Tokpavi, Magnin and Jay2008, Reference Tokpavi, Jay, Magnin and Jossic2009; Wachs & Frigaard Reference Wachs and Frigaard2016; Chaparian & Frigaard Reference Chaparian and Frigaard2017a , Reference Chaparian and Frigaardb ). However, parameter regimes characterised by weak inertia or dominated by large unyielded regions are not conducive to strong advective mixing.
The parameter space considered in this study is therefore complementary to the existing body of literature. We explore a broad range of Reynolds numbers within the laminar regime and yield-stress values away from the critical limit, where unyielded zones dominate and regularised numerical models are less accurate. This allows us to capture the mixing dynamics that involves vortex shedding and yield-stress-driven localisation.
The remainder of this paper is organised as follows: § 2 presents the model problem, governing equations and numerical methods. Section 3 discusses the flow dynamics and mixing regimes using representative cases and maps out mixing mechanisms and flow regimes at different mixing speeds and yield stresses. Finally, § 4 summarises our findings.
2. Problem set-up
2.1. Model problem
We investigate the stirring of a viscoplastic fluid (VPF) using a circular stirrer of diameter
$\hat {d}_s$
, which moves at a constant speed (
$\hat {r}_o \hat {\varOmega }$
) along a circular trajectory with radius
$\hat {r}_o = c \hat {d}_s$
(solid white line in figure 1). Here,
$c$
is the dimensionless stirring radius, fixed at
$c = 2$
in this study, and
$\hat {\varOmega }\gt 0$
denotes the angular velocity of the stirrer. Dimensional quantities are denoted by a
$\hat {.}$
symbol, while dimensionless quantities are written without it. To simulate mixing in an infinite domain, we employ a circular computational domain with a radius
$\hat {R}\gg \hat {r}_o$
(
$\hat {R} =34 \hat {r}_o$
), ensuring minimal boundary effects due to the domain’s sufficiently large size. To monitor mixing, the fluid in the bottom half of the domain is marked with a passive dye (shown in red in figure 1), where
$\alpha = 1$
, while the rest of the fluid is dye free (shown in blue), with
$\alpha = 0$
. The no-slip boundary condition is imposed on the walls of both the vessel and the stirrer.
Figure 1. Schematic of the domain geometry and initial conditions. The solid white line indicates the stirrer’s path. The red and blue colours indicate the dyed and dye-free regions. Note that the figure is not to scale.
The fluid is modelled using the Bingham model. The dimensionless form of the constitutive equation is given by
\begin{align} \left \{ \begin{array}{lr} \displaystyle \boldsymbol{\tau } = \left ( \frac {\textit{Bn}}{\dot {\gamma }} + 1 \right ) \dot {\boldsymbol{\gamma }} & \text{if} \quad \tau \geqslant Bn, \\[10pt] \dot {\boldsymbol{\gamma }} = 0 & \text{if} \quad \tau \lt Bn ,\end{array} \right . \end{align}
where
$\boldsymbol{\tau }$
and
$\boldsymbol{\dot {\gamma }}$
are the deviatoric stress and rate of strain tensors, respectively, and
\begin{align} \tau = \sqrt {\displaystyle \frac {1}{2}\ \boldsymbol{\tau }:\boldsymbol{\tau }} \quad \text{and} \quad \dot {\gamma } = \sqrt {\displaystyle \frac {1}{2}\ \boldsymbol{\dot {\gamma }}:\boldsymbol{\dot {\gamma }}} \end{align}
represent second invariants of these tensors. The Bingham number,
${\textit{Bn}}$
, is defined by
\begin{align} Bn = \frac {\hat {\tau }_{y}}{\hat {\mu }\hat {\varOmega }} ,\end{align}
where
$\tau _{y}$
and
$\hat {\mu }$
are the fluid’s yield stress and plastic viscosity.
The flow is governed by the Cauchy’s equations of motion and continuity while the dye concentration is described by an advection–diffusion equation
\begin{equation} \begin{aligned} \frac {\partial \boldsymbol{u}}{\partial {t}} + \boldsymbol{u} \boldsymbol{\cdot }\boldsymbol{\nabla } \boldsymbol{u} + \boldsymbol{\nabla } P &= \frac {1}{{\textit{Re}}} \boldsymbol{\nabla } \boldsymbol{\cdot }\boldsymbol{\tau } \\ \boldsymbol{\nabla } \boldsymbol{\cdot }\boldsymbol{u} &= 0 \\ \frac {\partial \alpha }{\partial t} + \boldsymbol{\nabla } \boldsymbol{\cdot }(\boldsymbol{u} \alpha ) &= \frac {1}{Pe} {\nabla} ^2 \alpha. \end{aligned} \end{equation}
Here,
$\boldsymbol{u}$
and
$P$
are the dimensionless velocity and pressure, respectively. The characteristic scales for length, time, velocity, shear stress and pressure are
$\hat {r}_o$
,
$\hat {\varOmega }^{-1}$
,
$\hat {r}_o \hat {\varOmega }$
,
$\hat {\mu } \hat {\varOmega }$
and
$\hat {\rho } \hat {r}_o^2 \hat {\varOmega }^2$
, where
$\hat{\rho}$
is the density of the fluid. Results, however, are presented in terms of the stirrer’s period,
$\hat {T}_{\textit{stirrer}}= 2\pi /\hat {\varOmega }$
\begin{align} T = \frac {\hat {t}}{\hat {T}_{\textit{stirrer}}} = \frac {t}{2\pi } .\end{align}
The Reynolds (
${\textit{Re}}$
) and Péclet numbers (
${\textit{Pe}}$
) are defined as
\begin{align} {\textit{Re}} &= \frac {\hat {\rho } \hat {\varOmega } \hat {r}^2_o}{\hat {\mu }} ,\end{align}
\begin{align} \textit{Pe} &= \frac {\hat {\varOmega } \hat {r}^2_o}{\hat {D}_m} ,\end{align}
where
$\hat {D}_m$
is the diffusion coefficient of dye. In this study, Péclet number is held constant at
${\textit{Pe}} = 10^3$
.
Definitions and ranges of the relevant dimensionless groups are summarised in table 1.
Table 1. Dimensionless groups governing the model problem.
To characterise the rate of mixing, a normalised variance of the dye,
$\sigma ^{2}_{R_{sd}}$
, is defined over a circular subdomain
$A_{R_{sd}}$
of radius
$R_{sd}$
that is concentric with the stirrer’s path
\begin{equation} \sigma ^{2}_{R_{sd}} = \frac {1}{A_{R_{sd}}} \int _{A_{R_{sd}}} \left (1 - \frac {\alpha }{\bar {\alpha }}\right )^{2} {\rm d}A .\end{equation}
Here,
$\bar {\alpha }$
is the average dye concentration over the subdomain, for all
$R_{sd}$
. For given values of the governing dimensionless parameters,
$R_{sd}$
is chosen to ensure the subdomain captures the region where dye concentration is affected by the stirring (within the timeframe of interest).
To quantify the kinetic energy,
$\textit{KE}$
is defined as
\begin{equation} {\textit{KE}} = \sqrt { \int _{A} |\boldsymbol{u}|^{2} {\rm d}A} ,\end{equation}
where
$\textit{KE}$
is the kinetic energy,
$|\boldsymbol{u}|$
is the speed and
$A$
represents the flow domain.
2.2. Numerical method
OpenFOAM – Open Source Computational Fluid Dynamics (2004) version 6 was used for numerical simulations. The twoLiquidMixingFoam solver, which employs the PIMPLE (PISO–SIMPLE) algorithm, which combines PISO (Pressure-Implicit with Splitting of Operators) and SIMPLE (Semi-Implicit Method for Pressure-Linked Equations) algorithms, to decouple pressure and velocity in the governing equations, was utilised. For temporal discretisation, we applied the second-order Crank–Nicolson scheme. The viscous terms are discretised using a second-order centred scheme, while the convective (inertial) terms are treated with a second-order upwind-biased scheme. The second-order vanLeer scheme is used for the concentration field transport. Adaptive time stepping was implemented using a constant Courant–Friedrichs–Lewy number (set to
$0.05$
to ensure sufficient accuracy in transient results).
In this study, we adopt a modified version of the bi-viscosity model originally developed by Tanner & Milthorpe (Reference Tanner and Milthorpe1983)
\begin{equation} \boldsymbol{\tau }(\dot {\gamma }) = \begin{cases} \displaystyle \left(1 + \dfrac {\textit{Bn}}{\dot {\gamma }_{\textit{cr}}}\right) \boldsymbol{\dot {\gamma }} & \text{if} \quad \dot {\gamma } \leqslant \dot {\gamma }_{\textit{cr}} ,\\[10pt] \displaystyle \left(1 + \dfrac {\textit{Bn}}{\dot {\gamma }}\right) \boldsymbol{\dot {\gamma }} & \text{if} \quad \dot {\gamma } \gt \dot {\gamma }_{\textit{cr}} .\end{cases} \end{equation}
Here,
$\dot {\gamma }_{\textit{cr}}$
is the critical strain rate. Because the numerical scheme involves regularisation, the accuracy of predictions for the unyielded regions decreases significantly at high
${\textit{Bn}}$
(Frigaard & Nouar Reference Frigaard and Nouar2005; Ahmadi, Olleik & Karimfazli Reference Ahmadi, Olleik and Karimfazli2021), making the scheme unsuitable for exploring the limit behaviour in that regime. For the parameter ranges considered here,
$\dot {\gamma }_{\textit{cr}}$
is chosen so that the error in the average radius of the yielded regions remains safely below 1 %; specifically,
$\dot {\gamma }_{\textit{cr}} = \hat {\dot {\gamma }}_{\textit{cr}} / \hat {\varOmega } = 10^{-4}$
.
To verify grid independence, five different mesh sizes were tested. Figure 2 shows the evolution of normalised variance and velocity norm, along with the corresponding relative errors for different mesh sizes. The finest mesh (
$1.3\times 10^5$
) was used to estimate the relative error. For the remainder of the simulations, we used a mesh size of
$8\times 10^4$
, which resulted in a relative error of less than 1 % for both the normalised variance and velocity norm. Further details on the benchmarking and validation of the numerical solver can be found in Daneshvar Garmroodi & Karimfazli (Reference Daneshvar Garmroodi and Karimfazli2024).
Figure 2. Time evolution of (a) normalised variance of dye concentration for different mesh sizes, (b) relative error in dye concentration variance,
$\displaystyle E_R(\sigma ^2_{6}) = ({\sigma ^2_{6} - \sigma ^2_{6,\ 1.3\times 10^5}})/({\sigma ^2_{6,\ 1.3\times 10^5}})$
, (c) kinetic energy for different mesh sizes and (d) relative error in kinetic energy,
$\displaystyle E_R({\textit{KE}}) =({{\textit{KE}} - {\textit{KE}}_{1.3\times 10^5}})/({{\textit{KE}}_{1.3\times 10^5}})$
, for
${\textit{Re}} = 300$
and
${\textit{Bn}} = 2$
.
3. Results and discussion
3.1. The Newtonian limit
Figure 3 shows snapshots of dye concentration in a Newtonian fluid (
${\textit{Bn}}=0$
) at
${\textit{Re}}=100$
. To facilitate the illustration of concentration development, the snapshots show subdomains of different radii (indicated as
$R_{sd}$
in the captions). The white and grey lines show the stirrer’s path and the streamlines, respectively. A normalised variance of the dye concentration,
$\sigma _{20}^2$
, is shown in figure 3(g). The markers in figure 3(g) indicate the time instances of the snapshots.
Figure 3. (a–f) Snapshots of the dye concentration field in a Newtonian fluid,
${\textit{Bn}} = 0$
,
${\textit{Re}}=100$
, at times
${T} = 1, 3, 10, 15, 30 \,\text{and}\,50$
. The white circle marks the stirrer’s path, while grey lines represent the streamlines. The radius of the field of view is provided in the caption. (g) Time evolution of the normalised variance, with circular markers indicating the time instances of the snapshots in (a–f).
When stirring begins, the stirrer crosses the dye interface periodically, stretching and folding it (see, e.g. figures 3
a and 3
b). This action creates a striated pattern near the stirrer’s path, from hereon referred to as the ‘central region’. The stretched interface and thin striations promote local diffusion. As a result, mixing proceeds rapidly during this stage, characterised by a sharp decline in
$\sigma ^2_{20}$
(see figure 3
g). Mixing slows down once the dye concentration becomes nearly uniform along the stirrer’s path (
$T\approx 10$
; see figure 3
c).
For
$10 \lesssim T \lesssim 20$
, the dye concentration appears approximately axisymmetric within the central region, and mixing is temporarily dominated by radial diffusion across the streamlines within this region and across the boundary of this region (figure 3
d). A well-mixed subdomain, where the dye concentration is
$\alpha \approx 0.5$
, thus develops in the central region. However, this region is not quiescent, as the flow is not axisymmetric. Figure 3(e) shows that the well-mixed region is slowly advected by the flow, following an approximately helical trajectory. Consequently, the dye interface is brought back into the stirrer’s path, where it is once again stretched and folded (figure 3
f). The acceleration of mixing during this stage is congruent with the downward concavity of
$\sigma _{20}^2$
for
$T \gtrsim 20$
(figure 3
g).
Another mechanism contributing to mixing during this phase (
$T \gtrsim 20$
) is the extensive stretching of the interface beyond the central region (see figure 3
f). As the well-mixed region drifts away from the centre, interface stretching extends beyond the stirrer’s direct influence. This contributes to the sharp decay rates observed in
$\sigma ^2_{20}$
during this period (
$T \gtrsim 20$
).
The development of the dye concentration is governed by the momentum balance in the flow field; i.e. the coupling between the dye concentration and momentum transfer is one way because the rheology and density are independent of the dye concentration. Figure 4 illustrates the development of the vorticity field in Newtonian fluids at
${\textit{Re}}=100$
, at the same time instances as figure 3. Similar to figure 3, the white circular line is the stirrer’s path and the grey lines are streamlines. To facilitate the illustration of the vorticity field, a scaled vorticity magnitude (
$\zeta$
) is defined as follows:
\begin{align} \zeta = \left \{ \begin{array}{lr} \displaystyle \log (\omega _z + 1) & \text{if} \quad \omega _z \geqslant 0 ,\\[6pt] -\log (|\omega _z| + 1) & \text{if} \quad \omega _z \lt 0. \end{array} \right . \end{align}
Figure 4. (a–f) Snapshots of the vorticity field in a Newtonian fluid,
${\textit{Bn}} = 0$
,
${\textit{Re}}=100$
, at times
${T} = 1, 3, 10, 15, 30\, \text{and}\,50$
. The white circle marks the stirrer’s path, while grey lines represent the streamlines. The radius of the field of view is provided in the caption.
Note that
$\zeta$
retains the sign of
$\omega _z$
; i.e. positive and negative
$\zeta$
indicate counter-clockwise (CCW) and clockwise (CW) rotation, respectively.
When stirring starts, two small attached eddies appear behind the stirrer. More vortices are shed as the stirrer moves along its path (see figure 4
a). The rotation of these vortices is consistent with the stirrer’s movement with CCW and CW vortices shed inside and outside the stirrer’s path, respectively. This is reminiscent of the archetypal problem of the flow around a cylinder of diameter
$\hat {d}_s$
that moves at the constant velocity of
$\hat {U}_o = \hat {r}_o\hat {\varOmega }$
. In the archetypal problem, laminar shedding in a Newtonian fluid of viscosity
$\hat {\mu }$
and density
$\hat {\rho }$
is expected at
$40 \lesssim {\textit{Re}}^* \lesssim 150$
where
${\textit{Re}}^*=\hat {\rho }\hat {U}_o\hat {d}_s/\hat {\mu }$
(Roshko Reference Roshko1954).
Figure 4(b) represents the early dynamics following the onset of stirring. The periodic passage of the stirrer disrupts the formation and advection of CCW vortices, while CW vortices remain outside the stirrer’s path, drifting slowly away (compare figures 4 b and 4 c). Meanwhile, as the attached eddies develop, the CCW eddy expands across the stirrer’s path. When the stirrer moves through this eddy, it generates additional CCW vortices (see figure 4 d). These CCW vortices also dissipate quickly before reaching beyond the central region.
Figure 4(a) also shows that the stirrer completes a period before the previously shed CW vortices advect away from the central region; e.g. during the second period (
$1\leqslant T \leqslant 2$
), the stirrer interacts with the CW vortices shed during the first period. This mechanism interferes with vortex shedding: temporarily (
$2 \leqslant T \lt 20$
), the stirrer’s path is surrounded by two CW vortices that interact with and follow the stirrer slowly moving away from the domain centre. The first two stages of mixing (stretching and folding of the interface followed by diffusion) take place concurrently with the vorticity dynamics discussed above.
At
$T \approx 20$
, the attached CCW eddy expands sufficiently to escape the stirrer’s path (figure 4
e). As this eddy gradually moves away from the central region, the approximate symmetry of the vorticity field, along with the corresponding symmetry in dye concentration in the central region (figure 3
e), is disrupted. From this point onward, vortices are shed away from the central region, travelling far across the flow domain (figure 4
f).
To better illustrate the stages of vortex development, we present an approximation of the locus of the vortex centres in figure 5. These centres, identified as the local minima and maxima of the vorticity field, correspond to the CCW and CW vortices, respectively. For simplicity and brevity, we refer to this collection of points as the vortex centres hereafter. When stirring begins, a pair of attached eddies forms and travels with the stirrer (figure 5
a). Simultaneously, two CW vortices are shed outside the stirrer’s path within the first period (
$T \leqslant 1$
). As shown in figure 5(b) , these shed vortices initially remain close to the central region, drifting away slowly. This quasi-periodic flow pattern leads to an approximately axisymmetric concentration field and a diffusion-dominated mixing phase. Finally, the vortices escape the central region, causing the well-mixed region to advect outward (figure 5
c) – a stage marked by interface stretching and folding beyond the central region. The kinetic energy of escaping vortices decays exponentially due to viscous dissipation, allowing them to advect indefinitely. Consequently, mixing continues unbounded as the vortices transport dye farther into the flow domain.
Figure 5. Time evolution of the vortex centres in a Newtonian fluid,
${\textit{Bn}} = 0$
,
${\textit{Re}}=100$
, over different time intervals, as indicated by the colour bar. Lighter shades correspond to earlier times within each panel. The white solid line represents the stirrer’s path, while the black dashed lines denote the subdomain boundary.
Comparing the development of mixing and the vorticity field, we identify three key mechanisms promoting mixing during different stages of the flow: (i) local stretching and folding of the dye interface with the movement of the stirrer, (ii) diffusion-dominated mixing when flow structure is approximately steady and (iii) the advection of dye and stretching of the interface with vortices that escape the central region.
3.2. Blending a viscoplastic fluid
Figure 6 shows the development of the dye concentration when the fluid has a small yield stress,
${\textit{Re}}=100$
and
${\textit{Bn}}=0.025$
. The dashed grey lines display the contours of
$\tau =1.01Bn$
, a conservative estimate of the boundary of the yielded region. As before, the white and grey solid lines represent the stirrer’s path and the streamlines, respectively.
Figure 6. Snapshots of the dye concentration field in a VPF with a low yield stress,
${\textit{Bn}} = 0.025$
,
${\textit{Re}}=100$
and
$R_{sd} = 10$
. The white circle indicates the stirrer’s path, grey lines represent streamlines and dashed grey lines show contours of
$\tau =1.01Bn$
.
The primary stages of mixing resemble those observed in the Newtonian case, beginning with the stretching and folding of the dye interface (figures 6 a and 6 b). This is followed by the formation of an approximately axisymmetric concentration profile and the emergence of a well-mixed region within the central area (figures 6 c and 6 d, respectively). As the well-mixed region drifts outward from the centre, the interface undergoes extensive stretching, accelerating the mixing process (figures 6 e and 6 f). A comparison between figures 3 and 6 further shows that the yield stress has little influence on the concentration field until the well-mixed region moves beyond the central area.
Figure 7 presents the evolution of the vorticity field at the same time instances as figure 6. Comparison with figure 4 illustrates that the vorticity fields are quite similar until vortices escape the central region. This indicates that, within the central region, the influence of the yield stress is negligible compared with purely viscous effects. However, outside this region, energy dissipation due to yield stress becomes increasingly significant with distance from the domain centre. Far enough from the centre, the fluid is expected to remain unyielded.
In the viscoplastic case, the escaped vortices predominantly drift in the azimuthal direction (see figures 7 e and 7 f). In contrast, in the Newtonian case, the radial displacement of the vortices is more pronounced. This highlights the first mechanism of mixing localisation in VPFs; when the fluid has a yield stress, escaped vortices are confined within a finite distance from the stirrer, effectively localising mixing. This behaviour contrasts with that of Newtonian fluids, where vortices can theoretically advect without bounds.
The effect of the yield stress on the vortex dynamics is further demonstrated in figure 8, which presents the time evolution of the approximate vortex centres over the same intervals as shown in figure 5. The similar development of vortices in the central region is evident in figures 8(a) and 8(b). However, a comparison of figures 5(c) and 8(c), reveals that the azimuthal movement of the vortex centres is more pronounced when
${\textit{Bn}}\gt 0$
.
Figure 7. Snapshots of the vorticity field in a Newtonian fluid, a VPF with a low yield stress,
${\textit{Bn}} = 0.025$
,
${\textit{Re}}=100$
and
$R_{sd} = 10$
. The white circle indicates the stirrer’s path, grey lines represent streamlines, and dashed grey lines show contours of
$\tau =1.01Bn$
.
Figure 8. Time evolution of the vortex centres in a VPF with a low yield stress,
${\textit{Bn}} = 0.025$
,
${\textit{Re}}=100$
, over different time intervals, as indicated by the colour bar. Lighter shades correspond to earlier times within each panel. The white solid line represents the stirrer’s path, while the black dashed lines denote the subdomain boundary. The radius of the field of view is provided in the caption.
Figure 9 illustrates the evolution of dye concentration at a moderate yield stress,
${\textit{Bn}} = 0.4$
. The approximate size of the yielded region is significantly smaller compared with the case of
${\textit{Bn}} = 0.025$
(see dashed grey lines), and mixing remains relatively confined to the central region. The initial two stages of mixing are qualitatively similar to the previous cases: stretching and folding of the interface are followed by diffusion-dominated mixing, leading to the formation of a well-mixed region. However, in this case, the well-mixed region remains largely quiescent. The deformation of the dye interface does not extend beyond the central region.
Figure 9. (a–f) Snapshots of the dye concentration field in a VPF with a moderate yield stress,
${\textit{Bn}} = 0.4$
,
${\textit{Re}}=100$
and
$R_{sd} = 6$
. The white circle indicates the stirrer’s path, grey lines represent streamlines and dashed grey lines show contours of
$\tau = 1.01Bn$
.
The evolution of the vorticity field at
${\textit{Bn}} = 0.4$
is shown in figure 10. The initial shedding pattern appears similar to the Newtonian case (figure 10
a). However, the shed CW vortices stay closer to the stirrer and the CCW eddy remains mostly confined within the stirrer’s path because the yield stress suppresses the advection of shed vortices and growth of attached eddies. Consequently, the size of the well-mixed region is reduced compared with the previous cases. This localisation is closely connected to the development of vortices. Figure 11 illustrates the time evolution of the approximate vortex centres. The formation of the attached eddies and vortex shedding during the first period is similar to the previous cases (figure 11
a). However, the shed vortices do not travel as far from the stirrer (figure 11
b). Instead, they remain trapped in the vicinity of the stirrer’s path (figure 11
c). As a result, mixing remains confined to the central area throughout the process (figure 9
f). This behaviour represents the second mechanism of mixing localisation due to yield stress; as shed vortices remain trapped near the stirrer, the stretching and folding of the interface are confined to this region, limiting the mixing rate.
Figure 10. (a–f) Snapshots of the vorticity field in a VPF with moderate yield-stress value,
${\textit{Bn}} = 0.4$
,
${\textit{Re}}=100$
,
$R_{sd} = 6$
. The white circle indicates the stirrer’s path, while the grey lines depict the streamlines. The dashed grey lines display the contours of
$\tau = 1.01Bn$
.
Figure 11. The time evolution of the vortex centres in the VPF with a moderate yield stress,
${\textit{Bn}} = 0.4$
,
${\textit{Re}}=100$
and
$R_{sd}=3$
, over different time intervals, as indicated by the colour bar. Lighter shades correspond to earlier times within each panel. The white solid line represents the stirrer’s path, while the black dashed lines denote the subdomain boundary.
As expected, at higher-yield-stress values, mixing becomes increasingly confined to the central region. An illustrative case for
${\textit{Bn}}=1$
is shown in figure 12, where the first and second rows display the evolution of dye concentration and vorticity fields, respectively. The dye concentration rapidly assumes an approximately axisymmetric distribution, after which mixing is primarily governed by radial diffusion. The two eddies remain attached to the stirrer, with no vortex shedding observed. This behaviour is further illustrated in figure 13, showing the time evolution of the approximate vortex centres. This represents the third mode of localisation; the complete suppression of vortex shedding, where interface deformation remains largely confined to the central region.
Figure 12. Snapshots of the dye concentration (a–d) and vorticity (e–h) fields in a VPF with a high yield stress,
${\textit{Bn}}=1$
,
${\textit{Re}}=100$
and
$R_{sd}=5$
. The white circle indicates the stirrer’s path, grey lines represent streamlines and dashed grey lines show contours of
$\tau = 1.01Bn$
.
Figure 13. The time evolution of the vortex centres in the VPF with a high yield stress,
${\textit{Bn}}=1$
,
${\textit{Re}}=100$
and
$R_{sd}=3$
. Lighter shades correspond to earlier times. The white solid line represents the stirrer’s path, while the black dashed line denotes the subdomain boundary.
Increasing
${\textit{Bn}}$
beyond
${\textit{Bn}} \sim 1$
does not introduce significant qualitative changes to the mixing dynamics. Figure 14 presents a representative case at
${\textit{Bn}}=100$
(
${\textit{Bn}} \gg 1$
), where the evolution of concentration and vorticity fields are shown in the first and second rows, respectively. Compared with
${\textit{Bn}}=1$
(see figures 12
a and 14
a), interface deformation after the first stirring period is markedly reduced. The stirrer fails to draw the interface deep into the dye-free region because the yielded zone is much smaller than in the
${\textit{Bn}}=1$
case. Indeed, the stirrer’s trajectory now partially traverses unyielded regions. The yielded zone quickly settles into a steady configuration relative to the stirrer. Both the size and shape of the yielded region resemble the asymptotic case of a cylinder translating through a VPF in the high-yield-stress limit (Tokpavi, Magnin & Jay Reference Tokpavi, Magnin and Jay2008; Supekar, Hewitt & Balmforth Reference Supekar, Hewitt and Balmforth2020). Initially, the interface within this region is periodically stretched by successive passages of the stirrer (figures 14
a and 14
b). However, mixing soon becomes diffusion dominated as a well-mixed core develops (figures 14
c and 14
d). Figure 14(e–h) highlights the rapid establishment of a flow field that is steady in the stirrer’s frame and strongly localised around it. Instead of extending along the stirrer’s trajectory, the attached eddies stretch predominantly in the radial direction, which explains the limited elongation of the interface along the path of the stirrer.
Figure 14. Snapshots of the dye concentration (a–d) with
$R_{sd}=3$
and vorticity field (e–h) with
$R_{sd}=2$
, in a VPF with a high yield stress,
${\textit{Bn}}=100$
and
${\textit{Re}}=100$
. The white circle indicates the stirrer’s path, grey lines represent streamlines and dashed grey lines show contours of
$\tau = 1.01Bn$
.
To compare mixing rates across different regimes, we present the normalised variance of dye concentration,
$\sigma _{15}^2$
, in figure 15. For reference, the solid blue line shows the purely diffusive case. As expected, mixing rates decrease with increasing
${\textit{Bn}}$
. At low
${\textit{Bn}}$
(e.g.
${\textit{Bn}}=0.05$
), initially the mixing curves closely match the Newtonian case, indicating that, at low yield stress, while the flow remains confined to the central region (
${T} \lesssim 10$
), yield stress has a negligible effect on local energy dissipation. However, the influence of yield stress becomes more pronounced in later stages, as the well-mixed region drifts outward (see
${T} \approx 20$
and
${T} \approx 30$
for
${\textit{Bn}} = 0.025$
and
${\textit{Bn}} = 0.05$
, respectively). This shift occurs when the dye interface enters a region where the energy dissipation and flow dynamics begin to be affected by the yield stress. From this stage onward, the
$\sigma _{15}^2$
curves diverge, as the advection of escaped vortices is increasingly suppressed by the yield stress.
Figure 15. Evolution of the normalised dye concentration variance for
$0 \leqslant Bn \leqslant 100$
. Darker shades of grey correspond to higher
${\textit{Bn}}$
values; the curves for
${\textit{Bn}}=10$
and
${\textit{Bn}}=100$
nearly overlap and are not visually distinguishable. The solid blue line denotes pure diffusion, and the red dashed lines show exponential fits during the diffusion-dominated stage.
${\textit{Re}}=100$
.
At moderate
${\textit{Bn}}$
, the divergence of
$\sigma _{15}^2$
from the Newtonian case becomes apparent earlier. This is because, as
${\textit{Bn}}$
increases, the dissipation caused by yield stress starts to influence smaller radii. The downward concavity in the variance is almost entirely suppressed, as vortices remain confined within the central region, preventing them from re-accelerating mixing.
The mixing rate decreases progressively with increasing
${\textit{Bn}}$
, reflecting both the overall slowdown of the flow and the suppression of vortices. For
${\textit{Bn}} \gtrsim 1$
, however, it shows little further variation. To assess these changes, we measure the mixing rate during the final, diffusion-dominated stage, where the concentration variance follows an exponential decay
\begin{align} \sigma _{15}^2 \approx \sigma _0^2 \exp \left (-2\lambda {T}\right ) .\end{align}
Here,
$\lambda$
is the decay constant. The red dashed lines in figure 15 depict the exponential fits. Following Christov & Homsy (Reference Christov and Homsy2009), we define the enhancement factor
$\eta _{\lambda } = \lambda /\lambda _p$
, where
$\lambda _p$
is the decay constant for the purely diffusive case. The intercept of the exponential fit,
$\sigma _0$
, can be viewed as an approximation of the standard deviation when the diffusion-dominated phase begins, providing a measure of mixing by the end of the advection-dominated stage. Here, a second enhancement factor is defined as
$\eta _\sigma = \sigma _{op}/\sigma _0$
, where
$\sigma _{op}$
is the intercept of the exponential fit for the purely diffusive case. This enhancement factor,
$\eta _\sigma$
, enables comparison of the extent to which mixing is achieved during the initial phase dominated by interface stretching.
Figure 16. Variation of the enhancement factors (a)
$\eta _\lambda$
and (b)
$\eta _\sigma$
with
${\textit{Bn}}$
at
${\textit{Re}}=100$
.
The variation of enhancement factors with
${\textit{Bn}}$
is shown in figure 16. Changes in
$\eta _\lambda$
are more pronounced than those in
$\eta _\sigma$
. This suggests that the long-term impact of yield stress on mixing can be more significant than its short-term effects; increasing the Bingham number by approximately three orders of magnitude reduces
$\eta _\sigma$
by approximately
$10\,\%$
, attributed to the suppression of dye interface stretching during the initial mixing stage. Meanwhile, the same increase in
${\textit{Bn}}$
results in an approximately
$80\,\%$
reduction in
$\eta _\lambda$
, as diffusion – the dominant mixing mechanism over longer time scales – becomes less effective when the interface area is reduced. Lastly, both enhancement factors drop relatively sharply for small
${\textit{Bn}}$
, but exhibit only a slow decline beyond
${\textit{Bn}} \sim 1$
. This behaviour is discussed in greater detail in § 3.3.
3.3. Flow and mixing regimes
In the previous section, we identified three distinct mixing regimes:
Regime SE (shedding, escaped vortices): at sufficiently low
${\textit{Bn}}$
, mixing is initially dominated by the stretching and folding of the interface. This is followed by a brief period where mixing is driven primarily by radial diffusion. In this regime, vortices escape the central region, carrying dye and stretching the interface deep into the domain. Once the vortices escape, mixing re-accelerates, driven by advection and further interface stretching.
Regime ST (shedding, trapped vorticed): at moderate
${\textit{Bn}}$
, mixing does not re-accelerate after the diffusion-dominated phase. In this regime, the shed vortices remain trapped in the central region, unable to escape. They are periodically influenced by the passage of the stirrer but remain confined, limiting further advection-driven mixing.
Regime NS (no shedding): at sufficiently high
${\textit{Bn}}$
, mixing is quickly dominated by diffusion. The yield stress completely suppresses vortex shedding, confining interface stretching to the immediate vicinity of the stirrer.
The mixing regimes described above are closely tied to the flow development and vortex dynamics. To facilitate comparison and distinction of the regimes, we examine the time evolution of kinetic energy (approximated by the velocity norm,
$||u||$
) and the size of the moving region. For the latter, we consider the general condition illustrated in figure 17. Here,
$r_{\textit{yo}}$
and
$r_{yi}$
mark the farthest and closest points where the line connecting the domain centre to the stirrer-path centre crosses the boundary of the quiescent unyielded region.
Figure 17. Schematic of the moving-region boundary (dashed grey line), with
$r_{yi}$
(red) and
$r_{\textit{yo}}$
(blue) marking the nearest and farthest intersections of the line from the domain centre to the stirrer-path centre with the boundary of the quiescent unyielded region.
When the fluid is Newtonian (
${\textit{Bn}}=0$
),
$r_{yi}=0$
and
$r_{\textit{yo}}\rightarrow \infty$
. On the other hand, if the fluid has a yield stress (
${\textit{Bn}}\gt 0$
), then
$r_{\textit{yo}}\lt \infty$
. Furthermore, the results presented above (e.g.
${\textit{Bn}}=0.4$
shown in figure 9) illustrate that, at sufficiently small
${\textit{Bn}}$
,
$r_{yi}=0$
. Nevertheless, at sufficiently large
${\textit{Bn}}$
(see, e.g.
${\textit{Bn}}=100$
shown in figure 14),
$r_{yi}\gt 0$
.
Figure 18(a) shows the time evolution of
$r_{\textit{yo}}$
at
${\textit{Re}}=100$
. In regimes SE and ST,
$r_{\textit{yo}}$
evolves on two distinct time scales: a shorter time scale (
$T_s \approx O(1)$
), corresponding to small oscillations in
$r_{\textit{yo}}$
, and a much longer time scale (
$T_l \approx O(10)$
), associated with the overall increase in
$r_{\textit{yo}}$
. The amplitude of the oscillations diminishes as
${\textit{Bn}}$
increases, eventually disappearing entirely in regime NS. The white dashed lines in figure 18(a) display the moving time average (
$\overline {r}_{\textit{yo}}$
) with an averaging window of
$T=1$
, which better illustrates variations over longer time scales.
Figure 18. (a) Time evolution of the outer radius of the yielded region,
$r_{\textit{yo}}$
, at different yield-stress values. Darker grey shades correspond to larger
${\textit{Bn}}$
. The white dashed lines display the moving time average (
$\overline {r}_{\textit{yo}}$
). (b) Time evolution of the velocity norm at different yield-stress values. The inset provides a magnified view over the range
$40\lt T\lt 44$
for regimes ST and NS.
${\textit{Re}}=100$
.
Immediately after stirring begins,
$r_{\textit{yo}}$
rises sharply, reflecting the onset of yielding. The oscillations in
$r_{\textit{yo}}$
observed in regimes SE and ST (e.g.
${\textit{Bn}}=0.025$
and
${\textit{Bn}}=0.4$
) are linked to vortex interactions, indicating that the moving region remains unsteady relative to the stirrer. By contrast, in regime NS, vortex shedding and escape from the central region do not occur; instead,
$r_{\textit{yo}}$
rapidly converges to a steady value, suggesting that the flow is steady in the stirrer’s frame of reference. This steady value depends only weakly on
${\textit{Bn}}$
: increasing
${\textit{Bn}}$
by two orders of magnitude reduces
$r_{\textit{yo}}$
by less than
$50\,\%$
.
At very large
${\textit{Bn}}$
, the flow configuration and moving region resemble the asymptotic case of a cylinder translating through a VPF in the high-yield-stress limit. For instance, at
${\textit{Bn}}=100$
the moving region corresponds to
$r_{\textit{yo}}\approx 1.86$
, only slightly larger than the plastic limit of
$r_{\textit{yo}}\approx 1.70$
for uniform flow past a cylinder (Supekar et al. Reference Supekar, Hewitt and Balmforth2020). This explains the very slow decrease of the mixing rate at
${\textit{Bn}}\gtrsim 1$
(see figure 16); in regime NS, the moving region very slowly approaches the plastic limit with
$\displaystyle \lim _{Bn\rightarrow \infty } r_{\textit{yo}}\sim 1.70$
. At high
${\textit{Bn}}$
, the radial extent of interface deformation is thus restricted, while stretching along the stirrer’s path diminishes with increasing
${\textit{Bn}}$
. We hypothesise that these two effects together account for the very slow decay of mixing efficiency when
${\textit{Bn}} \gtrsim 1$
(see figure 16).
In regime SE,
$r_{\textit{yo}}$
increases to a local maximum as the CW shed vortices travel away from the stirrer, orbiting the stirrer’s path (see
$T\approx 5$
and
${\textit{Bn}}=0.025$
in the figure).
$r_{\textit{yo}}$
decreases as these satellite vortices spread azimuthally and weaken. The final increase in
$r_{\textit{yo}}$
corresponds to the escape of vortices from the central region (e.g.
$T\approx 20$
for
${\textit{Bn}}=0.025$
). Although we expect
$r_{\textit{yo}}$
to approach an upper bound for
${\textit{Bn}}\geqslant 0$
, numerical estimation of this ultimate threshold in regime SE was not feasible within the computational domain used in this study. However, the approach of
$r_{\textit{yo}}$
to a steady average value,
$\overline {r}_{ys}$
, is clear in regimes ST and NS (see red dashed lines in the figure).
Finally, in regime ST, the relatively slow increase in
$r_{\textit{yo}}$
during the early stages (see
${\textit{Bn}}=0.4$
,
$T\lesssim 10$
) confirms that the yield stress delays the advection of shed vortices away from the stirrer. Additionally, no further increases in
$r_{\textit{yo}}$
are observed, as the vortices remain confined within the central region.
Figure 18(b) shows the time evolution of the kinetic energy. As with
$r_{\textit{yo}}$
, the kinetic energy exhibits two characteristic time scales: a short time scale (
$T_s \approx O(1)$
), corresponding to small oscillations in
$\textit{KE}$
, and a longer time scale (
$T_l \approx O(10)$
), reflecting the overall growth of
$\textit{KE}$
. The oscillation amplitude again diminishes with increasing
${\textit{Bn}}$
, disappearing entirely in regime NS. Here,
$\textit{KE}$
increases rapidly when stirring begins, before reaching an apparently quasi-steady state (see, e.g.
${\textit{Bn}}=0.05$
,
$10\lesssim {T} \lesssim 25$
). This corresponds to the time interval when the shed CW vortices advect within the central region. In regime SE, this stage is followed by a second, relatively rapid increase, corresponding to the escape of vortices from the central region (see, e.g.
${\textit{Bn}}=0.05$
,
$30 \lesssim {T} \lesssim 45$
). This transition is not observed in regime ST, where shed vortices remain trapped in the central region (see, e.g.
${\textit{Bn}}=0.4$
in the figure).
Two critical values of
${\textit{Bn}}$
may be identified to distinguish the three regimes.
\begin{align} \left \{ \begin{array}{lll} \text{Regime SE} & \text{if}& Bn \leqslant {\textit{Bn}}_{\textit{ce}} ,\\[5pt] \text{Regime ST} & \text{if} & {\textit{Bn}}_{\textit{ce}} \leqslant \textit{Bn} \leqslant {\textit{Bn}}_{\textit{ct}},\\[5pt] \text{Regime NS} & \text{if} & {\textit{Bn}}_{{\textit{ct}}} \leqslant Bn ,\end{array} \right . \end{align}
where
${\textit{Bn}}_{\textit{ce}}$
and
${\textit{Bn}}_{{\textit{ct}}}$
are the critical Bingham numbers marking the regime transitions.
To differentiate the regimes, we evaluated the Fourier spectrum of the oscillations,
${\textit{KE}}^{\prime}={\textit{KE}}-\overline {{\textit{KE}}}$
, where
$\overline {{\textit{KE}}} = \int _{{T}}^{{T}+1} {\textit{KE}} \, {\rm d}{T}$
. Figure 19 shows the frequency spectrum of
${\textit{KE}}^{\prime}$
at
${\textit{Re}}=100$
,
$0 \leqslant Bn \leqslant 100$
. All regimes have peaks at frequency one and its harmonics. This is the fundamental frequency of the system for all
${\textit{Bn}}$
considered. The prominence of the peaks diminishes as
${\textit{Bn}}$
increases. Regimes SE and ST have additional distinct peaks (see
$f\approx 1.6$
and
$1.8$
, respectively, in figure 19). These characteristic peaks vanish as mixing transitions between regimes. In contrast, regime NS does not display any peaks beyond the fundamental frequency and the harmonics. To estimate the critical values, we evaluated the strength of the characteristic peak for each regime and used the two closest data points to find, by extrapolation, the
${\textit{Bn}}$
value at which the characteristic peak disappears.
Figure 19. Fourier spectrum of
${\textit{KE}}^{\prime}$
. Here,
${\textit{Re}}=100$
.
Figure 20 illustrates the flow regime map in the (
${\textit{Re}},Bn$
) plane. The orange, green and blue zones correspond to regimes NS, ST and SE, respectively. Regimes SE and ST disappear at sufficiently small
${\textit{Re}}$
below which there is no shedding even in the Newtonian case. The square and circle markers indicate the estimates of
${\textit{Bn}}_ce$
and
${\textit{Bn}}_{\textit{ct}}$
, respectively.
Figure 20. Flow regime map in the (
${\textit{Re}},Bn$
) plane. Regimes NS, ST and SE, are indicated in orange, green and blue, respectively. The triangle markers are the data points. Square and circle markers display the estimates of
${\textit{Bn}}_{\textit{ce}}$
and
${\textit{Bn}}_{\textit{ct}}$
, respectively. The black dotted and dashed lines are linear fits to
${\textit{Bn}}_{\textit{ce}}$
and
${\textit{Bn}}_{\textit{ct}}$
.
The dashed and dotted lines in figure 20 illustrate that the variations of
${\textit{Bn}}_{\textit{ct}}$
and
${\textit{Bn}}_ce$
are well described by the following fits:
\begin{align} \begin{array}{l} {\textit{Bn}}_\textit{ce} = 0.00287 {\textit{Re}} - 0.14, \\ {\textit{Bn}}_{\textit{ct}} = 0.0111 {\textit{Re}} - 0.53 .\end{array} \end{align}
The linear variation of the critical Bingham numbers with
${\textit{Re}}$
is reminiscent of the effective Reynolds number,
${\textit{Re}}_e$
, in flows of VPFs
\begin{align} {\textit{Re}}_{e} = \frac {\hat {\rho } \hat {U}_o^2}{\hat {\tau }_y+\hat {\mu } \hat {\dot {\gamma }}_c} ,\end{align}
where
$\hat {\dot {\gamma }}_c$
is the characteristic strain rate. The primary challenges in identifying the characteristic strain rate in flows of VPFs are twofold: firstly, the physically representative
$\hat {\dot {\gamma }}_c$
may vanish as yield stress increases and the fluid becomes quiescent. Secondly, the relationship between
$\hat {\dot {\gamma }}_c$
and yield stress is not known a priori; see Thompson & Soares (Reference Thompson and Soares2016) and Ahmadi, Olleik & Karimfazli (Reference Ahmadi, Olleik and Karimfazli2022) for more details.
In this context, the linear critical equations presented in (3.4) reveal the effective Reynolds numbers at which regime transitions occur. Assume that
$\hat {U}_o= \hat {r}_o\hat {\varOmega }$
and
$ \hat {\dot {\gamma }}_c = a \hat {\varOmega }$
, where
$a$
is a constant, then (3.5) yields
\begin{align} & {\textit{Re}}_e = \frac {{\textit{Re}}}{a+Bn} .\end{align}
Comparing this with (3.4), we find
\begin{align} \begin{array}{lll} &\displaystyle {\textit{Re}}_{e,ce} = \frac {{\textit{Re}}}{{\textit{Bn}}_{{\textit{ct}}} + 0.53} \approx 350, & a_{ce}\approx 0.14, \\[3pt] &\displaystyle {\textit{Re}}_{e,{\textit{ct}}} = \frac {{\textit{Re}}}{{\textit{Bn}}_{\textit{ce}} + 0.14} \approx 90, & a_{\textit{ct}}\, \approx 0.53 .\end{array} \end{align}
In the limit corresponding to Newtonian fluids (
${\textit{Bn}} \rightarrow 0$
), the effective Reynolds number is given by
${\textit{Re}}_e = {\textit{Re}}/a$
. On the other hand, as
$c \rightarrow \infty$
, the critical Reynolds number beyond which the wake behind a cylinder in a free stream of Newtonian fluids becomes unstable is (Roshko Reference Roshko1954)
\begin{align} {\textit{Re}}_{\textit{cr}}^* = \frac {{\textit{Re}}_{\textit{cr}}}{c} \approx 40 ,\end{align}
which implies that, for
$c \rightarrow \infty$
,
${\textit{Re}}_{\textit{cr}} \approx 80$
, or alternatively
${\textit{Re}}_{e,{\textit{ct}}} \approx {\textit{Re}}_{\textit{cr}}/a_{\textit{ct}} \approx 151$
, marks the onset of wake instability. Comparing these values with
${\textit{Re}}_{e,{\textit{ct}}}$
when
$c=2$
(see 3.7) suggests that the wake behind the cylinder becomes more unstable as the path curvature decreases.
The two distinct definitions of the effective Reynolds number shown in (3.7) confirm that, for a given flow set-up, the characteristic strain rate changes with
${\textit{Bn}}$
. It follows that defining a unique effective Reynolds number that fully describes the hydrodynamics observed in the (
${\textit{Re}},Bn$
) plane, e.g. help distinguish the different regimes at different values of the Bingham number, is not feasible.
To quantitatively compare the extent of localisation in different regimes, figure 21(a) displays the steady time-averaged radius of the moving region,
$\overline {r}_{ys}$
, in regimes ST and NS for a wide range of
${\textit{Re}}$
. Distinct colours and marker shapes represent various
${\textit{Re}}$
values. Empty and filled markers indicate regimes ST and NS, respectively.
Figure 21. (a) Steady values of the time-averaged radius of the yielded region,
$\bar {r}_{ys}$
, at different
${\textit{Bn}}$
and
${\textit{Re}}$
. Each symbol (with unique colour) corresponds to a specific
${\textit{Re}}$
; empty and filled symbols denote regimes ST and NS, respectively. (b) Evolution of
$\bar {r}_{ys}$
in regime ST as a function of the modified Bingham number,
$\displaystyle {\textit{Bn}}_m = ({Bn-{\textit{Bn}}_{\textit{ct}}})/({{\textit{Bn}}_{\textit{ce}}-{\textit{Bn}}_{\textit{ct}}})$
.
In regime
$\textrm{ST}$
,
$\bar {r}_{ys}$
decreases with decreasing
${\textit{Re}}$
and increasing
${\textit{Bn}}$
, confirming the influence of both yield stress and purely plastic stresses on the flow field. Figure 21(b) shows that, in this regime,
$\bar {r}_{ys}$
values approximately collapse onto the same curve when plotted against a scaled Bingham number,
${\textit{Bn}}_m$
,
\begin{align} {\textit{Bn}}_m=\frac {{\boldsymbol{Bn}}-{\boldsymbol{Bn}}_{\textit{ct}}}{{\boldsymbol{Bn}}_{\boldsymbol{ce}}-{\boldsymbol{Bn}}_{\textit{ct}}} .\end{align}
Additionally, figure 21(a) reveals that
$\bar {r}_{ys}$
is approximately independent of
${\textit{Re}}$
in regime NS. Furthermore, in line with the limited changes in the mixing rate at high Bingham (
${\textit{Bn}} \gtrsim 1$
),
$\bar {r}_{\textit{yo}}$
changes very slowly with
${\textit{Bn}}$
within this range.
Finally, figure 22 illustrates the variation of enhancement factors with
${\textit{Bn}}$
across different
${\textit{Re}}$
. Different values of
${\textit{Re}}$
are distinguished by different marker shapes. Empty and filled symbols indicate regimes ST and NS, respectively. The initial rapid decrease in
$\eta _\lambda$
and
$\eta _\sigma$
, similar to the trend observed in figure 16, corresponds to regime ST. The enhancement factors show little variation in regime NS. This suggests that the suppression of advected shed vortices plays a significantly larger role in mixing localisation and reduction than the contraction of the sheared layer around the stirrer, the sole localisation mechanism active in regime NS. Moreover, the sharpest decline occurs near the transition to regime NS, indicating that the complete suppression of shedding marks a key threshold in the localisation process.
Figure 22. Variation of enhancement factors, (a)
$\eta _{\lambda }$
and (b)
$\eta _{\sigma }$
with
${\textit{Bn}}$
and
${\textit{Re}}$
. Each symbol with unique colour shows a specific
${\textit{Re}}$
. The empty and filled symbols indicate regimes ST and NS, respectively.
It is also evident that the influence of a small yield-stress value is more pronounced at lower
${\textit{Re}}$
; for instance, comparing
${\textit{Re}}=100$
and
${\textit{Re}}=300$
in the figure. At higher
${\textit{Re}}$
, the influence of
${\textit{Bn}}$
on the enhancement factors is less pronounced at low
${\textit{Bn}}$
, and a very sharp decline in the enhancement factors is observed near
${\textit{Bn}} \lesssim {\textit{Bn}}_{\textit{ct}}$
.
4. Summary
In this study, we investigated a canonical two-dimensional mixing set-up to establish a mechanistic understanding of how fluid mechanics drives transitions in the mixing regimes of yield-stress fluids.
Our numerical simulations model an infinite, two-dimensional domain filled with a quiescent VPF described by the Bingham model. In this set-up, a cylinder moves at a constant speed along a circular path, stirring a fluid with uniform density and rheological properties. The stirrer diameter is fixed at half the stirring radius throughout the study. The bottom half of the domain is initially marked with a passive dye, with the centre of the stirrer’s path aligned with the dye interface. This representative model allowed us to explore the fundamental features of yield-stress fluid mixing in two-dimensional settings.
To decouple the influence of the flow dynamics from dye concentration, we considered only one-way coupling, where the flow is unaffected by dye concentration changes. This approach simplifies the isolation of causal relationships between fluid mechanics phenomena and mixing events.
In the Newtonian case, we identified three primary mixing mechanisms when stirring a fluid with heterogeneous dye distribution: (i) the stretching and folding of the dye interface within the central region as the stirrer initiates movement, (ii) diffusion-dominated mixing when dye distribution becomes approximately uniform along streamlines and (iii) enhanced mixing due to vortex shedding, where shed vortices transport dyed regions into dye-free areas (or vice versa), significantly extending the dye interface.
For the laminar regimes studied, the fluid dynamics evolves on two time scales: the average energy of the system changes on a slow time scale, roughly an order of magnitude slower than the stirrer period, while energy oscillations, related to vortex interactions, occur on a faster time scale. As yield stress increases, the average yielded region size, kinetic energy and energy oscillations decrease, eventually leading to oscillation suppression. We identified three mechanisms by which yield stress localises mixing.
-
(i) Finite vortex advection: in yield-stress fluids, advection of shed vortices is suppressed, with vortices travelling within a finite radius from the stirrer due to energy decay. This defines a maximum distance for dye transport and leads to localised mixing.
-
(ii) Entrapment of shed vortices: at moderate yield stresses, shed vortices cannot escape the stirrer’s vicinity, resulting in periodic interactions between the stirrer and previously shed vortices, which promote mixing localisation.
-
(iii) Suppression of vortex shedding: at high yield stresses, vortex shedding ceases entirely, confining mixing to the interface stretching caused by direct interaction between the stirrer and the dye interface.
We classified the observed mixing regimes based on these mechanisms: in regime SE, shed vortices escape the central region, causing dye variance evolution similar to that in Newtonian fluids, with mixing initiated by interface stretching and folding, followed by a diffusion-dominated phase and then acceleration by escaping eddies. Regime ST is defined by the entrapment of vortices near the stirrer, limiting mixing to its immediate vicinity, while regime NS shows no vortex shedding, with mixing characterised by initial interface stretching followed by diffusion across streamlines and across the boundaries of the well-mixed region.
Using a fast Fourier transform on energy oscillations, we distinguished these regimes, all featuring a fundamental mode at frequency one. Additional spectral peaks in regimes SE and ST, introduced by shed vortices, further enabled us to identify critical transition criteria. In the
$({\textit{Re}},Bn)$
plane, we observed distinct separation of the three regimes along lines that relate critical Bingham numbers to Reynolds numbers. This relationship allows us to define two unique, effective Reynolds numbers, each capturing a transition between two regimes at a constant value. This supports the hypothesis that fluid mechanics phenomena underlying mixing regime transitions are closely linked to the bluff-body flow dynamics, traditionally described by threshold Reynolds numbers marking stability and flow transition modes.
Comparing the decay of dye concentration variance,
$\sigma ^2$
, across regimes, we hypothesise that, among the various localisation mechanisms, vortex entrapment near the stirrer has the most significant impact on mixing. This is supported by the contrasting long-term decay behaviours observed in regimes SE and ST: in regime ST, where vortices remain entrapped,
$\sigma ^2$
shows an exponential decay, while in regime SE, it decays at an accelerating rate due to the influence of vortices that escape into the outer regions. To compare regimes ST and NS, we introduced an enhancement factor,
$\eta _\lambda$
, which quantifies the relative decay rate of
$\sigma ^2$
during the final, diffusion-dominated stage of mixing, compared with the purely diffusive case. In regime ST,
$\eta _\lambda$
declines sharply with increasing
${\textit{Bn}}$
, but in regime NS, it becomes largely independent of both
${\textit{Re}}$
and
${\textit{Bn}}$
. This suggests that vortex entrapment is indeed the primary driver of mixing localisation. Moreover, once vortex shedding is fully suppressed (in regime NS), increasing stirrer speed has minimal effect on enhancing mixing.
Finally, mixing localisation has far-reaching effects beyond the initial localisation of interface stretching: while the change in the degree of mixing by the end of the advection-dominated phase may not be very significant (e.g. approximately 10 % for an increase of three orders of magnitude in
${\textit{Bn}}$
), the subsequent change in the decay rate of
$\sigma ^2$
can be much more significant (more than 80 %).
The fundamental principle of laminar mixing in stirred tanks is the periodic movement of a stirrer along a closed path. In this work, we model this process in its most simplified form: a single stirrer moving along a circular path across an initial concentration gradient, within an effectively infinite domain. The dynamics is governed by
${\textit{Re}}$
and
${\textit{Bn}}$
and
$c$
. For a fixed stirring radius, we observe a gradual transition from attached eddies to vortex shedding, with vortices either remaining trapped around the stirrer or escaping the central region.
The vortex dynamics identified in this study is closely linked to the classical problem of vortex shedding behind bluff bodies. The Strouhal number associated with vortex shedding behind bluff bodies, defined as
$St = \hat {f}\hat {d}_s / (\hat {r}_o\hat {\varOmega })$
, increases monotonically with
${\textit{Re}}$
and asymptotically approaches a maximum value,
$St_m$
. Here,
$\hat {f}$
is the shedding frequency. In contrast,
$St$
is expected to decrease with increasing
${\textit{Bn}}$
, as the effective Reynolds number is reduced. Accordingly, variations in
$c$
are expected to produce more complex vortex structures, since the maximum number of vortices shed per stirrer period scales with
$cSt_m$
. A detailed investigation of the influence of
$c$
lies beyond the scope of the present study. Nevertheless, we hypothesise that, in the absence of wall effects, the primary mixing modes and yield-stress-driven localisation mechanisms identified here remain relevant. Qualitatively, the archetypal regimes fall into three categories:
$\textrm{NS}$
(no shedding),
$\textrm{ST}$
(shedding with trapped vortices), and
$\textrm{SE}$
(shedding with escaping vortices).
For passive tracers, the transitions between the mixing regime observed here correlate with distinct effective Reynolds numbers based on the stirrer size, shape and speed, alongside the Bingham number. Alternatively, the critical Bingham numbers marking these transitions are expected to scale linearly with the Reynolds number based on purely viscous stresses. Given the lack of an a priori definition for effective Reynolds numbers in such flows, we recommend further investigation of the linear critical criteria observed here.
Funding
We gratefully acknowledge financial support from Natural Sciences and Engineering Research Council (NSERC), Canada, through the Discovery Grant program. This research was enabled in part by support provided by Calcul Québec (www.calculquebec.ca) and Compute Canada (www.computecanada.ca).
Declaration of interests
The authors report no conflict of interest.
Open access
