1. Introduction
The state-of-the-art development of microfluidic devices has gained attention from both academia and industry because of its widespread applications in biomedicine and biochemistry (Verpoorte & De Rooij Reference Verpoorte and De Rooij2003; Stone, Stroock & Ajdari Reference Stone, Stroock and Ajdari2004). It is worth mentioning that, unlike conventional procedures followed in the biochemical and biomedical industries, miniaturization techniques provide critical advantages such as reduced sample volumes, rapid processes, portability and a high surface-area-to-volume ratio (Djenidi & Moghtaderi Reference Djenidi and Moghtaderi2006; Das et al. Reference Das, Banik, Chen, Sinha and Mukherjee2015). It is quite common for biofluids such as blood, synovial fluid, serum, saliva, nasal fluid, DNA solution, etc., which typically exhibit non-Newtonian behaviour, to experience laminar flow within such small fluidic channels. Notably, effective mixing, quick transportation, rapid testing and efficient diagnostics of samples/solutes with non-Newtonian fluid rheology are necessary for many biochemical processes and pathological diagnostics (Cosentino et al. Reference Cosentino, Madadi, Vergara, Vecchione, Causa and Netti2015; Huang et al. Reference Huang, Chen, Wong and Liow2016; Gaikwad, Mondal & Wongwises Reference Gaikwad, Mondal and Wongwises2018; Herreman et al. Reference Herreman, Nore, Cappanera and Guermond2021). As reported in the literature, researchers have employed both elastic and inelastic models to mimic the rheological and transport characteristics of such biofluids (Vagner & Patlazhan Reference Vagner and Patlazhan2019; Nwani et al. Reference Nwani, Merhaben, Gates and Benneker2022). It is important to add here that significant efforts have been made towards developing effective micromixers through geometrical modulation of flow configurations and the applications of external fields (Alekseenko et al. Reference Alekseenko, Kuibin, Okulov and Shtork1999; Krishnaveni et al. Reference Krishnaveni, Renganathan, Picardo and Pushpavanam2017; Mehta, Pati & Mondal Reference Mehta, Pati and Mondal2021; Shyam, Mondal & Mehta Reference Shyam, Mondal and Mehta2021). In order to achieve proper mixing in on-chip microfluidic devices, simultaneously ensuring effective fluidic transport operations therein, researchers have developed numerous methodologies consistent with either passive methods or active techniques; all while articulating the underlying processes (Cetegen & Mohamad Reference Cetegen and Mohamad1993; Djenidi & Moghtaderi Reference Djenidi and Moghtaderi2006; Hadigol et al. Reference Hadigol, Nosrati, Nourbakhsh and Raisee2011; Rezk et al. Reference Rezk, Qi, Friend, Li and Yeo2012; Chen et al. Reference Chen, Lv, Wei and Li2022).
To achieve effective and efficient advective mixing through passive methods, it is necessary to generate vorticial flow in the mixing process (Afzal & Kim Reference Afzal and Kim2014; Matsunaga & Nishino Reference Matsunaga and Nishino2014; Madana & Ashraf Ali Reference Madana and Ashraf Ali2020). However, such vortices tend to dissipate unless they are continuously sustained by the flow configuration/geometry. The dissipation of vortex or swirl has captivated researchers over the past century. Establishing swirl flow in the fluidic pathway requires implementation of a conducive flow configuration to instigate a swirling motion in the fluid. It may be mentioned here that a swirl flow environment can be developed through either a twisted tape set-up or by utilizing a swirl generator. In the much-celebrated experimental study by Taylor (Reference Taylor1950), it was shown that the employment of a swirl generator at the inlet of the chosen fluidic pathway facilitates tangential fluid injection into the channel, which, in turn, leads to the generation of a primary vortex with a characteristic shape, known as the Rankine vortex (Reader-Harris Reference Reader-Harris1994). Important to mention is that this swirling flow indicates the presence of a tangential velocity along the axial flow direction (Parchen & Steenbergen Reference Parchen and Steenbergen1998).
Early works in this domain primarily focused on investigating the decay of swirl in pipes, both through theoretical and experimental approaches (Stokes et al. Reference Stokes, Graham, Lawson and Boger2001a,Reference Stokes, Graham, Lawson and Bogerb; Binagia et al. Reference Binagia, Phoa, Housiadas and Shaqfeh2020). In the early 1950s, Talbot (Reference Talbot1954) conducted theoretical and experimental studies on laminar swirling flows and revealed that the experimental swirl decay rates closely matched the theoretical predictions. Additionally, an unsteady theoretical flow analysis by Sibulkin (Reference Sibulkin1962) was conducted for the vortex tube, and it was observed that the radial distributions of velocity and temperature at different axial locations qualitatively agreed with experimental results reported by Lay (Reference Lay1959). Subsequently, researchers examined the decay of swirl (Kreith & Sonju Reference Kreith and Sonju1965; Steenbergen & Voskamp Reference Steenbergen and Voskamp1998) in fully developed turbulent flow environments through their analytical solution and identified the dependence of swirl velocity on Reynolds number (Re) and axial location. The decay in swirl, specified by swirl intensity, was found to follow an exponential trend, attributed to changes in parameters such as the wall friction factor, Reynolds number, transition radius and inlet swirl number (Morton Reference Morton1969; Kitoh Reference Kitoh1991; Yao & Fang Reference Yao and Fang2012; Kumar, Shakya & Kaushik Reference Kumar, Shakya and Kaushik2020; Kumar et al. Reference Kumar, Gaikwad, Kaushik and Mondal2023). Numerical simulations, as demonstrated by Kiya, Fukusako & Arif (Reference Kiya, Fukusako and Arif1971), also assumed a Rankine vortex as the inlet swirl velocity profile. This assumption was based on experimental studies (Alekseenko et al. Reference Alekseenko, Kuibin, Okulov and Shtork1999) which showed that the swirl velocity in tubes is essentially a combination of the forced and free vortex components (Reader-Harris Reference Reader-Harris1994). More recently, the decay of swirl intensity in laminar flow of a non-Newtonian fluid was found to be influenced by the rheological nature of the fluid (Kathail, Pranav & Kaushik Reference Kathail, Pranav and Kaushik2019).
However, recent developments in micro/mini fluidic applications within the medical industries have sparked a growing interest in the enhancement of transport and mixing of non-Newtonian fluids. We believe that harnessing the swirl velocity as a passive mixing technique for non-Newtonian fluids holds promise, particularly in laminar flow conditions. Therefore, in this current study, we aim to investigate the fundamental physics of decaying swirl and its influence on the mixing biofluids. Specifically, we seek to analyse how fluid rheology, imposed by vorticial (swirl) flow within a narrow cylindrical tube, affects the complete mixing process between two similar inelastic non-Newtonian fluids. To achieve this, we first present an analytical solution for swirl velocity profiles and swirl decay, laying the groundwork for our investigation. Then, we utilize the analytical swirl velocity solution to numerically solve the species transport equation using an in-house finite volume-based code. With the results in hand, we discuss both the qualitative and quantitative aspects of mixing between two fluids in the context of decaying swirl flow in laminar regime.
2. Problem formulation and mathematical model
In this study, we execute a theoretical analysis to obtain the analytical expression for velocity fields of a non-Newtonian fluid, whose rheology is represented by the Ostwald–de-Waele power-law model. We consider the underlying flow through a narrow fluidic cylindrical pipe, driven by the inlet swirl (Taylor Reference Taylor1950). On using the analytically obtained velocity fields, we solve the species transport equation essentially for the concentration distribution in the chosen fluidic pathway.
2.1. Flow configuration: geometry and description
We consider the transport of non-Newtonian fluids through a cylindrical microfluidic pipe of radius R and length 
$L$, as shown in figure 1. The coordinate system describing the flow with coordinates r, 
$\theta $ and 
$z$, corresponding to velocity fields 
${u_r}$, 
${u_\theta }$ and 
${u_z}$ respectively, is attached to the inlet (left centre) of the channel. The inlet cross-sections of the entrances of fluids A and B are designated by 
$\theta = [0,{\rm \pi} ]$ and 
$[{\rm \pi} ,2{\rm \pi} ]$, respectively. It may be mentioned here that the chosen configuration seems to mimic an application of blood flow getting infused with other biological species/reagents or the case of a DNA carrier fluid being mixed with the chemical reagents or fluorescent tracer particles (Cosentino et al. Reference Cosentino, Madadi, Vergara, Vecchione, Causa and Netti2015). The concentrations of these fluids, such as blood or DNA carrier fluids, are identified by zero concentration 
$({C_0} = 0)$, whereas the concentration of tracer particles or chemical reagents is identified by higher concentration 
$({C_0} = 1)$ (see figure 1). The selection of such similar kinds of fluids also justifies the assumption of considering constant thermophysical properties, as well as non-Newtonian fluid behaviour in the underlying analysis.
Figure 1. Schematic diagram describing the flow of non-Newtonian fluids with initial concentrations 1 and 0 in the upper and lower domains at the inlet of a pipe. A swirl motion consistent with the Rankine vortex is imposed at the pipe inlet. The coordinate system (
$r - \theta - z$) is attached at the centre of the pipe inlet.
Note that the fluids are considered to have same thermophysical properties and, therefore, one set of governing equations is solved for the flow field. Also, the characteristic length scale of the fluidic configurations, typical for a biomicrofluidic set-up, places the underlying flow in the fully developed laminar regime (Re ~ 1 to 100) (Matsunaga & Nishino Reference Matsunaga and Nishino2014; Kaushik, Shyam & Mondal Reference Kaushik, Shyam and Mondal2022).
2.2. Momentum transport: description of flow field
The governing equations (Deen Reference Deen2016) describing the steady, incompressible and laminar flow of non-Newtonian fluid(s) through the chosen fluidic configuration, as shown schematically in figure 1, are the mass, momentum and species transport equations in a cylindrical coordinate system. It is worth mentioning here that we do not consider any body force in the present analysis. Below, we write the continuity and momentum transport equations for describing the underlying flow in the vectorial form
\begin{gather}\boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{u} = 0,\end{gather}
\begin{gather}\rho (\boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{\nabla })\boldsymbol{u} ={-} \boldsymbol{\nabla }p + \boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{\tau }.\end{gather}
Note that, in the above (2.1)–(2.2), 
$\rho$ is the fluid density, p is the fluid pressure, the velocity vector 
$\boldsymbol{u}(r,\theta ,z) = {u_r}{\boldsymbol{i}_r} + {u_\theta }{\boldsymbol{i}_\theta } + {u_z}{\boldsymbol{i}_z}$ and 
$\boldsymbol{\nabla } = (\partial /\partial r){\boldsymbol{i}_r} + (\partial /r\partial \theta ){\boldsymbol{i}_\theta } + (\partial /\partial z){\boldsymbol{i}_z}$.
The components of the deviatoric stress tensor 
$(\boldsymbol{\tau })$ in a cylindrical coordinate system for an incompressible power-law fluid are given as (Bird, Stewart & Lightfoot Reference Bird, Stewart and Lightfoot2006; Deen Reference Deen2016)
\begin{equation}\left. {\begin{array}{r@{}} {{\tau_{rr}} = {\mu_e}\left( {2\dfrac{{\partial {u_r}}}{{\partial r}}} \right),\quad {\tau_{\theta \theta }} = {\mu_e}\dfrac{2}{r}\left( {\dfrac{{\partial {u_\theta }}}{{\partial \theta }} + {u_r}} \right),}\\ {{\tau_{zz}} = {\mu_e}\left( {2\dfrac{{\partial {u_z}}}{{\partial z}}} \right),\quad {\tau_{r\theta }} = {\tau_{\theta r}} = {\mu_e}\left( {r\dfrac{\partial }{{\partial r}}\left( {\dfrac{{{u_\theta }}}{r}} \right) + \dfrac{1}{r}\dfrac{{\partial {u_r}}}{{\partial \theta }}} \right),}\\ {{\tau_{\theta z}} = {\tau_{z\theta }} = {\mu_e}\left( {\dfrac{{\partial {u_\theta }}}{{\partial z}} + \dfrac{1}{r}\dfrac{{\partial {u_z}}}{{\partial \theta }}} \right),\quad {\tau_{zr}} = {\tau_{rz}} = {\mu_e}\left( {\dfrac{{\partial {u_z}}}{{\partial r}} + \dfrac{{\partial {u_r}}}{{\partial z}}} \right).} \end{array}} \right\}\end{equation}
Here, 
${\mu _e}$ is the effective viscosity for the power-law fluid and it is given as 
${\mu _e} = {\mu _0}{(\sqrt {0.5(\boldsymbol{\mathsf{D}}:\boldsymbol{\mathsf{D}})} )^{n - 1}}$, where 
${\mu _0}$ is the flow consistency index, n is the power-law index and 
$\boldsymbol{\mathsf{D}}$ stands for the deformation tensor defined as 
$\boldsymbol{\mathsf{D}}: = \boldsymbol{\nabla }\boldsymbol{u} + {(\boldsymbol{\nabla }\boldsymbol{u})^\textrm{T}}$.
To solve (2.1)–(2.2) analytically, we consider that the flow is axisymmetric with negligible body force and becomes fully developed along the axial direction i.e. 
${u_z} = {u_z}(r)$ (Kumar et al. Reference Kumar, Shakya and Kaushik2020). Considering the aforementioned assumptions, the component of the flow velocity in the radial direction is calculated to be constant using the continuity (2.1) and proven to be zero 
$({u_r} = 0)$, consistent with the no-slip condition at the wall. Proceeding further with this condition and using (2.3), the momentum transport equation in the radial direction reduces to 
$\rho (u_\theta ^2/r) = \partial p/\partial r$, yielding the relation between the radial pressure gradient and the tangential velocity component. It shows that the radial variation of pressure simply supplies the force necessary to keep the fluid elements moving through a circular path within the channel.
Consequently, to solve the momentum transport equation (2.2), the physically justified boundary conditions, based on the aforementioned discussion and assumptions, in compact form are given as: 
${u_r} = 0$; 
$\partial {u_z}/\partial r{|_{r = 0}} = 0$, 
${u_z}(r){|_{r = R}} = 0$; 
${u_\theta }(r,z){|_{r = 0}} = 0$, 
${u_\theta }(r,z){|_{r = R}} = 0$. At the inlet (z = 0) of the considered fluidic configuration, a Rankine vortex (Yao & Fang Reference Yao and Fang2012; Kumar et al. Reference Kumar, Shakya and Kaushik2020) is imposed to create the swirl, i.e. 
${u_\theta }(r,z){|_{z = 0}} = \{ ({u_{\theta ,i,max}})(r/{r_t}),\;r \le {r_t}$ and 
$({u_{\theta ,i,max}})({r_t}(R - r)/r(R - {r_t})),\;r \ge {r_t}$. Here, 
${u_{\theta ,i,max}}$ and 
${r_t}$ represent the maximum inlet swirl velocity and transition radius, respectively. The dimensional transition radius, denoted as 
${r_t}$, signifies the point where the transition from a forced to a free vortex occurs at the channel inlet. Now, in order to solve the flow velocity component in the azimuthal direction, first, we solve the axial momentum, which will then be superimposed to find the tangential velocity component.
2.2.1. Analytical solution
We write the reduced form of the momentum transport equation (2.2) in the axial direction as follows:
\begin{equation}0 ={-} \frac{{\partial p}}{{\partial z}} + \left[ {\frac{1}{r}\frac{\partial }{{\partial r}}(r{\tau_{rz}})} \right] ={-} \frac{{\partial p}}{{\partial z}} + \left[ {\frac{1}{r}\frac{\partial }{{\partial r}}\left( {r{\mu_e}\left( {\frac{{\partial {u_z}}}{{\partial r}}} \right)} \right)} \right].\end{equation}
To obtain (2.4), we refer to the components of the deviatoric stress tensor through (2.3) and use axisymmetric 
$(\partial /\partial \theta = 0)$ as well as fully developed flow 
$({u_z} = {u_z}(r))$ assumptions.
In addition, with the order of magnitude analysis, the effective viscosity is derived as 
${\mu _e} = {\mu _0}{|{\partial {u_z}/\partial r} |^{n - 1}}$ (Sarma, Gaikwad & Mondal Reference Sarma, Gaikwad and Mondal2017). Substituting the value of 
${\mu _e}$ and employing the boundary conditions, i.e. 
$\partial {u_z}/\partial r{|_{r = 0}} = 0$ and 
${u_z}(r/R = 1) = 0$, respectively, we solve (2.4) for the expression of axial velocity, which reads as
\begin{align}{u_z} &= \left( {\frac{n}{{n + 1}}} \right){\left( {\frac{{{R^{n + 1}}|\mathrm{\Delta }p|}}{{2{\mu_0}L}}} \right)^{1/n}}\left[ {1 - {{\left( {\frac{r}{R}} \right)}^{(n + 1)/n}}} \right]\nonumber\\ & = \left( {\frac{{3n + 1}}{{n + 1}}} \right){u_{av}}\left[ {1 - {{\left( {\frac{r}{R}} \right)}^{(n + 1)/n}}} \right].\end{align}
In (2.5), 
${u_{av}} = \int {{u_z}\,\textrm{d}A} /\int {\textrm{d}A}=(n/(3n + 1)){(({R^{n + 1}}|\Delta p|)/(2{\mu _0}L))^{1/n}}$ represents the average axial flow velocity.
We write the simplified form of (2.2) using the aforementioned assumptions to obtain the tangential momentum equation as follows:
\begin{align}{u_z}\frac{{\partial {u_\theta }}}{{\partial z}} &= \frac{1}{\rho }\left[ {\frac{1}{{{r^2}}}\frac{\partial }{{\partial r}}({r^2}{\tau_{r\theta }}) + \frac{\partial }{{\partial z}}({\tau_{\theta z}})} \right] \nonumber\\ &= \frac{1}{\rho }\left[ {\frac{1}{{{r^2}}}\frac{\partial }{{\partial r}}\left( {{r^2}{\mu_e}\left( {r\frac{\partial }{{\partial r}}\left( {\frac{{{u_\theta }}}{r}} \right)} \right)} \right) + \frac{\partial }{{\partial z}}\left( {{\mu_e}\left( {\frac{{\partial {u_\theta }}}{{\partial z}}} \right)} \right)} \right].\end{align}
Employing the assumptions discussed before and using the deviatoric stress components as delineated in (2.3), we make an effort to write (2.6) in the reduced form. Note that (2.6) is the reduced form of the tangential momentum equation and is obtained by employing a few assumptions, i.e. axisymmetric flow 
$(\partial /\partial \theta = 0)$ and fully developed flow 
$({u_z} = {u_z}(r))$ together with the consideration of the continuity equation 
$({u_r} = 0)$, as discussed in § 2.2. As is apparent from (2.6), we have discarded the axial diffusion term i.e. 
${\partial ^2}{u_\theta }/\partial {z^2}$. We perform an order of magnitude analysis to justify the omission of the axial diffusion term from (2.6) as follows. In this analysis, we have imposed swirl motion on the constituent fluids at the inlet to obtain enhanced mixing. Now, even in the presence of an imposed swirl motion, we see that the magnitude of the tangential velocity (
${u_\theta }$) is less than the axial velocity (
${u_z}$) even in the regime very close to the pipe inlet (here, 
${u_{\theta ,max}}/{u_{avg}}{|_{Re = 10,100}}\sim 0.2,\; 0.7$; 
${u_{z,max}}/{u_{avg}} = 2$; cf. figure 2a). Certainly, at further downstream locations of the pipe, the tangential velocity will be even smaller because of the swirl decay compared with the fully developed axial flow velocity. Important to mention is that we have considered 
$L = 120R$ in this analysis, signifying that 
$R \ll L$. Thus, a relatively smaller magnitude of tangential velocity together with larger axial length allow us to discard the axial diffusion term, which is 
${\partial ^2}{u_\theta }/\partial {z^2}\sim {u_\theta }/{L^2} \ll 1$. As such, the magnitude of the axial diffusion term is an order less than the diffusion in the radial direction. Consistent with this order of magnitude analysis, we have discarded the axial diffusion term in (2.6). Substituting the effective viscosity 
${\mu _e} = {\mu _0}{|{\partial {u_z}/\partial r} |^{n - 1}}$ (Sarma et al. Reference Sarma, Gaikwad and Mondal2017) in (2.6), we obtain the dimensionless form of (2.6) in the following form:
\begin{equation}{K_0}U{r^2}\frac{{\partial W}}{{\partial Z}} = \left[ {\frac{\partial }{{\partial r}}({r^{(3n - 1)/n}})\left( {\frac{{\partial W}}{{\partial r}} - \frac{W}{r}} \right)} \right].\end{equation}
In (2.7), the term 
${K_0} = Re{( - 1)^{n - 1}}{(n/(3n + 1))^{n - 1}}$, where the Reynolds number, 
$Re = \rho u_{av}^{2 - n}{R^n}/{\mu _0}$. Here, 
${K_0}$ is analysed under the consideration of real values while varying the power-law index, n. The dimensionless variables appearing in (2.7) are 
$U(r) = {u_z}/{u_{av}}$; 
$W(z,r) = {u_\theta }/{u_{av}}$; 
$r\sim {r^\ast } = r/R$; 
$z\sim {z^\ast } = z/R$, 
${r_t}\sim r_t^\ast{=} {r_t}/R$; where 
$r_t^\ast $ represents the non-dimensional transition radius between a forced and a free vortex. It is worth mentioning here that the chosen scales (velocity and length) in this endeavour are consistent with those used in the seminal works reported in the literature (Kreith & Sonju Reference Kreith and Sonju1965; Ingham et al. Reference Ingham, Keen, Heggs and Morton1990; Tumin Reference Tumin1996; Rouquier, Pothérat & Pringle Reference Rouquier, Pothérat and Pringle2019; Kim Reference Kim2024).
Figure 2. Validation of present analytical swirl velocity profile (a) with Yao & Fang (Reference Yao and Fang2012) at Re = 10 and 100 for Newtonian fluid 
$(n = 1.0)$ and (b) with the three-dimensional numerical simulation (ANSYS) results for non-Newtonian fluid at power-law index, 
$n = 0.8,\;1.0,\;1.2$. The other parameters considered for validation are Re = 100, axial location, z = 1 and transition radius, 
${r_t} = 0.9$. Panel (c) represents the validation of existing experimental results of the axial velocity profile and present numerical model at Re = 26 with limiting case for a Newtonian fluid, n = 1.
To solve (2.7), which is a nonlinear partial differential equation, we apply axisymmetric and no-slip boundary conditions with 
$W(z,r = 0) = 0$ and 
$W(z,r = 1) = 0$, respectively. Note that a Rankine vortex is imposed at the inlet (z = 0) to create the swirl motion therein. The dimensionless form of the Rankine vortex (Yao & Fang Reference Yao and Fang2012; Kumar et al. Reference Kumar, Shakya and Kaushik2020) is given by
\begin{equation}W(z = 0,r) = \left\{ {\frac{{{u_{\theta ,i,max}}}}{{{u_{av}}}}\frac{r}{{{r_t}}},\;r \le {r_t}\;\textrm{and}\;\frac{{{u_{\theta ,i,max}}}}{{{u_{av}}}}\frac{{{r_t}({1 - r} )}}{{r({1 - {r_t}} )}},\;r \ge {r_t}.} \right. \;\end{equation} In this endeavour, we look for the analytical solution of the tangential velocity component. Employing a technique consistent with the separation of variables method and defining 
$W(z,r) = F(z)G(r)$, we solve the nonlinear partial differential equation (2.7) by reducing it to two ordinary differential equations (ODEs) essentially to obtain the tangential velocity component. Below, we write the ODEs with the consideration of 
$\lambda $ as positive constant
\begin{gather}{K_0}\frac{{\textrm{d}F}}{{\textrm{d}z}} + {\lambda ^2}F = 0,\end{gather}
\begin{gather}({r^{(3n - 1)/n}})\frac{{{\textrm{d}^2}G}}{{\textrm{d}{r^2}}} + (2n - 1/n)\frac{{\textrm{d}G}}{{\textrm{d}r}} + ({\lambda ^2}U{r^2} - (2n - 1/n){r^{n - 1}})G = 0.\end{gather} Further, to obtain the analytical series solution, we consider 
$\lambda $ as constant positive real eigenvalue and make use of the axisymmetric boundary condition (i.e. 
$W(z,0) = 0$). We obtain the swirl flow velocity 
$W(z,r)$ in terms of eigenvalues (
${\lambda _m}$; where 
$m = 1,2, \ldots \infty $) and by using symbolic notation for function WhittakerM as
\begin{gather}\textrm{Whit}{\textrm{M}_{m,n}} = \textrm{WhittakerM}\left[ {\frac{{n{\lambda_m}\sqrt {3n + 1} }}{{(2{{(n + 1)}^{3/2}})}},\;\frac{{3n - 1}}{{(2n + 2)}},\;\frac{{(2n{\lambda_m}){r^{((n + 1)/n)}}\sqrt {3n + 1} }}{{{{(n + 1)}^{3/2}}}}} \right].\end{gather} The swirl velocity is obtained by combining the solutions of the aforementioned ODEs (2.9) and (2.10). Below, we present the swirl velocity profile by utilizing the above notation for the Whittaker function (WhittakerM) as ‘
$\textrm{Whit}{\textrm{M}_{m,n}}$’ as follows:
\begin{equation}W(z,r) = \sum\limits_{m = 1}^\infty {{C_m}\,\textrm{exp}( - \lambda _m^2z/{K_0})(\textrm{Whit}{\textrm{M}_{m,n}}/r)} .\end{equation}
Here, 
$\textrm{WhittakerM}(a,b,r)$ is the special function for the solution of Whittaker's equation. Note that 
$\textrm{WhittakerM}(a,b,r)$ is the modified form of the confluent hypergeometric equation and it is defined as (Whittaker Reference Whittaker1903; Abramowitz & Stegun Reference Abramowitz and Stegun1970) 
$\textrm{WhittakerM}(a,b,r) = (\textrm{exp}( - 0.5r))({r^{b + 0.5}})M(b - a + 0.5,\;1 + 2b,\;r)$, where 
$\textrm{M}(a,b,r) = \sum\nolimits_{k = 0}^\infty {({a^k})(k)/({b^k})(k!)}$ is Kummer's confluent hypergeometric function.
Substituting the value of the power-law index 
$n = 1$, representing a special case of a Newtonian fluid, the solution can be written as
\begin{equation}W(z,r) = \sum\limits_{m = 1}^\infty {{C_m}\,\textrm{exp}( - \lambda _m^2z/Re){{\left( {\textrm{WhittakerM}\left[ {\frac{{{\lambda_m}}}{{2\sqrt 2 }},\frac{1}{2},\sqrt 2 {\lambda_m}{r^2}} \right]} \right)} / r}} .\end{equation} We consider the first 30 eigenvalues of (2.12), obtained using the no-slip boundary condition, i.e. 
$W(z,r = 1) = 0$, to get a convergence of the order 10−12 for different values of the power-law index (n = 0.8, 1.0, 1.2). Employing this specified boundary condition, i.e. 
$W(z,r = 1) = 0$
\begin{equation}\sum\limits_{m = 1}^\infty {\textrm{WhittakerM}\left[ {\frac{{n{\lambda_m}\sqrt {3n + 1} }}{{(2{{(n + 1)}^{3/2}})}},\;\frac{{3n - 1}}{{(2n + 2)}},\frac{{(2n{\lambda_m})\sqrt {3n + 1} }}{{{{(n + 1)}^{3/2}}}}} \right]} = 0.\end{equation}We derived the above expression to obtain the eigenvalues. We consider a numerical method, consistent with the iterative approach, to determine the eigenvalues of the function
\begin{equation}f = \sum\limits_{m = 1}^\infty {\textrm{WhittakerM}\left[ {\frac{{n{\lambda_m}\sqrt {3n + 1} }}{{(2{{(n + 1)}^{3/2}}}},\;\frac{{3n - 1}}{{(2n + 2)}},\;\frac{{(2n{\lambda_m})\sqrt {3n + 1} }}{{{{(n + 1)}^{3/2}}}}} \right]} = 0.\end{equation}
The method is applied for different values of n, assuming that f is a continuous function. The characteristic function f is considered to possess m real roots, as denoted by λ 1, λ 2, … and λm. The roots of the real-valued function 
$f = 0$ are determined using the bisection method (Henderson & Aldridge Reference Henderson and Aldridge1992; Grassia Reference Grassia2020). In this method, the interval limits [a, b] for first iteration are initially set to the lower bound a = 1 and upper bound b = 4. The subsequent finite eigenvalues (λi, where i = 1,2,…30) are calculated by ensuring that the condition 
$(f(a) \times f(b)) < 0$ is satisfied, indicating that the roots are located within the specified interval. Here, the bisection method is chosen for its simplicity, reliability and effectiveness in finding isolated single roots in a fixed range. In table 1, we present only the first ten (10) eigenvalues for different values of the power-law index (n = 0.8, 1.2), including n = 1.0 (Newtonian fluid). It may be added here that we undertake an effort to cross-verify our solution methodology, and in doing so, we compare the calculated eigenvalues for n = 1 with the reported results of Yao & Fang (Reference Yao and Fang2012), included in table 1 as well. Important to mention is that the calculated eigenvalues for n = 1.0 conform to those reported by Yao & Fang (Reference Yao and Fang2012). We, however, will make all 30 eigenvalues available to the readers upon request. By appealing to Sturm–Liouville theorem (Kaplan Reference Kaplan1981), we calculate the value of 
${C_m}$, a coefficient appearing in (2.12), by applying the orthogonality condition of the eigenfunctions. We establish orthogonality conditions by considering distinct eigenvalues, denoted as 
${\lambda _m}$ and 
${\lambda _{m1}}$. Specifically, when 
${\lambda _m} \ne {\lambda _{m1}}$, we rely on the Sturm–Liouville theorem (Kaplan Reference Kaplan1981), employing a weight function defined as 
$r(1 - {r^{((n + 1)/n)}})$, within the interval [0,1]. This weight function ensures orthogonality of the eigenfunctions. The orthogonality condition is expressed as follows:
\begin{align}& \int_{r = 0}^{r = 1} {\dfrac{{\textrm{WhittakerM}\left[ {\dfrac{{n{\lambda_m}\sqrt {3n + 1} }}{{(2{{(n + 1)}^{3/2}})}},\dfrac{{3n - 1}}{{(2n + 2)}},\dfrac{{(2n{\lambda_m}){r^{((n + 1)/n)}}\sqrt {3n + 1} }}{{{{(n + 1)}^{3/2}}}}} \right]}}{r}} \nonumber\\ & \quad \times \dfrac{{\textrm{WhittakerM}\left[ {\dfrac{{n{\lambda_{m1}}\sqrt {3n + 1} }}{{(2{{(n + 1)}^{3/2}})}},\dfrac{{3n - 1}}{{(2n + 2)}},\dfrac{{(2n{\lambda_{m1}}){r^{((n + 1)/n)}}\sqrt {3n + 1} }}{{{{(n + 1)}^{3/2}}}}} \right]}}{r}\nonumber\\ &\quad \times r(1 - {r^{((n + 1)/n)}})\,\textrm{d}r = 0. \end{align}Table 1. The first ten eigenvalues for the generalized Whittaker function having values of power-law index, n = 0.8, 1.0 and 1.2.
Utilizing the Sturm–Liouville theorem and the orthogonality of eigenfunctions at 
${\lambda _m} = {\lambda _{m1}}$, with the inlet condition (2.8), the coefficient 
${C_m}$ mentioned in (2.12) can be obtained as
\begin{align}{C_m} &= \int_{r = 0}^{r = 1} {W(0,r)\left( {\frac{{(\textrm{Whit}{\textrm{M}_{m,n}})}}{r}} \right)r(1 - {r^{((n + 1)/n)}})\,\textrm{d}r} /\nonumber\\ &\quad {\left( {\int_{r = 0}^{r = 1} {{{\left[ {\frac{{(\textrm{Whit}{\textrm{M}_{m,n}})}}{r}} \right]}^2}r(1 - {r^{((n + 1)/n)}})\,\textrm{d}r} } \right)} .\end{align} In order to quantify the intensity of swirl along the channel, we define below in (2.18) a dimensionless swirl number 
$S(z)$ as the ratio of the axial flux of angular momentum to the axial flux of axial momentum at a particular cross-section (Reader-Harris Reference Reader-Harris1994; Alekseenko et al. Reference Alekseenko, Kuibin, Okulov and Shtork1999; Maddahian et al. Reference Maddahian, Kebriaee, Farhanieh and Firoozabadi2011)
\begin{align}S(z )= \int_0^R {{{u_z^\ast u_\theta ^\ast {r^{{\ast} 2}}\,\textrm{d}{r^\ast }} / {\left( {R\int_0^R {u_z^{{\ast} 2}{r^\ast }\,\textrm{d}{r^\ast }} } \right)}}} = \int_{r = 0}^{r = 1} {{{UW{r^2}\,\textrm{d}r} / {\left( {\int_{r = 0}^{r = 1} {{U^2}r\,\textrm{d}r} } \right)}}} .\end{align} Using (2.8) and (2.12) in (2.18), the ratio 
$S(z)/S(0)$, i.e. the ratio of the swirl number 
$S(z)$ at any axial location to the that at the inlet, can be obtained as
\begin{align}\frac{{S(z)}}{{S(0)}} = \frac{{\sum\limits_{m = 1}^\infty {{C_m}\,\textrm{exp}( - \lambda _m^2z/{K_0})} \int_{r = 0}^{r = 1} {(\textrm{Whit}{\textrm{M}_{m,n}})r(1 - {r^{((n + 1)/n)}})\,\textrm{d}r} }}{{\int_{r = 0}^{r = 1} {W(0,r){r^2}(1 - {r^{((n + 1)/n)}})\,\textrm{d}r} }}.\end{align}Having established the analytical expression of both the axial and tangential velocity profiles, we quantify mixing of two fluids/solutes as elaborated in the next sub-section. For this part, we first obtain the spatial distribution of the species concentration in the chosen fluidic pathway and then calculate the underlying mixing efficiency for different values of the power-law index.
2.2.2. Benchmarking of analytical method and selection of parameters
To substantiate the efficacy of the proposed theoretical framework, we consider the triple benchmarking strategy as graphically presented in figure 2(a–c). In figure 2(a), we plot the swirl velocity profile obtained from the present theoretical framework in the limiting case of a Newtonian fluid, i.e. n = 1, for two different values of Re (= 10, 100) and compare the same with the reported profile of Yao & Fang (Reference Yao and Fang2012). The other parameters considered for this plot are z = 1, 
${r_t} = 0.9$. It is noteworthy that, to obtain the analytical solution of the swirl velocity using (2.10), the values of λm, crucial in satisfying the no-slip boundary condition i.e. 
$W(z,r = 1) = 0$, are independent of the Reynolds number – a point emphasized by Kumar et al. (Reference Kumar, Shakya and Kaushik2020) and Yao & Fang (Reference Yao and Fang2012). It may be mentioned here that analytical solutions for both no-slip and slip cases have been established in the literature as well (Kumar et al. Reference Kumar, Shakya and Kaushik2020). To verify the convergence of our results, we performed a validation analysis in figure 2(a) by comparing the eigenvalues obtained from our calculation with those reported by Yao & Fang (Reference Yao and Fang2012) for Re = 10 and 100. This benchmarking endeavour ensures the accuracy and consistency of our results. A closer match between the present and published results, as witnessed in figure 2(a), vouches for the efficacy of the proposed analytical model.
Effort has been taken in figure 2(b) to compare our analytical solutions of the flow velocity for different non-Newtonian fluids, including both shear-thinning (n = 0.8) and shear-thickening (n = 1.2) fluids with the corresponding full-scale simulated results. For this validation, we performed three-dimensional simulations employing the finite volume framework of ANSYS Fluent and considered an identical flow configuration as chosen in this endeavour. For the plots depicted in figure 2(b), the other parameters are z = 1, 
${r_t} = 0.9$ at Re = 100. To ascertain the credibility of the theoretical framework developed in our analysis, or more precisely, to ascertain the correctness of the analytical solutions, obtained considering a few assumptions, we performed another validation to compare the analytically obtained flow velocity for non-Newtonian fluids (shear-thinning and shear-thickening fluids) with the full-scale simulated results. We wish to emphasize that the finite volume method framework of ANSYS Fluent is employed in this study to obtain simulated results and subsequently to validate them with the analytical results as presented in figure 2(b). We used the Semi-Implicit Method for Pressure-Linked Equations (SIMPLE) algorithm of the ANSYS Fluent solver for pressure–velocity coupling and employed a second-order upwind scheme to discretize the momentum transport equations (Patankar Reference Patankar1980). We obtain a closer match between analytical swirl velocity profiles and the corresponding numerical results for three different values of n of 0.8, 1.0 and 1.2, as shown in figure 2(b). The comparison analysis presented in figure 2(b) underscores that our analytical solutions, obtained by considering a few physically justified assumptions, faithfully capture the full-scale (three-dimensional) simulated results. It is imperative to note that the validation of the analytical results through a comparison with full-scale numerical simulations using ANSYS Fluent underlines the credibility of the proposed analytical framework. For the sake of completeness, we mention here that the simulations were performed considering a total 12081477 control volumes and assigning a residual criterion of 10−7 for all the transport variables. As evident from figure 2(b), analytical calculations obtained for all the values of n (= 0.8, 1.0, 1.2) match in a fairly accurate manner with the simulated results. Furthermore, we undertake an effort in figure 2(c) to compare the simulated axial velocity profile with the experimental results available in this paradigm (Stokes et al. Reference Stokes, Graham, Lawson and Boger2001b). It is worth mentioning here that the experimental validation is limited to a Newtonian fluid (n = 1.0), while the parameters considered are as follows: Re = 26 and 
${r_t} = 0.5$. Consistent with the experimental configuration as considered by the authors in their study (Stokes et al. Reference Stokes, Graham, Lawson and Boger2001b), we consider only the Rankine vortex with an angular velocity of 13.04 rad s−1 at the inlet. Important to mention is that the swirl specified at the inlet has two different configurations over the domain as follows: force vortex structure up to transition radius and free vortex regime afterwards. Depicted variations in figure 2(c) vouch for a closer as well as consistent match between the experimental and full-scale simulated results without considering the axisymmetric and fully developed flow.
We observe that, for the considered values of Re, our analytical framework developed to solve for the swirl momentum transport can predict the swirl velocity profile reported by Yao & Fang (Reference Yao and Fang2012). Quite notably, the present analytical solutions closely match with full-scale three-dimensional simulated results as well. Notably, this model benchmarking endeavour justifies the accurateness of our modelling framework, which in turn, ascertains the credibility of the analytical solutions as well. As endorsed by the benchmarking analyses, as presented above, we consider the analytically derived velocity profile to obtain the concentration field in the fluidic pathway essentially to save computational time and effort.
For the present analysis, we consider the ranges of Re ~ 1–102 and n ~ 0.6–1.4 (Matsunaga & Nishino Reference Matsunaga and Nishino2014; Cortes-Quiroz, Azarbadegan Zangeneh Reference Cortes-Quiroz, Azarbadegan and Zangeneh2017; Majhi, Nayak & Weigand Reference Majhi, Nayak and Weigand2023). It may be mentioned here that, to analyse the mixing performance in this endeavour, we vary the transition radius ratio from 0.7 to 0.9 by keeping the Péclet number, 
$Pe( = 2600)$ constant (Rezk et al. Reference Rezk, Qi, Friend, Li and Yeo2012).
2.2.3. Description of flow field
Researchers have concluded that the generation of a vortex leads to the chaotic nature between fluid layers and helps to attain efficient mixing of two fluids/solutes in a shorter length. To understand the vortex-assisted mixing of non-Newtonian fluids in a narrow fluidic channel, which is the prime focus of this endeavour, we develop a mathematical model. Our modelling framework employs a semi-analytical technique for solving the momentum transport equation, coupled with the finite volume-based numerical method for the species transport equation, as discussed in § 2.2. In the present analysis, we focus on the rate of change of Rankine vortex strength and its influence on the underlying mixing of two constituent non-Newtonian fluids, investigated for a window of pertinent parameters such as Reynolds number, power-law index and transition radius. We consider an inlet swirl number equal to one throughout this investigation in order to avoid disparaging statements regarding the supremacy of angular momentum and axial momentum in the field. We begin with the description of the swirl velocity profile as demonstrated in figure 2(a,b). As seen in figure 2(a,b), the swirl velocity profile starts to decrease in the outer region along the outward radial direction after it reaches its maximum value.
Figure 3 shows the variation of the swirl velocity and swirl decay profile with a change in the power-law index. For fluids with higher n values, as shown in figure 3, more time is needed for the viscous force to impede the swirl velocity. The magnitude of peak swirl velocity decreases more strongly for fluids with increasing n due to the stronger viscous effect (caused by larger effective viscosity), as shown in figure 3(a). Also, the variation of the swirl velocity profile for different values of transition radius 
$({r_t} = 0.7,\;0.9)$, obtained at Re = 100, z = 1, is shown in figure 3(a). For the shear-thickening fluids with larger n value, a relatively higher effective viscosity promotes the transmission of the zero momentum at the wall deeper into the pipe, as witnessed by the offset of the point of the peak velocity radially inward towards the axis of the channel. This phenomenon further lowers the magnitude of the swirl velocity peak, as confirmed in figure 3(a). Figure 3(a) illustrates how the swirl velocity profile deviates significantly for larger values of 
${r_t}$ with the change in power-law index close to the inlet of the channel (z = 1). This is directly caused by the fact that the higher velocity gradients lead to an increase in the effective viscosity of a shear-thickening fluid. Additionally, given the Rankine vortex is imposed at the inlet, which is evident from the presence of a radial pressure gradient from core to annuls, the transition occurs at 
$r = {r_t}$ is a function of swirl velocity and radius of the channel. The formation of boundary layer ensures that there is a velocity gradient in the flow field. It may be added here that the shear effect predominates closer to the wall and gradually diminishes along the radial direction towards the axis of the channel, which affects the free vortex flow. Although there is a smooth transition from a forced to a free vortex, it is established that, closer to the inlet, the peak value always occurs at radii less than the transition radius. This is due to the frictional effects of the wall on the free vortex region. With decreasing transition radius, the magnitude of the swirl velocity increases in order to retain the momentum in the radial direction for a fixed value of Re.
Figure 3. Plot showing the analytical swirl velocity profiles for different values of power-law index 
$n = 0.8,\;1.0,\;1.2$ and at an axial location z = 1. (a) Variation with pipe radius for two different values of transition radii (
${r_t}$) of 0.7 and 0.9. (b)Variation of maximum swirl velocity with Reynolds number (1 to 100) with transition radius 
${r_t} = 0.7$. (c) Axial variation of maximum swirl velocity for Re = 100, 
${r_t} = 0.7$. (d) Plots depicting the axial variation of swirl intensity for n = 0.8 and 1.2, considering other parameters as 
${r_t} = 0.7$, 0.9 and Re = 100. (e) Qualitative prediction of swirl velocity decay, as shown by the path lines obtained from numerical (ANSYS Fluent) solutions at axial planes (z = 20 and 35) for n = 0.8 and 1.2, while the other parameters considered are Re = 100, 
${r_t} = 0.7$ and 0.9.
Figure 3(b) plots the variation of swirl velocity versus Re in the range 
$1 \le Re \le 100$. As seen from figure 3(b), the magnitude of the swirl velocity increases substantially on increasing the value of Re, while for a given Re, the swirl velocity becomes higher for shear-thinning fluids. The primary flow parameter that plays crucial role in swirl-driven flows is the Re, typically calculated based on the average axial velocity. The higher Re accelerates the swirl's motion into the channel. Consequently, as shown in figure 3(b), the maximum magnitude of the swirl velocity, represented by 
${W_{max}}$, increases with an increase in Re from 1 to 100 for 
${r_t} = 0.7$, calculated at an axial location z = 1. The highest magnitude of swirl tends to remain in its radial location following the interaction between swirl modulated advective effects and fluid viscous effects. As a result, the influence of the channel surface on the wall viscosity causes the swirl velocity to travel as closely as possible along the pipe axis. Also, a more viscous nature of the constituent fluids, as realized by an increasing the shear-thickening effect, results in a reduction of 
${W_{max}}$ with an increase in Re.
We show, in figure 3(c), the decay of the maximum swirl velocity 
$({W_{max}})$ along the axial direction for three different values of n (= 0.8, 1.0 and 1.2), the other parameters being Re = 100 and 
${r_t} = 0.7$. The decay of 
${W_{max}}$ is slower for the shear-thinning fluids compared with the shear-thickening fluids. It is important to observe from figure 3(c) that, for the shear-thinning fluids, the swirl velocity travels much deeper into the channel inside before dissipation. This happens due to the larger velocity gradient near the channel wall that tends to increase the shear stress therein for shear-thickening fluids, which in turn, accelerates the swirl velocity decay compared with shear-thinning fluids.
Figure 3(d) demonstrates the intensity of swirl decay along the axial direction for both the shear-thinning (n = 0.8) and shear-thickening, (n = 1.2) fluids. The plots are depicted for two different values of transition radii 
${r_t} = 0.7$, 0.9 and obtained at Re = 100. It is apparent from the description made in the preceding section(s) that, for a given value of Re, the magnitude of the swirl becomes higher and at smaller transition radius. This finding is consistent with the observation we established in our earlier discussion pertaining to the relationship between the swirl velocity and Reynolds number. This can be explained by the fact that, as the transition radius increases, the swirl velocity gradient increases close to the channel wall. A larger velocity gradient near the channel wall tends to increase the shear stress therein and speed up the swirl velocity's decay. Furthermore, irrespective of the magnitude of Re, the increased viscous effect at higher transition radii tends to speed up the swirl decay for higher values of the power-law index. As can be seen from figure 3(d), the effect of the transition radius is more prominent on the swirl momentum transport of shear-thinning fluids 
$(n < 1)$ compared with shear-thickening fluids 
$(n > 1)$. Also, figure 3(d) witnesses that the swirl velocity decay is somewhat higher for larger transition radii. We attribute this observation to the increased shear stress developed at higher transition radii owing to the larger velocity gradients. It can be inferred from the foregoing discussion as follows: fluids that thin out under shear would have an enhanced transport capability at smaller transition radii than fluids that thicken under shear. Shear-thickening (dilatant) fluids should not be used in applications that demand an extended swirl effect since this special class of fluids promotes swirl to decay more quickly than shear-thinning (pseudoplastic) fluids.
In figure 3(e), we show a qualitative description of swirl transport and its decay for different fluids to facilitate a better understanding of the results at hand. The other parameters considered for this plot are Re = 100, 
${r_t} = 0.7$ and 0.9. The streamlines emerging from the channel entry are portrayed in figure 3(e) to provide greater insights into the swirl decay. Note that the plots depicted are obtained from simulations. Figure 3(e) witnesses (see supplementary movies 1–4 available at https://doi.org/10.1017/jfm.2024.792 for a clear insight) that swirl decay occurs nearly at an axial location z = 35 and 20 for shear-thinning fluids (n = 0.8) and shear-thickening fluids (n = 1.2), respectively, from the inlet of the chosen fluidic configuration. We would like to discuss two important aspects, as observed from figure 3(e). First, for the shear-thinning fluids 
$(n = 0.8)$, the swirl momentum penetrates a greater distance along the channel length for a given strength of inlet swirl, Re (= 100) and 
${r_t}$. This observation on complying with our earlier explanation signifies the impact of rheology modulated fluid viscosity on the swirl transport. Second, the inlet swirl induces an azimuthal motion within the fluid layers, enhancing the contact area and chaotic interaction between the participating fluids. This enhancement in convective-based mixing is particularly significant irrespective of the fluid rheology, emphasizing the enduring impact of the inlet swirl in promoting efficient mixing within the system. This qualitative representation of swirl transport in figure 3(e), derived from full-scale simulations, harmonizes closely with the quantitative analysis of swirl decay presented in figure 3(d). The coherence observed between our analytical solutions and full-scale simulated results once more underlines the credibility of the analytical technique proposed in this analysis.
2.3 Species transport: description of concentration field
Consistent with the assumptions considered in this analysis, except for the axisymmetric one (as the constituent two fluids occupy either half of the pipe inlet cross-section), the simplified species transport equation for zero radial flow velocity can be written in the form as given below
\begin{equation}{u_\theta }\frac{1}{r}\frac{{\partial {C^\ast }}}{{\partial \theta }} + {u_z}\frac{{\partial {C^\ast }}}{{\partial z}} = {D_0}\left[ {\frac{1}{r}\frac{\partial }{{\partial r}}\left( {r\frac{{\partial {C^\ast }}}{{\partial r}}} \right) + \frac{1}{{{r^2}}}\frac{{{\partial^2}{C^\ast }}}{{\partial {\theta^2}}} + \frac{{{\partial^2}{C^\ast }}}{{\partial {z^2}}}} \right].\end{equation}
Here, 
${C^\ast }$ is solute concentration and 
${D_0}$ is the diffusion coefficient of non-Newtonian fluids. To distinguish the species concentration of constituent fluids/solutes in the flow domain, we assign 
${C^\ast } = 1$ for stream A (red colour) and 
${C^\ast } = 0$ for stream B (blue colour), as shown in figure 1.
To non-dimensionalize equation (2.20), we use the velocity and length scales already defined in § 2.2.1, while we define dimensionless concentration 
$C = {C^\ast }/{C_0}$. Here, 
${C_0}$ is the initial concentration of stream A. It is worth mentioning here that, considering a few factors, typically used by the researchers in this paradigm (Matsunaga & Nishino Reference Matsunaga and Nishino2014; Cortes-Quiroz et al. Reference Cortes-Quiroz, Azarbadegan and Zangeneh2017; Majhi et al. Reference Majhi, Nayak and Weigand2023), we select the aforementioned initial species concentration of the candidate fluids. For the sake of the completeness of our ongoing discussion, we here mention a few of those factors as follows: maintaining the desired level of homogeneity or stratification in a laminar flow regime, controlling swirl intensity and direction for fluids with comparable effective viscosities and enabling a sufficiently reduced time to achieve the desired level of convective-based efficient mixing (Matsunaga & Nishino Reference Matsunaga and Nishino2014; Cortes-Quiroz et al. Reference Cortes-Quiroz, Azarbadegan and Zangeneh2017; Majhi et al. Reference Majhi, Nayak and Weigand2023). The non-dimensional form of species transport equation is given as
\begin{equation}\frac{W}{r}\frac{{\partial C}}{{\partial \theta }} + U\frac{{\partial C}}{{\partial z}} = \frac{1}{{Pe}}\left[ {\frac{1}{r}\frac{\partial }{{\partial r}}\left( {r\frac{{\partial C}}{{\partial r}}} \right) + \frac{1}{{{r^2}}}\frac{{{\partial^2}C}}{{\partial {\theta^2}}} + \frac{{{\partial^2}C}}{{\partial {z^2}}}} \right].\end{equation}
In (2.21), the Péclet number 
$Pe( = {u_{av}}R/{D_0})$ is defined as the ratio of the convection strength to the diffusion strength.
To obtain the mixing performance of the constituent fluids in the chosen fluidic device under the modulation of vortical flows, we solve (2.21) using our in-house developed finite volume-based code (Patankar Reference Patankar1980). The analytically derived flow velocity profiles (both tangential and axial velocities) are used to solve the species transport (2.21). The non-dimensional boundary conditions for the species transport (2.21) are outlined next in (2.22a–c)
\begin{gather}\textrm{Along}\;\textrm{the}\;\textrm{radial}\;\textrm{direction:}\quad {\left. {\frac{{\partial C}}{{\partial r}}} \right|_{r = 0}} = 0\textrm{;}\quad {\left. {\frac{{\partial C}}{{\partial r}}} \right|_{r = 1}} = 0,\end{gather}
\begin{gather}\textrm{Along}\;\textrm{the}\;\textrm{azimuthal}\;\textrm{direction:}\quad {C_{(\theta = 0)}} = {C_{(\theta = 2{\rm \pi} )}};\quad {\left. {\frac{{\partial C}}{{\partial \theta }}} \right|_{\theta = 0}} = {\left. {\frac{{\partial C}}{{\partial \theta }}} \right|_{\theta = 2{\rm \pi} }},\end{gather}
\begin{gather}
\textrm{Along}\;\textrm{the}\;\textrm{axial}\;\textrm{direction:}\quad
{C_{(z\; = \; 0)}} = \left\{ {\begin{array}{*{20}{l}}
{1,}&{0 \leq \theta < {\rm \pi}}\\ {0,}&{{\rm \pi} \leq \theta <
2{\rm \pi} } \end{array}} \right.;\quad {\left. {\frac{{\partial
C}}{{\partial z}}} \right|_{z = L}} =
0.\end{gather}
We briefly discuss the numerical framework employed in this study to solve (2.21) using the aforementioned boundary conditions (2.22a–c) in the forthcoming section.
2.3.1. Numerical analysis
As mentioned before, we solve the species transport (2.21) numerically by employing our in-house developed finite volume code. We use a power-law scheme for the discretization of (2.21) to obtain converged solutions for a higher value of the Péclet number (Pe) using the boundary conditions given in (2.22a–c).
The mixing performance is typically assessed by analysing the variation of mixing efficiency with respect to the mixing methods (active or passive) and by varying the involved parameters such as Péclet number, characteristic length scale of the device etc., within their permissible range. In the present endeavour, to elucidate the mixing performance, we calculate mixing efficiency 
${\eta _m}$ by using the expression given below (Wang et al. Reference Wang, Wang, Feng and Lin2015; Gaikwad, Kumar & Mondal Reference Gaikwad, Kumar and Mondal2020; Shyam et al. Reference Shyam, Mondal and Mehta2021; Kaushik et al. Reference Kaushik, Shyam and Mondal2022)
\begin{gather}{\eta _m}(z) = \left[ {1 - \left( {{{\int_A {|{{C_i}(r,\theta ,z) - {C_\infty }} |\,\textrm{d}r\,\textrm{d}\theta } } / {\left( {\int_A {|{{C_0}(r,\theta ) - {C_\infty }} |\,\textrm{d}r\,\textrm{d}\theta } } \right)}}} \right)} \right] \times 100\,\%.\end{gather}
In (2.23), 
${C_\infty }$ represents the dimensionless concentration of species at a perfectly mixed state, which is usually taken as 0.5.
As evident from (2.23), 
${\eta _m}$ gives the efficiency at each axial location of the channel. Therefore, we can estimate the mixing length from this expression, which, in turn, will allow us to predict the required size of the proposed mixing assay. We mention here that, at the channel inlet, 
${C_0} = 1$, 0 in the upper half and lower half, respectively, as shown in figure 1. Note that 
${\eta _m} \in (0,1)$ in (2.23) can take any value between zero and one. Important to mention here is that the zero value implies no mixing of the constituent components, while a value of unity indicates a completely mixed state of the participating components.
2.3.2. Grid performance analysis
Just to ascertain that the mixing efficiency, estimated numerically in the present analysis, does not have any grid resolution bias, we perform a grid performance test for a set of other parameters as n = 0.8 at Pe = 2600, 
${r_t} = 0.9$ and Re = 100, as shown in figure 4. In doing so, we calculate the grid convergence index (GCI), as shown in figure 4(a).
Figure 4. (a) Plot showing the GCI for three different grid refinements, defined with a dummy variable for 
$\varDelta = {1.2^0},\;{1.2^1},\;{1.2^2}$. (b) The mixing efficiency at the pipe outlet is plotted for three distinct grid refinements in all directions as considered for GCI analysis. The other parameters considered for these plots (a,b) are Reynolds number (Re = 100), transition radius (
${r_t} = 0.9$), Péclet number (Pe = 2600) and power-law index (n = 0.8). (c) Typical grid structure and distribution are shown for the fluidic configuration considered here with an axial distance of 120 times the radius (z = 120R).
To conduct this analysis, we considered a grid number with a dummy variable Δ, proportional to the dimension of the chosen model, in three mutually perpendicular directions. Here, we vary Δ with a spacing ratio of 1.2 in three mutually perpendicular directions having finer grid distributions for the 
$\varDelta = {1.2^1}$ and 
$\varDelta = {1.2^2}$ cases than that of the case 
$\varDelta = {1.2^0}$. Note that a typical grid structure corresponding to 
$\varDelta = {1.2^0}$ pertaining to the chosen fluidic configuration with an axial distance of 120 times the pipe radius (L = 120R) is shown in figure 4(c). From the depicted results of GCI for two sets of grid refinements in figure 4(a), we find that 
$\textrm{GC}{\textrm{I}_{23}}/{1.2^P}\textrm{GC}{\textrm{I}_{12}} \approx 1$ with an order of convergence 
$P = 1.129$, where 
$\textrm{GC}{\textrm{I}_{12}} = 4.98\,\%$ and 
$\textrm{GC}{\textrm{I}_{23}} = 6.12\,\%$. This calculation of GCI suggests that the solutions indeed lie comfortably within the asymptotic range of convergence, affirming the reliability of the results of this endeavour.
Following this GCI plot, we show in figure 4(b) the axial variation of mixing efficiency for three distinct grid refinements in all directions, as considered for GCI analysis. The parameters considered for this plot are n = 0.8, Pe = 2600, 
${r_t} = 0.9$ and Re = 100. Remarkably, with a change in grid refinement from 
$62.40 \times {10^6}$ to 
$107.83 \times {10^6}$ elements, we observe an insignificant variation in the mixing efficiency, with a percentage change below 1% at a specific axial location, z = 100. Based on these observations from figure 4(a,b), we consider 
$\varDelta = {1.2^1}$, equivalently 100 × 260 × 2400 (= 62.4 × 106) grids for the subsequent analysis, as discussed in the upcoming sections.
2.3.3. Solute mixing: prediction, transition and efficiency
From the discussion above in § 2.2.3, it is apparent that the shear-thinning fluids have a significant influence on the swirl velocity. Also, we understand from the foregoing discussion that, for the shear-thinning fluids (n < 1), convection has a dominant influence over molecular diffusion on the underlying mixing for higher values of Re, as discussed next.
The variation of mixing efficiency along the axial direction for different values of power-law index is depicted in figure 5. The other parameters considered for this plot are Pe = 2600, Re = 100 and 
${r_t} = 0.7$. We can see from figure 5(a) that, with a decreasing value of the power-law index, the mixing efficiency increases along the axial direction. It may be mentioned here that, for the shear-thinning fluids (n < 1), the effect of convection on the underling mixing is better realized for a decreasing value of the transition radius and increasing magnitude of Re. The higher strength of rotational inertia for Re = 100 and more shear-thinning behaviour of the participating fluids for n = 0.8 lead to an increase of the tangential flow velocity at a given strength of inlet swirl. The higher tangential velocity is responsible for the augmented chaotic nature of the candidate fluids due to an increase in the contact area between them, which, in turn, promotes convective mixing along the axial direction. It is because of this reason that the mixing efficiency of shear-thinning fluids becomes higher, as witnessed in figure 5(a). We show, in figure 5(b), the influence of transition radii on the axial variation of mixing efficiency. The inference obtained from the depicted variation of the swirl velocity profile versus transition radius (see figure 3a) suggests that the swirl intensity is higher for a smaller transition radius, i.e. 
${r_t} = 0.7$. As can be seen from figure 5(b), for the smaller value of 
${r_t}$, the mixing efficiency at the channel outlet becomes higher, attributed primarily to a higher intensity of swirl. Notably, for a given 
${r_t}$, the mixing efficiency at any axial location becomes higher for the shear-thinning fluids (n < 1) compared with the shear-thickening fluids (n > 1), as seen in figure 5(b). For a specified strength of inlet swirl, the rate of swirl decay is slower for shear-thinning fluids (see figure 3d), which renders better mixing efficiency for this class of fluids at any axial location.
Figure 5. The plot of mixing efficiency along the axial direction: (a) with change in value of the power-law index (n = 0.8, 1.0, 1.2), where the other parameters considered for the plot are Re = 100, Pe = 2600 and 
${r_t} = 0.7$; and (b) at two different transition radii (
${r_t} = 0.7$, 0.9) for shear-thinning (n = 0.8) and shear-thickening fluids (n = 1.2), where the other parameters considered for the plot are Re = 100, Pe = 2600.
As the prime focus of this endeavour is to investigate the mixing performance of non-Newtonian fluids/solutes, we make an attempt in figure 6 to plot the variation of mixing efficiency 
$({\eta _m})$ with flow behaviour index (n), characterizing the rheological behaviour of the inelastic non-Newtonian fluids considered in this analysis. We may mention here that the mixing transition is a consequence of the appearance and development of small-scale three-dimensional swirl structures in the flow field. The large excursions corresponding to large vortex structures typically depend on Re, and are responsible for convective-based mixing in shorter length (Breidenthal Reference Breidenthal1981). Precisely by depicting figure 6, we would like to demonstrate how the transition in mixing occurs in a swirl-driven flow environment as modulated by the fluid rheology. The other parameters are Re = 100, Pe = 2600 and 
${r_t} = 0.7$.
Figure 6. The plot shows qualitative aspect of mixing efficiency using concentration contour with change in the power-law index (n = 0.6, 0.8, 1.0, 1.2, and 1.4) at two different axial locations, z = 30 and 120. The indicated concentration contours are used to show the qualitative aspect of advective mixing for a shear-thinning fluid at n = 0.6 compared with diffusive mixing for a shear-thickening fluid at n = 1.4. It is based on the complete rotation of fluid at n = 0.6 compared with n = 1.4. The other parameters considered for this analysis are Re = 100, Pe = 2600 and 
${r_t} = 0.7$.
The higher internal convection of the constituent fluids, which follows a nonlinear trend with increasing n, as demonstrated by the concentration contours in figure 6, leads to better mixing being achieved in swirl-driven flow environments. It may be mentioned here that the concentration contours give a good indication of the impact of swirl on vortex formation. At an axial location z = 30, the concentration contour for n = 0.6 (see concentration contour of figure 6) exhibits a complete bulk fluid rotation compared with the reference vortex formed for n = 1.4 (see supplementary movies 5–7 for a clear insight). This observation indicates that convection-based mixing dominates over molecular diffusion for n = 0.6. Depicted contours in figure 6 suggest that, for more shear-thickening behaviour of the candidate fluids (n = 1.4), mixing based on molecular diffusion is clearly visible to the naked eye at an axial location of z = 30 compared with z = 120. As can be verified from figure 6, with increase of the shear-thinning nature of the constituent components, the underlying mixing is seen to become more pronounced regardless of the magnitude of Re, transition radius and Pe. A lesser apparent viscosity of constituent fluids having a more shear-thinning nature for a given deformation rate (here, strength of inlet swirl) inhibits faster attenuation of swirl for its given strength, which, in turn, ensures better mixing performance at any axial location, as witnessed in figure 6.
It may be mentioned here that the transmission of a given strength of the inlet swirl through the participating fluids largely relies on the momentum transport, which, in turn, depends on the Reynolds number. Thus, in a swirl-driven flow environment, underlying mixing of constituent fluids will be governed by both the inevitable molecular diffusion and chaotic advection. Considering this aspect, we plot in figure 7(a) the variation of mixing efficiency with Re for three different values of n (= 0.8, 1.0 and 1.2). The following parameters are considered for this plot as Pe = 2600, z = 100 and 
${r_t} = 0.7$. We discuss two important observations as follows. First, the mixing efficiency is seen to be higher for shear-thinning fluids 
$(n < 1)$ compared with shear-thickening fluids 
$(n > 1)$. We attribute this observation to the lesser apparent viscosity of shear-thinning fluids, which prevents faster swirl decay. Second, a higher inertial effect for higher Reynolds number accelerates the swirl motion penetrating much deeper into the channel before attenuation. However, at lower Re, impeding viscous resistance leads to swirl decay over shorter distances from the inlet. This observation is justifiable qualitatively from the mixing concentration contours depicted in figure 7(b). It is because of the diminishing effect of inlet swirl at low Re that we observe a linear trend of mixing in regime-I of figure 7(a), attributed primarily to the diffusion-based mixing. The contours depicted in figure 7(b) suggest that, for the smaller values of Re (= 1, 10), the effect of inlet swirl is not translated into the downstream locations of the fluidic configuration. This observation is attributed to the weaker rotational fluid motion developed at the pipe inlet for smaller values of Re. It is because of this reason that we obtain diffusion-based mixing for smaller Re even at a further downstream location of the pipe, as witnessed in figure 7. The role of convention-based mixing is shown in regime-II, which follows the nonlinear trend. For Re greater than 10, bulk fluids rotate more than 
$90^\circ$. This flow structure helps to increase the contact surface area between the fluid layers and enhances convection-based mixing, as supported by the contours shown in figure 7(b).
Figure 7. (a) The plot shows mixing efficiency (left side) as a function of Reynolds number (1 to 100), as well as shear thinning and thickening (n = 0.8 and 1.2, respectively), in addition to Newtonian fluid (n = 1.0) at an axial location, z = 100. (b) The qualitative aspect of mixing efficiency (right side) using concentration contours with change in the same n at Re = 1, 10, 40, 80 and 100 is shown. On increasing the Reynolds number, advective mixing (regime II) is obtained at Re = 100 compared with diffusive (regime I) at Re = 1 and 10 based on complete rotation of fluid compared with only twist, respectively. The other parameters considered for this analysis are Pe = 2600 and 
${r_t} = 0.7$.
The contours illustrating the mixing efficiency in figures 6 and 7 establish the dependency of the underlying mixing on the power-law index 
$(n \in 0.6\hbox{--}1.4)$ and Reynolds number 
$(Re \in {10^0}\hbox{--}{10^2})$. The discussion made above in the context of figures 6 and 7 underscores that convection-dominated mixing becomes prominent for shear-thinning fluids with higher Re. To substantiate this observation, we next make an effort to establish a correlation among the mixing efficiency 
$({\eta _m})$, power-law index and Reynolds number. The developed correlation, represented by a second-order polynomial, exhibits a high degree of fit with a coefficient of determination value of 0.9935. The developed explicit, quadratic form of 
${\eta _m}$ in terms of the independent variables, i.e. power-law index and Reynolds number at the pipe outlet (z = 120) can be expressed as follows:
\begin{equation}{\eta _m} = {a_1} + {a_2}n + {a_3}Re - {a_4}{n^2} - {a_5}nRe + {a_6}R{e^2}.\end{equation}
Here, the constants 
${a_1},\;{a_2},\;{a_3},\;{a_4},\;{a_5}$ and 
${a_6}$ are 23.31, 15.92, 0.8756, 10.9, 0.3503 and 0.0002416, respectively. Note that the proposed correlation in (2.24) corresponds to the data set obtained for 
${r_t} = 0.7$ and 
$Pe\; = \; 2600$.
On using (2.24), we obtain the following insights, as discussed next. The mixing efficiency at the channel outlet (z = 120) becomes higher with increasing magnitude of Re and for an enhanced shear-thinning nature of the fluid. This inference effectively supports the findings depicted in figures 6 and 7. Hence, it becomes evident that the convective-based mixing efficiency is significantly influenced by both the rheology of the constituent fluids and the flow Reynolds number.
3. Summary, perspective and outlook
In the present analysis, we investigate the effect of swirl flow (vorticial flow) on the mixing of non-Newtonian fluids in a narrow cylindrical pipe under laminar flow conditions. We analytically derive the expression for the swirl velocity profile by superimposing the fully developed flow and applying the Rankine vortex condition at the channel inlet. We incorporate analytically obtained flow velocities (both swirl (tangential) and axial velocities) into the species transport equation to obtain the concentration distribution of the constituent fluids/solutes along the cross-section of the chosen fluidic configuration. We examine the influence of inlet swirl by varying pertinent parameters such as flow behaviour (power-law) index (
$n = 0.6\hbox{--}1.4)$, Reynolds number 
$(Re = {10^0}\hbox{--}{10^2})$ and transition radius 
$({r_t} = 0.7\hbox{--}0.9)$.
Our findings reveal that an increase in the value of the power-law index, signifying the change in behaviour of the constituent fluids from shear thinning to shear thickening, results in an increase in the fluids’ apparent viscosity for a given inlet swirl. We observe that the decay of swirl is critically dependent on the magnitude of the power-law index or flow behaviour index. We unveil through this analysis that both the power-law index and the Reynolds number significantly impact the length of swirl decay. This increased apparent viscosity of the constituent fluids with higher values of power-law index leads to a reduction in the magnitude of both axial and tangential velocities. It is demonstrated that decreasing the value of the transition radius and increasing the magnitude of the Reynolds number leads to an enhancement of swirl intensity to act over a greater extent of the flow configuration, attributed primarily to the combined effects of a higher value of the radial pressure gradient and a reduced wall shear stress. Based on the findings obtained from our simulations, we elucidate how swirl alters the flow structure, inducing rotation and enhancing the contact surface area of the participating fluids. Our study demonstrates that chaotic convection, or the swirl-driven rotation of the fluid, plays a significant role in enhancing mixing efficiency at a higher Reynolds number. We believe that the insights gained from this study may strengthen the fundamental idea of vortex-assisted mixing for non-Newtonian fluids in a narrow fluidic confinement. The current work may be particularly valuable for researchers aiming to enhance transport capabilities and convective mixing through the use of swirl flow, especially by altering the rheological properties of base fluids, like water through a polymer dilution, to change their rheological behaviour.
Supplementary movies
Supplementary movies are available at https://doi.org/10.1017/jfm.2024.792.
Acknowledgements
D.K. gratefully acknowledges funding from the Prime Minister Research Fellowship (PMRF), Ministry of Education, Govt. of India to carry out this research work. The authors sincerely thank the anonymous reviewers for their insightful comments to improve the quality of the manuscript.
Declaration of interests
The authors report no conflict of interest.
Author contributions
D.K.: formal analysis, data curation, conceptualization, investigation, methodology, software, visualization, validation, writing – original draft. P.K.M.: conceptualization, funding acquisition, methodology, project administration, resources, supervision, validation, writing – review, and editing.
