1. Introduction
A solid sphere translating in a cylindrical tube filled with a quiescent liquid has been one of the classical and benchmark problems at the forefront of transport phenomena for many decades. It represents an idealization of many industrially relevant processes, examples being fluidized and fixed bed reactors, slurry reactors, falling ball viscometers, equipment for separating solid–liquid mixtures in mining and petrochemical industries, processing of polymer melts, etc. Not only of practical significance, this problem is also of fundamental interest in its own right. As a result, this problem has been extensively investigated in the research community, and much has been written about it in the literature for both Newtonian and non-Newtonian fluids (McKinley Reference McKinley2002; Chhabra Reference Chhabra2006; Michaelides Reference Michaelides2006). Earlier investigations of this problem were restricted to simple Newtonian fluids like water, and it was then gradually extended to complex non-Newtonian fluids like polymer solutions and melts due to their overwhelming applications in scores of industrial settings like food, petrochemicals, personal care products, etc. (Chhabra Reference Chhabra2006). Earlier investigations revealed that both the blockage ratio (ratio of sphere to tube diameter) and nonlinear rheological properties of fluids like shear-thinning, shear-thickening, viscoplasticity, etc., greatly influenced the flow characteristics like the drag force, wake length, etc., in comparison to an unconfined situation and for Newtonian fluids. In addition to the investigations carried out for generalized Newtonian fluids, many studies have also been presented concerning viscoelastic fluids. Some typical and complex flow features were seen in these fluids as compared to those seen either in Newtonian fluids or any generalized Newtonian fluid. This complexity was not only observed in the variation of integral parameters like the drag force but also seen in the flow fields near the sphere. For instance, a downward or upward shifting in the axial velocity profile along the upstream or downstream axis of the sphere has been observed both experimentally and numerically for viscoelastic fluids in comparison to that seen in Newtonian fluids (Bush Reference Bush1994; Arigo et al. Reference Arigo, Rajagopalan, Shapley and McKinley1995; Arigo & McKinley Reference Arigo and McKinley1998). Additionally, a flow reversal phenomenon and/or the presence of a ‘negative wake’ downstream of the sphere have also been observed in viscoelastic fluids (Bisgaard Reference Bisgaard1983; Harlen Reference Harlen2002).
The next generation of studies on this benchmark problem has considered a solution comprised of various types of surfactant molecules. When these molecules dissolve in water above a critical concentration they spontaneously self-assemble into large and flexible aggregates of micelles of different shapes like spherical, ellipsoidal, wormlike or lamellar (Moroi Reference Moroi1992). The rheological properties of these wormlike micellar (WLM) solutions were found to be more complex than those seen for polymer solutions or melts (Rothstein Reference Rothstein2003, Reference Rothstein2008). This is due to the fact that these wormlike micelles can undergo continuous scission and reformation in an imposed shear or extensional flow field, unlike polymer molecules that are unlikely to break due to the presence of a strong covalent backbone. Because of extensive applications over a wide range of industrial settings, a considerable number of studies have also been performed on the falling-sphere problem in these fluids in the creeping flow regime. For instance, Jayaraman & Belmonte (Reference Jayaraman and Belmonte2003) conducted an experimental investigation of this problem in cetyltrimethylammonium bromide (CTAB)/sodium salicylate (NaSal) WLM solution. They found an unsteady motion of the sphere in the direction of its sedimentation. They proposed that the cause of this instability was due to the destruction of the flow-induced structure formed in the sphere's vicinity. However, in a later study with the same WLM solution, Chen & Rothstein (Reference Chen and Rothstein2004) claimed that this instability was due to the sudden rupture of the long micelles downstream of the sphere. This reason was further established in a later study with cetylpyridinium chloride (CPyCl)/NaSal WLM solution by Wu & Mohammadigoushki (Reference Wu and Mohammadigoushki2018). The unsteady motion of a falling sphere was also observed in the study by Kumar et al. (Reference Kumar, Majumdar, Sood, Govindarajan, Ramaswamy and Sood2012) with cetyl trimethyl ammonium 
$p$-toluenesulphonate (CTAT)/ sodium chloride micellar solution and in a recent study by Wang et al. (Reference Wang, Wang, Xu, Dou and Su2020) with octadecyl trimethyl ammonium chloride (OTAC)/NaSal micellar solution. To characterize the onset of this unsteady motion of the falling sphere, Mohammadigoushki & Muller (Reference Mohammadigoushki and Muller2016) and Zhang & Muller (Reference Zhang and Muller2018) found a criterion by calculating the local extensional Weissenberg number downstream of the sphere based on the local maximum extension rate. This criterion is found to be universally valid as they discovered that it does not depend on the micelle chemistry and solution rheological behaviours. For instance, it does not depend on whether the solution shows shear-banding phenomena or not. Furthermore, their predictions for the unsteady motion were in line with that predicted by Chen & Rothstein (Reference Chen and Rothstein2004) and Wu & Mohammadigoushki (Reference Wu and Mohammadigoushki2018).
Therefore, most studies have proposed that the unsteady motion of a sedimenting sphere in WLM solutions is due to the presence of strong extensional flow downstream of the sphere, causing the sudden rupture of highly aligned and stretched micelles in this region (Rothstein Reference Rothstein2008). However, there was no direct evidence presented for this, and it was only indirectly proved using flow-induced birefringence and particle image velocimetry experiments (Chen & Rothstein Reference Chen and Rothstein2004; Wu & Mohammadigoushki Reference Wu and Mohammadigoushki2018). The present study aims to establish this hypothesis using numerical simulations based on the Vasquez–Cook–McKinley (VCM) constitutive model (Vasquez, McKinley & Cook Reference Vasquez, McKinley and Cook2007) for a WLM solution. However, it should be mentioned here that the problem considered in this study is not the exact representation of the prior experimental settings wherein the sphere is allowed to sediment in a tube due to its own weight (Jayaraman & Belmonte Reference Jayaraman and Belmonte2003; Chen & Rothstein Reference Chen and Rothstein2004; Wu & Mohammadigoushki Reference Wu and Mohammadigoushki2018; Zhang & Muller Reference Zhang and Muller2018). The sphere may rotate or undergo lateral motion during the sedimentation or even it may not reach a terminal velocity (Mohammadigoushki & Muller Reference Mohammadigoushki and Muller2016). Therefore, in actual experiments, the flow may become three-dimensional and non-axisymmetric. To realize the corresponding experimental conditions accurately, one has to solve numerically the full governing field equations, namely continuity, momentum and micellar constitutive equations in a three-dimensional computational domain along with an equation of the sphere motion. In the present simulations, we consider a problem wherein the sphere is translating steadily along the axis of a tube, and this can be a situation in the corresponding experiments of the falling-sphere problem when the sphere will reach a terminal velocity. Although this is not the case in actual experiments, by using this simplified problem, we aim to show that this unsteady motion downstream of the sphere is, indeed, caused by the breakage of micelles. Therefore, this will further establish the hypothesis for the unsteady motion of a sphere falling in WLM solutions, as observed in prior experiments (Jayaraman & Belmonte Reference Jayaraman and Belmonte2003; Chen & Rothstein Reference Chen and Rothstein2004; Wu & Mohammadigoushki Reference Wu and Mohammadigoushki2018; Zhang & Muller Reference Zhang and Muller2018).
To prove the aforementioned hypothesis, as stated above, the present study plans to use the two-species VCM constitutive model for characterizing the rheological behaviour of WLM solutions. This model considers the micelles as elastic segments composed of Hookean springs, which all together form an elastic network. The breakage and reformation dynamics were incorporated in this model based on Cate's original reversible breaking theory for wormlike micelles (Cates Reference Cates1987). For different viscometric flows, a very good agreement has been found between the predictions obtained with the VCM model and the corresponding experimental results (Pipe et al. Reference Pipe, Kim, Vasquez, Cook and McKinley2010; Mohammadigoushki et al. Reference Mohammadigoushki, Dalili, Zhou and Cook2019) whereas for a non-viscometric complex flow, a good qualitative correspondence has been seen in recent studies (Kalb, Villasmil U. & Cromer Reference Kalb, Villasmil U. and Cromer2017, Reference Kalb, Villasmil U. and Cromer2018; Khan & Sasmal Reference Khan and Sasmal2020; Sasmal Reference Sasmal2020). Therefore, this VCM model's capability of predicting the flow behaviour of WLM solutions in various flow fields is well established. We also use the single-species Giesekus constitutive equation in our analysis to show the importance of breakage and reformation of micelles for the onset of the unsteady motion.
2. Problem formulation and governing equations
The problem considered herein is the study of the flow characteristics of a sphere of diameter 
$d$ translating steadily along the axis of a cylindrical tube of diameter 
$D$ filled with an incompressible WLM solution in the creeping flow regime, as schematically shown in figure 1(a). The present problem is solved in an Eulerian reference frame wherein the coordinate system is centred on and travelling with the sphere. In this coordinate system, the velocity vector is assumed to be zero on the sphere surface. Furthermore, at the inlet and tube walls, the dimensionless fluid axial velocity is set to be unity, and the radial velocity is set to be zero (discussed in detail in the subsequent section) as schematically shown in figure 1(a). Furthermore, the flow is assumed to be two-dimensional and axisymmetric in nature.
Figure 1. (a) Schematic of the present problem with both Cartesian and spherical coordinates. (b) Different mesh densities used in the present study with a zoomed view near the sphere surface. (c) Implementation of the wedge boundary condition to approximate the two-dimensional and axisymmetric condition of the present problem in OpenFOAM.
Two values of the blockage ratio (the ratio of the sphere diameter to the tube diameter, i.e. 
$d/D$), namely 0.33 and 0.1, are considered in this study, and the upstream 
$(L_{u})$ and downstream 
$(L_{d})$ lengths of the tube are chosen as 65
$d$. These values are sufficiently high so that the end effects are negligible. This was further confirmed by performing a systematic domain independence study.
2.1. Flow equations
Under the circumstances mentioned above, the flow field will be governed by the following equations written in their dimensionless forms.
Equation of continuity:
\begin{equation} \boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{U} = 0. \end{equation}Cauchy momentum equation:
\begin{equation} El^{{-}1} \frac{\textrm{D}\boldsymbol{U}}{\textrm{D}t} ={-}\boldsymbol{\nabla} P + \boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{\tau}. \end{equation}
In the above equations, 
$\boldsymbol {U}$, 
$t$ and 
$\boldsymbol {\tau }$ are the velocity vector, time and total extra stress tensor, respectively, and 
$El$ is the elasticity number defined at the end of this section. For an intertialess flow, the left-hand side of (2.2) is essentially zero. The total extra stress tensor, 
$\boldsymbol {\tau }$, for a WLM solution is given as
\begin{equation} \boldsymbol{\tau} = \boldsymbol{\tau_{w}} + \boldsymbol{\tau_{s}}, \end{equation}
where 
$\boldsymbol {\tau _{w}}$ is the non-Newtonian contribution from the wormlike micelles and 
$\boldsymbol {\tau _{s}}$ is the contribution from the Newtonian solvent which is equal to 
$\beta \dot {\boldsymbol {\gamma }}$. Here the parameter 
$\beta$ is the ratio of the solvent viscosity to the zero-shear-rate viscosity of the WLM solution and 
$\dot {\boldsymbol {\gamma }} = \boldsymbol {\nabla } \boldsymbol {U} + \boldsymbol {\nabla } \boldsymbol {U}^{T}$ is the strain-rate tensor. For the two-species VCM model, the total extra stress tensor is given by
\begin{equation} \boldsymbol{\tau} = \boldsymbol{\tau}_{w}^{VCM} + \boldsymbol{\tau_{s}} = (\boldsymbol{A} + 2\boldsymbol{B}) - \left(n_{A} + n_{B}\right)\boldsymbol{I} + \beta_{VCM}\dot{\boldsymbol{\gamma}}. \end{equation}
Here 
$n_{A}$ and 
$\boldsymbol {A}$ are the number density and conformation tensor of the long worm A, respectively, and 
$n_{B}$ and 
$\boldsymbol {B}$ are those of the short worm B in the two-species model. The temporal and spatial evaluations of the number density and conformation tensor of worms are considered in the following subsection. For the single-species Giesekus model, this is given by
\begin{equation} \boldsymbol{\tau} = \boldsymbol{\tau}_{w}^{G} + \boldsymbol{\tau_{s}} = ( \boldsymbol{A} - \boldsymbol{I}) + \beta_{G}\dot{\boldsymbol{\gamma}}. \end{equation}
Note that here all the lengths, velocity, time and conformation tensors are non-dimensionalized using 
$d$, 
$d/\lambda _{eff}$, 
$\lambda _{eff}$ and 
$G_{0}^{-1}$, respectively, where 
$\lambda _{eff} = {\lambda _{A}}/({1+c_{Aeq}^{\prime}\lambda _{A}})$ is the effective relaxation time in the two-species VCM model, 
$G_{0}$ is the elastic modulus and 
$\lambda _{A}$ and 
$c_{Aeq}^{\prime}$ are the dimensional relaxation time and equilibrium breakage rate of the long chain A, respectively. In the case of the single-species model, 
$\lambda _{eff}$ is replaced by the Maxwell relaxation time 
$\lambda$ during the non-dimensionalization. The elasticity number is defined as 
$El = {Wi}/{Re}$, where 
$Wi_{S} = {\lambda _{eff}U}/{d}$ is the shear Weissenberg number and 
$Re = {d U \rho }/{\eta _{0}}$ is the Reynolds number.
2.2. Two-species VCM constitutive equations
The VCM constitutive equations provide the species conservation equations for the long 
$(n_{A})$ and short 
$(n_{B})$ worms along with the equations for the evolution of the conformation tensors of the long 
$(\boldsymbol {A})$ and short 
$(\boldsymbol {B})$ worms. According to this model, the equations for the variation of 
$n_{A}$, 
$n_{B}$, 
$\boldsymbol {A}$ and 
$\boldsymbol {B}$ are given in their non-dimensional forms as follows (Vasquez et al. Reference Vasquez, McKinley and Cook2007):
\begin{gather} \mu\frac{\textrm{D}n_{A}}{\textrm{D}t} - 2\delta_{A} \nabla^{2}n_{A} = \frac{1}{2} c_{B} n_{B}^{2} - c_{A}n_{A}, \end{gather}
\begin{gather} \mu\frac{\textrm{D}n_{B}}{\textrm{D}t} - 2\delta_{B} \nabla^{2}n_{B} ={-} c_{B} n_{B}^{2} + 2 c_{A}n_{A}, \end{gather}
\begin{gather} \mu \boldsymbol{A}_{(1)} + A -n_{A} \boldsymbol{I} -\delta_{A} \nabla^{2}\boldsymbol{A} = c_{B} n_{B} \boldsymbol{B} - c_{A} \boldsymbol{A}, \end{gather}
\begin{gather} \epsilon \mu \boldsymbol{B}_{(1)} + B -\frac{n_{B}}{2} \boldsymbol{I} -\epsilon\delta_{B} \nabla^{2}\boldsymbol{B} ={-}2\epsilon c_{B} n_{B} \boldsymbol{B} + 2 \epsilon c_{A} \boldsymbol{A}. \end{gather}
Here the subscript 
$( )_{(1)}$ denotes the upper-convected derivative which is given as 
${\partial ()}/{\partial t} + \boldsymbol {U}\boldsymbol {\cdot } \boldsymbol {\nabla } () - ( (\boldsymbol {\nabla } \boldsymbol {U})^{T} \boldsymbol {\cdot } () + ()\boldsymbol {\cdot } \boldsymbol {\nabla } \boldsymbol {U})$. The non-dimensional parameters 
$\mu$, 
$\epsilon$ and 
$\delta _{A,B}$ are given as 
${\lambda _{A}}/{\lambda _{eff}}$, 
${\lambda _{B}}/{\lambda _{A}}$ and 
${\lambda _{A} D_{A,B}}/{d^{2}}$, respectively, where 
$\lambda _{B}$ is the relaxation time of the short chain 
$B$ and 
$D_{A, B}$ are the dimensional diffusivities of the long and short species A and B, respectively. Furthermore, according to the VCM model, the non-dimensional breakage rate 
$(c_{A})$ of the long chain A into two equally sized short chains B depends on the local state of the stress field, and it is given by the expression 
$c_{A} = c_{Aeq} + \mu ({\xi }/{3})( \dot {\boldsymbol {\gamma }}: {\boldsymbol {A}}/{n_{A}})$, whereas the reforming rate of the long chain A from the two short chains B is assumed to be constant which is given by the equilibrium reforming rate, i.e. 
$c_{B} = c_{Beq}$. Here the nonlinear parameter 
$\xi$ is the scission energy required to break a long micelle chain into two shorter chains. The significance of this parameter is that as its value increases, the amount of stress needed to break the chain increases.
2.3. Single-species Giesekus constitutive equation
In the single-species constitutive equation, the number density of the wormlike micelles remains constant due to the absence of breakage and reformation, and hence one does not need to solve any species conservation equation, as solved in the two-species VCM model. However, one has to solve the equation for the evaluation of the polymer conformation tensor (which is related to the stresses, as mentioned above) as follows:
\begin{equation} \boldsymbol{A}_{(1)} + \boldsymbol{A} - \boldsymbol{I} ={-}\alpha \left( \boldsymbol{A} - \boldsymbol{I} \right) \boldsymbol{\cdot} \left( \boldsymbol{A} - \boldsymbol{I} \right). \end{equation}
The dimensionless parameter 
$\alpha$ is known as the Giesekus mobility factor, and for 
$0 < \alpha <1$, the above equation is known as the Giesekus constitutive equation. This equation is derived based on the kinetic theory of closely packed polymer chains, and the mobility factor 
$\alpha$ is introduced in this model in order to take into account the anisotropic hydrodynamic drag on the polymer molecules (Giesekus Reference Giesekus1982).
3. Numerical details
All the governing equations, namely mass, momentum, Giesekus and VCM constitutive equations, have been solved using the finite volume method based open-source computational fluid dynamics code OpenFOAM (Weller et al. Reference Weller, Tabor, Jasak and Fureby1998). In particular, the recently developed rheoFoam solver available in rheoTool (Pimenta & Alves Reference Pimenta and Alves2016) has been used in the present study. A detailed discussion of the present numerical set-up and its validation has been presented in our recent studies (Khan & Sasmal Reference Khan and Sasmal2020; Sasmal Reference Sasmal2020), and hence it is not repeated here. The following boundary conditions were employed in order to solve the present problem. On the sphere surface, standard no-slip and no-penetration boundary conditions, for the velocity, i.e. 
$U_{x} = U_{y} = 0$, are imposed whereas a no-flux boundary condition is assumed for both the stress and micellar number density, i.e. 
$\boldsymbol {n} \boldsymbol {\cdot } \boldsymbol {\nabla } \boldsymbol {A} = \boldsymbol {n} \boldsymbol {\cdot } \boldsymbol {\nabla } \boldsymbol {B} = 0$ and 
$\boldsymbol {n} \boldsymbol {\cdot } \boldsymbol {\nabla } n_{A} = \boldsymbol {n} \boldsymbol {\cdot } \boldsymbol {\nabla } n_{B} = 0$. It should be mentioned here that micelles may undergo a slip flow at the sphere surface, particularly if the sphere surface is roughened in nature (Mohammadigoushki & Muller Reference Mohammadigoushki and Muller2018). However, in the present study, the sphere is assumed to be solid with a smooth surface, and hence the application of the no-slip boundary condition is justified at this stage. On the tube wall, 
$U_{x} = U$ and 
$U_{y} = 0$, and again no-flux boundary conditions for the stress and micellar number density are imposed. At the tube outlet, a Neumann type of boundary condition is applied for all variables except for the pressure for which a zero value is assigned here. A uniform velocity of 
$U_{x} = U$, a zero gradient for the pressure and a fixed value for the micellar number density are employed at the tube inlet. Furthermore, the whole computational domain was subdivided into eight blocks in order to mesh it, as shown in figure 1(b). Three different meshes of hexagonal block structure, namely M1 , M2 and M3 with different numbers of cells on the sphere surface 
$(N_{s})$ as well as in the whole computational domain 
$(N_{t})$, were created for each blockage ratio. A schematic of three different mesh densities is shown in figure 1(b) for 
$BR = 0.33$. In creating any mesh density, the cells were further compressed towards the sphere surface in order to capture the steep gradients of velocity, stress and micellar concentration. After performing the standard mesh independent study, the mesh M2 (with 
$N_{s} = 240$ and 
$N_{t} = 74\,200$ for 
$BR = 0.33$ and 
$N_{s} = 240$ and 
$N_{t} = 78\,600$ for 
$BR = 0.1$) was found to be adequate to capture the flow physics for the whole range of conditions encompassed here for both blockage ratios. Similarly, a time step size of 
$\Delta t = 0.000055$ was found to be suitable to carry out the present study. Finally, the two-dimensional and axisymmetric problem is realized in OpenFOAM by applying the standard wedge boundary condition (with wedge angle less than 
$5^{\circ }$) on the front and back surfaces of the computational domain, as schematically shown in figure 1(c). The computational domain is kept one cell thick in the 
$\theta$ direction, and the axis of the wedge lies on the 
$x$ coordinate, as per the requirement for applying the wedge boundary condition (OpenFOAM 2020). Such simplification does not compromise the accuracy of the results as long as the flow is two-dimensional and axisymmetric, and also markedly reduces the computational cost compared to a full three-dimensional simulation.
4. Results and discussion
The VCM model parameters chosen in the present study are: 
$\mu = 5.7$, 
$\epsilon = 4.5 \times 10^{-4}$, 
$\beta _{VCM} = 6.8 \times 10^{-5}$, 
$\xi = 0.7$, 
$n_{B}^{0} = 1.13$, 
$\delta _{A} = \delta _{B} = \delta = 10^{-3}$. These values are obtained by fitting the experimental results for small-amplitude oscillatory shear and step-strain experiments for a mixture of cetylpyridinium chloride/NaSal added to water (Pipe et al. Reference Pipe, Kim, Vasquez, Cook and McKinley2010; Zhou, McKinley & Cook Reference Zhou, McKinley and Cook2014). The rheological characteristics of the WLM solution with these parameter values in homogeneous shear and uniaxial extensional flows are shown in figure 2. One can clearly see in this figure the two typical properties of a WLM solution, namely the shear thinning in shear flows and the extensional hardening and subsequent thinning in extensional flows. Additionally, the present WLM solution also shows a shear-banding phenomenon. The corresponding parameters for the single-species Giesekus model are chosen as 
$\beta _{G} = 4.98 \times 10^{-3}$ and 
$\alpha = 0.2$ and 0.8. For the single-species model, at 
$\alpha =0.8$, the solution shows the shear banding and extensional thinning properties.Simulations were carried out for a shear Weissenberg number 
$(Wi_{S})$ of up to 2 for both the two-species VCM and single-species Giesekus models in the creeping flow regime. Up to a shear Weissenberg number of 0.6 (not shown here), the streamlines are attached to the sphere surface, and they follow an orderly path without crossing each other for both the single-species Giesekus and two-species VCM models. Hence there is a perfect fore-and-aft symmetry present in the streamline patterns, and also the flow remains steady up to this value of the Weissenberg number.
Figure 2. Variations of the non-dimensional shear stress with the shear rate (a) and non-dimensional first normal stress difference with the extension rate (b) in homogeneous shear and extensional flows, respectively. The insets show the corresponding variations in the shear and extensional viscosities.
However, as the Weissenberg number gradually starts to increase, clear differences are observed in the flow patterns obtained with the Giesekus and VCM models. For instance, at 
$Wi_{{S}} = 2.0$ (figure 3a), the streamlines are still attached to the sphere surface for the single-species Giesekus model, thereby suggesting no boundary layer separation happens for this model. Furthermore, the flow remains in the steady state at this value of the Weissenberg number. This is confirmed by plotting the temporal variation of the streamwise velocity at a probe location downstream of the sphere 
$(X = 1.0, Y = 0)$, as shown in figure 3(e). One can see that the velocity reaches a steady value with time. On the other hand, at the same Weissenberg number, for the two-species VCM model, separation of the boundary layer happens, and as a result, a small recirculation region is seen to form downstream of the sphere (figure 3b). As time further progresses, the wake detaches from the sphere surface, and its size becomes small (figure 3c), and ultimately it disappears as can be seen in figure 3(d). This formation and disappearance of the wake is seen to repeat with time. It should be mentioned here that a weak recirculation region was seen beside the sphere but not downstream of it for the falling-sphere problem (Chen & Rothstein Reference Chen and Rothstein2004). All these observations suggest that the flow becomes unsteady at this Weissenberg number. This is further confirmed by plotting the streamwise velocity in figure 3(e) at the same probe location as that obtained for the Giesekus model, and one can clearly see how the velocity fluctuates with time. In fact, the unsteadiness in the flow field appears at a much lower Weissenberg number of approximately 
$Wi_{S} = 0.6$ for the VCM model. Furthermore, a region of very high velocity magnitude is seen to appear at about one sphere diameter away downstream of the sphere at 
$t = 50.7$ (figure 3b). This indicates the presence of a ‘negative wake’ downstream of the sphere. As time gradually increases, the magnitude of this region progressively decreases, and it is further shifted downstream of the sphere (figure 3c). Finally, it has vanished at a time 
$t = 50.9$ (figure 3d). On further increasing the time, it appears again, and these processes of appearance and disappearance of the negative wake downstream of the sphere repeat with time. This was also observed in experiments for the falling-sphere problem by Chen & Rothstein (Reference Chen and Rothstein2004).
Figure 3. Representative streamlines and velocity magnitude plots for the Giesekus model (a) and for the VCM model at three different times, namely 
$t = 50.7$ (b), 
$t = 50.8$ (c) and 
$t = 50.9$ (d) at 
$Wi_{S} = 2.0$. Temporal variation of the streamwise velocity component 
$(U_{X})$ for both Giesekus and VCM models (e). Power spectrum plots of the velocity fluctuations obtained with the VCM model at 
$Wi_{S} = 1.0$ (f) and 
$Wi_{S} = 2.0$ (g).
To further analyse the nature of the unsteadiness in the flow field, the power spectrum of these velocity fluctuations is plotted in figures 3(f) and 3(g) at shear Weissenberg numbers 1.0 and 2.0, respectively, for the VCM model. At 
$Wi_{S} = 1.0$, the velocity fluctuations are characterized by a single and large narrow peak in the frequency spectrum, thereby suggesting the flow field to be a periodic one. On the other hand, at 
$Wi_{S} = 2.0$, the frequency spectrum of the velocity fluctuations is characterized by a broad range with significant contribution from higher-order frequencies (figure 3g). This suggests that the flow becomes quasi-periodic at this value of the Weissenberg number. Therefore, it can be seen that as the shear Weissenberg number gradually increases, the flow transits from steady to periodic and then from periodic to quasi-periodic for the two-species VCM model. This transition in the flow regime with the shear Weissenberg number is entirely in line with that observed experimentally for the falling-sphere problem by Zhang & Muller (Reference Zhang and Muller2018). Their experimental observations further found an irregular flow pattern after the quasi-periodic flow regime on further increasing the Weissenberg number. However, due to the numerical stability problem at high Weissenberg number, we were unable to run simulations at Weissenberg numbers beyond 2, and hence this irregular flow pattern was not observed in our simulations within the ranges of conditions encompassed in this study.
Next, we explore the reason behind this unsteady motion, which is observed for the two-species VCM model but not seen for the single-species Giesekus model over the ranges of conditions considered in this study. The explanation goes as follows. At high values of the Weissenberg number, for instance at 
$Wi_{S} = 2.0$, and at 
$t = 50.7$, the long micelles tend to break into smaller ones downstream of the sphere due to the presence of high flow strength in this region. This can be seen in figure 4(a) where the surface plot of number density of long micelles is presented at the same three times as used for showing the streamline and velocity magnitude plots in figure 3. To present it more quantitatively, the variation of number density of long micelles along the downstream axis is plotted in figure 4(d). It can be seen that the number density of long micelles shows a minimum at approximately 
$X = 0.5d$. The corresponding variation of the streamwise velocity component along the downstream axis of the sphere is depicted in figure 4(e). Due to the breakage of long micelles, the extensional load, which was carried by the long micelles, is unable to be carried by the small micelles. This results in the formation of a small recirculation region downstream of the sphere. Due to the formation of this recirculation region, the velocity first gradually decreases, attains a minimum and then gradually starts to increase. It then shows a peak and ultimately reaches the velocity far away downstream of the sphere. This stage of the unsteady motion was termed as the acceleration stage in case of the falling-sphere problem (Chen & Rothstein Reference Chen and Rothstein2004; Mohammadigoushki & Muller Reference Mohammadigoushki and Muller2016). During this stage, a negative wake was formed, as shown in figure 3(b). The occurrence of a velocity overshoot in figure 4(e) further confirms the existence of this negative wake downstream of the sphere. The region of velocity overshoot is seen just beside the region where the concentration of long micelles is minimum downstream of the sphere. As time further progresses, new long micelles come into the downstream region, and the extensional stresses again start to develop. The velocity downstream of the sphere gradually starts to decrease, and the position at which the velocity changes its sign further moves downstream, as can be seen in figure 4(e) at 
$t = 50.8$. Ultimately, at 
$t = 50.9$, the velocity gradually reaches the far-downstream velocity without showing any minimum or peak in its profile. This means that at this time, no recirculation region as well as no negative wake are present in the flow field (which are also evident in the streamline and velocity magnitude plots shown in figure 3), and under these conditions, the fluid behaves like a highly elastic Boger fluid. This was called the deceleration stage in the unsteady motion for the falling-sphere problem (Chen & Rothstein Reference Chen and Rothstein2004; Mohammadigoushki & Muller Reference Mohammadigoushki and Muller2016). Due to the deceleration in the flow field at later times, the breakage of long micelles also decreases downstream of the sphere, as can be seen in figure 4(e). On further increasing the time, the acceleration stage again appears, and it repeats with time. Therefore, this study provides evidence that the acceleration and deceleration motion past a steadily translating sphere in WLM solutions is solely caused by the breakage of long micelles downstream of the sphere. This is further confirmed by the fact that the single-species Giesekus model shows a steady flow field under otherwise identical conditions. This explanation presented herein was also proposed by Chen & Rothstein (Reference Chen and Rothstein2004) in their experimental investigation of the unsteady motion of a falling sphere in WLM solutions. The present study further confirms their hypothesis about the onset of this instability using the two-species VCM model.
Figure 4. Surface plot of the number density of long micelles at 
$Wi_{S} = 2.0$ and at three different times, namely 50.7 (a), 50.8 (b) and 50.9 (c). Variation of the number density of long micelles (d) and streamwise velocity component along the downstream axis of the sphere (e).
Some studies associated with the falling-sphere problem revealed that the onset of the sphere's unsteady motion is directly associated with the strong extensional flow behaviour downstream of the sphere, which ultimately leads to the breakage of long micelles (Mohammadigoushki & Muller Reference Mohammadigoushki and Muller2016; Zhang & Muller Reference Zhang and Muller2018). Therefore, we also calculate the extension rate along the downstream axis of the steadily translating sphere, defined as 
$\dot {\epsilon }_{XX} = {\partial U_{X}}/{\partial X}$. This is presented in figure 5(a) again at three different times, namely 
$t = 50.7$, 50.8 and 50.9, and at 
$Wi_{S} = 2.0$. One can see that there is a large temporal variation present in the extension rate downstream of the sphere up to a distance of around 
$X = 1d$ from the rear stagnation point of the sphere, and beyond that it becomes almost zero. At 
$t = 50.7$, the variation in the strain rate is seen to be higher as compared to that seen at later times. At this time, the maximum breakage of long micelles occurs downstream of the sphere, as shown in figure 4(a). The temporal variation of the maximum strain rate 
$(\dot {\epsilon }_{M})$ along the downstream axis of the sphere is shown in figure 5(b). It can be clearly seen that the maximum strain rate varies quasi-periodically, and a large variation in its value is present. This is again, at least qualitatively, in line with that observed experimentally for the falling-sphere problem (Mohammadigoushki & Muller Reference Mohammadigoushki and Muller2016; Zhang & Muller Reference Zhang and Muller2018). Similar to the experimental investigations of the falling-sphere problem, we also define an extensional Weissenberg number based on the time-averaged values of this maximum strain rate as 
$Wi_{{Ext}} = \lambda _{eff} \dot {\epsilon }_{M}$. The variation of this extensional Weissenberg number with the corresponding shear Weissenberg number is shown in figure 5(c). From this figure, it is seen that the value of the extensional Weissenberg number increases with the shear Weissenberg number, and the transition from steady to unsteady periodic is marked by a sharp increase in the value of the extensional Weissenberg number, which is again in line with the corresponding experimental observations for the falling-sphere problem (Mohammadigoushki & Muller Reference Mohammadigoushki and Muller2016; Zhang & Muller Reference Zhang and Muller2018).
Figure 5. (a) Variation of the extension rate along the downstream axis of the sphere at three different times and at 
$Wi_{{S}} = 2.0$. (b) Temporal variation of the maximum extension rate downstream of the sphere at 
$Wi_{{S}} = 2.0$. (c) Variation of the extensional Weissenberg number 
$(Wi_{{Ext}})$ versus 
$(Wi_{{S}})$.
Finally, the effect of the ratio of sphere diameter to tube diameter on the onset and generation of this unsteady motion past the sphere is discussed. Simulations were carried out at another value of 
${d}/{D} = 0.1$ in order to compare with flow characteristics as discussed above at 
${d}/{D} = 0.33$ under otherwise identical conditions. Figure 6(a) shows the temporal variation of the streamwise velocity at two ratios of sphere to tube diameter at a probe location (
$X = 1, Y = 0$) downstream of the sphere and at a shear Weissenberg number of 
$Wi_{S} = 2.0$. One can see that the velocity field again shows a quasi-periodic nature at 
${d}/{D} = 0.1$ similar to that seen at 
${d}/{D} = 0.33$. However, the magnitude of the velocity during the acceleration stage decreases, whereas it increases during the deceleration stage. Moreover, the magnitude of the velocity fluctuations slightly decreases with decreasing values of the ratio of sphere to tube diameter. This is clearly evident in the power spectrum plot (figure 6b) wherein the amplitude of the maximum frequency spectrum is seen to be slightly higher at 
${d}/{D} = 0.33$ than that seen at 
${d}/{D} = 0.1$. At 
${d}/{D} = 0.1$, a similar transition in the velocity field is seen to that observed at 
${d}/{D} = 0.33$, i.e. it transits from steady to unsteady periodic to unsteady quasi-periodic upon increasing the Weissenberg number. However, the onset of the unsteady motion is slightly delayed as the ratio of sphere to tube diameter decreases. For instance, at 
${d}/{D} = 0.33$, the unsteady motion starts at 
$Wi_{S} = 0.6$ whereas it starts at around 
$Wi_{S} = 1.0$ for a ratio of sphere to tube diameter of 0.1.The reason behind this can be explained as follows. As the ratio of sphere to tube diameter decreases, the extensional flow strength downstream of the sphere also decreases due to the presence of less confinement. This can be seen in figure 6(c) wherein the temporal variation of the maximum extension rate downstream of the sphere is seen to be less at 
${d}/{D} = 0.1$ in comparison with that seen at 
${d}/{D} = 0.33$. As a result, the time-averaged extensional Weissenberg number is also found to be less for the former case in comparison with the latter under otherwise identical conditions (figure 6d). This lowering in the extensional flow strength downstream of the sphere tends to decrease the breakage of the micelles in this region, which in turn delays the tendency of appearance of the unsteady motion. This also further confirms the hypothesis that the unsteady motion past a sphere translating steadily in a micellar solution is caused by the breakage of stretched micelles downstream of the sphere.
Figure 6. Effect of the ratio of sphere diameter to tube diameter or blockage ratio (BR) on the temporal variation of the streamwise velocity (a), power spectrum plot (b) and the temporal variation of the maximum extension rate (c) downstream of the sphere at 
$Wi_{S} = 2.0$. Effect of BR on the variation of the extensional Weissenberg number with the shear Weissenberg number (d).
5. Conclusions
This study presents an extensive numerical investigation of the flow characteristics past a sphere translating steadily along the axis of a cylindrical tube filled with a WLM solution. For doing so, the present study uses the two-species VCM and single-species Giesekus constitutive models for representing the rheological behaviour of WLM solutions. Over the ranges of conditions encompassed in this study, a transition of the flow field downstream of the sphere from a steady to unsteady periodic to unsteady quasi-periodic regime is seen as the shear Weissenberg number gradually increases. A similar transition in the velocity field was also observed in experiments on the sedimentation of a sphere in WLM solutions. The onset of this unsteady motion is marked by a steep increase in the value of the extensional Weissenberg number, defined downstream of the sphere based on the maximum extension rate, once again in accordance with the experiments on the falling-sphere problem. Due to this increase in the extensional flow strength downstream of the sphere, the breakage of long micelles occurs, thereby causing the unsteady motion in the flow field downstream of the sphere. This is further confirmed as the single-species Giesekus model for the WLM solution predicts a steady velocity field under otherwise identical conditions. This explanation is in line with that proposed in the earlier experimental investigations reported in the literature for the sedimentation of a sphere in WLM solutions. Furthermore, it is seen that the onset of this unsteady motion is delayed as the ratio of sphere diameter to tube diameter decreases due to the decrease in the extensional flow strength downstream of the sphere.
Although a very good qualitative agreement is found between the present numerical predictions for the steadily translating sphere and the experimental findings for the sedimentation of a sphere, it should be mentioned here that the present simulation set-up does not mimic the exact experimental settings for the falling-sphere problem. In the experiments on the falling-sphere problem, the sphere may rotate or undergo lateral motion or even may not reach a terminal velocity. Therefore, to realize the experimental settings for the falling-sphere problem, the governing equations, namely continuity, momentum and micelle constitutive equations, need to be solved in full three-dimensional numerical settings along with an equation for the sphere motion. In contrast to this, in the present study, the sphere is assumed to be translating steadily along the axis of the tube, and the problem is solved in a coordinate system that is centred on and travelling with the sphere. This could be a situation in the corresponding experiments on the falling-sphere problem when the sphere reaches its terminal velocity without rotation and lateral motion. Although this is hardly a case in actual experiments, with this simplified problem, for the first time, we have provided evidence that the unsteady motion past a sphere is caused by the breakage of long micelles, resulting from an increase in the extensional flow strength downstream of the sphere. We believe that the analysis shown in this study will further support the hypothesis presented earlier for the unsteady motion of a falling sphere in micellar solutions. In our future study, we plan to carry out full three-dimensional numerical simulations with exactly the same settings as those used for the experiments for the sedimentation of spheres in micellar solutions. This will give us the opportunity to conduct a more accurate and quantitative comparison with the experiments carried out for the falling-sphere problem.
Acknowledgements
The author would like to thank Mr A. Chauhan for his help in the meshing of the geometry in OpenFOAM.
Funding
The author would like to thank IIT Ropar for providing funding through the ISIRD research grant (Establishment1/2018/IITRPR/921) to carry out this work.
Declaration of interests
The author reports no conflict of interest.
