1. Introduction
Mass transfer in two-phase systems has several applications in the chemical engineering field, such as the design of efficient and sustainable reactors for the production of pharmaceutical and agrochemical compounds. The development of continuous flow reactors, which are characterised by a continuous flow of reactants and products, has recently attracted the attention of researchers due to promising performance compared with standard batch devices. A successful design of continuous flows reactors based on the Taylor–Couette (TC) flow was proposed recently for electrochemical (Love et al. Reference Love, Lee, Gennari, Jefferson-Loveday, Pickering, Poliakoff and George2021; Lee et al. Reference Lee, Love, Mansouri, Waldron Clarke, Harrowven, Jefferson- Loveday, Pickering, Poliakoff and George2022) and photochemical (Lee et al. Reference Lee, Amara, Clark, Xu, Kakimpa, Morvan, Pickering, Poliakoff and George2017, Reference Lee, Sharabi, Jefferson-Loveday, Pickering, Poliakoff and George2020) applications involving both single- and two-phase (gas–liquid) reactions. The success of such design is mainly due to the excellent mixing properties of TC flows and the optimal bubble size distribution within the reaction vessel.
A TC (Couette Reference Couette1890; Taylor Reference Taylor1923) apparatus consists of two coaxial rotating cylinders and the flow behaviour within the gap is a well-studied configuration that exhibits several consecutive states during the transition from the laminar regime (low rotating speeds) to a fully turbulent flow. In the last few decades, TC flow has captured the attention of both scientists active in the study of laminar to turbulent transition (see, for example, Gollub & Swinney Reference Gollub and Swinney1975; Smith & Townsend Reference Smith and Townsend1982; Townsend Reference Townsend1984) as well as engineers involved in the design of rotating devices, such as rotating machinery (Nicoli, Johnson & Jefferson-Loveday Reference Nicoli, Johnson and Jefferson-Loveday2022) or chemical reactors (Schrimpf et al. Reference Schrimpf, Esteban, Warmeling, Fárber, Behr and Vorholt2021). Extensive literature has been published on the characterisation of TC flows and the interested reader is referred to the works of Di Prima & Swinney (Reference Di Prima and Swinney1981), Andereck, Liu & Swinney (Reference Andereck, Liu and Swinney1986), Wang (Reference Wang2015), Grossmann, Lohse & Sun (Reference Grossmann, Lohse and Sun2016) and the references therein for a detailed review. In this work, only the configuration where the inner cylinder is rotating and the outer one is kept stationary is considered, but similar behaviours can be observed in the more generic case of counter-rotating walls.
The majority of studies concerning disperse bubbly flows in TC apparatuses is mainly devoted to the analysis of drag reduction mechanisms (such as bubble deformability or effective compressibility of the flow, Ferrante & Elghobashi Reference Ferrante and Elghobashi2004) on the rotating walls (Murai, Oiwa & Takeda Reference Murai, Oiwa and Takeda2005; Van den Berg et al. Reference Van den Berg, Luther, Lathrop and Lohse2005; Murai, Oiwa & Takeda Reference Murai, Oiwa and Takeda2008; Sugiyama, Calzavarini & Lohse Reference Sugiyama, Calzavarini and Lohse2008; Van Gils et al. Reference Van Gils, Narezo Guzman, Sun and Lohse2013; Murai Reference Murai2014; Wang, Wang & Liu Reference Wang, Wang and Liu2022), as well as bubble accumulation patterns and their interaction with the flow structures (Shiomi et al. Reference Shiomi, Kutsuna, Akagawa and Ozawa1993; Djeridi, Gabillet & Billard Reference Djeridi, Gabillet and Billard2004; Murai et al. Reference Murai, Oiwa and Takeda2005; Climent, Simonnet & Magnaudet Reference Climent, Simonnet and Magnaudet2007; Mehel, Gabillet & Djeridi Reference Mehel, Gabillet and Djeridi2007; Ymawaki et al. Reference Ymawaki, Hosoi, Takigawa, Noui-Mehidi and Ohmura2007; Chouippe et al. Reference Chouippe, Climent, Legendre and Gabillet2014; Gao, Kong & Vigil Reference Gao, Kong and Vigil2015b, Reference Gao, Kong and Vigil2016). On the other hand, a comprehensive understanding of gas–liquid mass transfer in TC flows is missing in the literature and the available studies mainly focus on the experimental quantification of mass transfer coefficients (
$k_m$) through the measurement of dissolved gaseous concentration in the liquid solution (Ramezani et al. Reference Ramezani, Kong, Gao, Olsen and Vigil2015; Qiao et al. Reference Qiao, Yan, Teoh, Tong and Wang2018). Since the interfacial gas–liquid area (
$a_\varSigma$) is difficult to measure experimentally, the product 
$k_m a_\varSigma$ is generally provided instead and correlation formulae of the type 
$Sh=f(Re, Sc)$ are proposed, where 
$Sh$ is the Sherwood number and 
$Sc$ is the Schmidt number. The Sherwood number compares the mass transfer coefficient against the characteristic velocity of diffusion and is defined as 
$Sh = {k_m L_{ref}}/{D_c}$ (where 
$L_{ref}$ is the reference length and 
$D_c$ is the diffusion coefficient of the gaseous species in the liquid (continuous) phase), whereas the Schmidt number is 
$Sc={\nu _c}/{D_c}$, and 
$\nu _c$ is the kinematic viscosity of the liquid. Gao et al. (Reference Gao, Kong, Ramezani, Olsen and Vigil2015a) combine a theoretical model, based on the penetration theory of Higbie (Reference Higbie1935), with an Euler–Euler numerical framework to quantify mass transfer in TC reactors. The Euler–Euler approach does not resolve the gas–liquid interface and allows disperse bubbly flows to be modelled in large domains. However, its accuracy relies on the choice of appropriate closure models, which typically depend on the specific application for the exchange of interfacial mass and momentum. To the best of the authors’ knowledge, no studies have been published on the modelling of bubbles in TC flows by means of fully resolved interfacial simulations. The present study (of which part of the material is based on one of the present author's thesis, Gennari Reference Gennari2023) contributes to the understanding of gas–liquid mass transfer in TC flows by deploying a fully resolved and state-of-the-art numerical volume-of-fluid (VOF) framework to capture both the fluid flow and mass transfer occurring at the interface and investigate how bubble dissolution is affected by the different regimes of TC flows. In the rest of this section, the main features of TC flows are briefly discussed along with a review of the available numerical methodologies for fully resolved two-phase flows with mass transfer.
The non-dimensional groups generally used for the characterisation of this flow configuration take into account both the geometry of the apparatus (see figure 1), which is defined by the radius ratio 
$(\eta = {r_{in}}/{r_{out}} )$ and the aspect ratio 
$( \varGamma = {L_z}/ (r_{out} - r_{in}) )$, as well as the Reynolds number 
$Re = {\rho _c U_{in} (r_{out} - r_{in} )}/{\mu _c}$, where the subscript 
$c$ is used to refer to the liquid (continuous phase) within the reactor. The inner and outer radii are 
$r_{in}$ and 
$r_{out}$, respectively, whereas 
$L_z$ is the axial extension of the device and 
$U_{in}=r_{in} \omega _{in}$ is the peripheral speed of the inner rotor.
Figure 1. Geometrical parameters of a TC apparatus and representation of counter-rotating Taylor vortices.
The first instability that occurs in a (planar, time-independent and axisymmetric) cylindrical Couette flow, when the rotating speed exceeds a critical value, consists of pairs of counter-rotating vortices (also known as Taylor cells) superimposed on the main flow; this flow regime is referred to as Taylor vortex flow (TVF) and the cells have a characteristic toroidal-like shape. The flow is periodic in the axial direction, axisymmetric and time-independent. The Reynolds number at which this instability occurs is referred to as critical Reynolds (
$Re_{cr}$) and the expected wavelength 
$\lambda$ (i.e. the axial extension of two consecutive Taylor cells; see figure 1) is approximately twice the gap between the cylinders. As the rotating speed is further increased beyond the critical Reynolds, a second instability is observed, which causes the vortices to travel along the azimuthal direction, following a wavy trajectory. The boundaries between two adjacent Taylor cells have a sinusoidal shape (wave) and the flow is no longer time-independent. The waves are periodic along the azimuthal direction and this configuration is referred to as wavy vortex flow (WVF). In this regime, the flow can exhibit multiple states, i.e. different number of Taylor cells and azimuthal waves for the same 
$Re$ (Coles Reference Coles1965). A third instability occurs for larger Reynolds numbers and it is characterised by the appearance of two sharp frequencies in the power spectra of the velocity field. The first is still associated with the travelling of azimuthal waves (as in the WVF regime), whereas the second is related to a modulation of the amplitude and the frequency of such waves (Gorman & Swinney Reference Gorman and Swinney1982). This configuration is generally referred to as modulated wavy vortex flow (MWVF). For the first three regimes (i.e. TVF, WVF and MWVF), Koschmieder (Reference Koschmieder1979) reported that the axial wavelength increases with the rotating speed up to approximately 
$Re = 10Re_{cr}$ (for an apparatus with 
$\eta =0.896$), after which 
$\lambda$ is found to be independent of the rotating speed. For 
$Re>10Re_{cr}$ the azimuthal waves progressively disappear and the flow transitions towards a turbulent regime. This is the last state of TC flow and is generally referred to as turbulent Taylor vortex flow (TTVF). From visual observations, the flow is still structured into azimuthal cells, although the velocity field is no longer well-organised into a toroidal pattern, due to the presence of strong velocity fluctuations.
When a soluble gas is introduced in a liquid solution, the system reaches an equilibrium state where part of the gas is dissolved into the liquid according to the partial pressure exerted by the gas on the interface between the phases. In the present work, the interface is always assumed saturated and Henry's law (see § 2.1) is used to compute the concentration jump between the disperse phase (i.e. the gas) and the continuous phase (i.e. the liquid). Whenever the dissolved concentration in the continuous phase (
$c_{bulk}$) is below the interfacial saturated value 
$( c_c )_\varSigma$, i.e. the saturation ratio 
$\zeta = c_{bulk}/( c_c )_\varSigma < 1$, a diffusion-driven process that depends on the local concentrations at the interface (Groß & Pelz Reference Groß and Pelz2017) takes place and redistributes gas molecules from the disperse phase into a concentration boundary layer 
$\delta _c$ on the liquid side of the interface, leading to bubble dissolution. Assuming a uniform concentration within 
$\delta _c$ and no species initially dissolved in the continuous domain (i.e. 
$c_{bulk}=0$), the concentration boundary layer thickness can be estimated as 
$\delta _c = D_b/Sh$, where 
$D_b$ is the bubble diameter. In actual cases of dissolving rising bubbles, 
$\delta _c$ is a local quantity that varies around the interface and is determined by an advection–diffusion process. The relative importance of these two transport mechanisms is estimated by the Péclet number, defined as 
$Pe=Re_bSc$. For large Reynolds number, Levich (Reference Levich1962) used the potential flow theory to approximate the flow field around a moving spherical particle and derived the well-known formula 
$Sh=( 2/\sqrt {{\rm \pi} } ) \sqrt {Pe}$. A similar functional relationship is also found in other theoretical formulations, such as Oellrich, Schmidt-Traub & Brauer (Reference Oellrich, Schmidt-Traub and Brauer1973), as well as experimental correlation models (Takemura & Yabe Reference Takemura and Yabe1998). By combining this relationship with the hydrodynamic boundary layer theory 
$\delta _h \approx D_b/\sqrt {2Re_b}$ (Levich Reference Levich1962), the ratio of concentration to hydrodynamic boundary layer thicknesses evolves as 
$\delta _c/\delta _h \propto 1/\sqrt {Sc}$ (Weiner & Bothe Reference Weiner and Bothe2017).
One of the limiting factors of fully resolved numerical simulations of interfacial flows with mass transfer is due to the small scales that occur at large 
$Sc$ and 
$Pe$ numbers; additional challenges are given by the discontinuities that characterise the interface in terms of velocity (whenever mass transfer between two phases with different density occurs) and concentration of soluble species. Different approaches have been developed in the past to address these points, such as neglecting volume changes for highly dilute species (Bothe & Fleckenstein Reference Bothe and Fleckenstein2013; Farsoiya, Popinet & Deike Reference Farsoiya, Popinet and Deike2021) or smearing the interfacial mass transfer term to improve stability (Hardt & Wondra Reference Hardt and Wondra2008). Other methodologies adopt the ghost fluid method (Fedkiw et al. Reference Fedkiw, Aslam, Merriman and Osher1999) to deal with interfacial velocity jumps (Nguyen, Fedkiw & Kang Reference Nguyen, Fedkiw and Kang2001; Sussman Reference Sussman2003; Tanguy, Ménard & Berlemont Reference Tanguy, Ménard and Berlemont2007; Tanguy et al. Reference Tanguy, Sagan, Lalanne, Couderc and Colin2014), whereas recent works have focused on techniques to derive a divergence-free velocity formulation at the interface to advect the indicator function in a VOF framework (Guo Reference Guo2020; Scapin, Costa & Brandt Reference Scapin, Costa and Brandt2020; Malan et al. Reference Malan, Malan, Zaleski and Rousseau2021; Gennari, Jefferson-Loveday & Pickering Reference Gennari, Jefferson-Loveday and Pickering2022; Boyd & Ling Reference Boyd and Ling2023; Cipriano et al. Reference Cipriano, Saufi, Frassoldati, Faravelli, Popinet and Cuoci2024). Specific numerical schemes have been developed to preserve the jump between the concentration values at the interface and can be divided into two families, namely one-scalar (one transport equation per species) and two-scalar (two equations per species) methods. One-scalar approaches include the work of Bothe et al. (Reference Bothe, Koebe, Wielage, Prúss and Warnecke2004) and are further extended in (Haroun, Legendre & Raynal Reference Haroun, Legendre and Raynal2010; Marschall et al. Reference Marschall, Hinterberger, Schúler, Habla and Hinrichsen2012; Deising, Marschall & Bothe Reference Deising, Marschall and Bothe2016; Maes & Soulaine Reference Maes and Soulaine2018); examples of these methods coupled with algebraic VOF frameworks can be found in Maes & Soulaine (Reference Maes and Soulaine2020) for competing mass transfer and in Zanutto et al. (Reference Zanutto, Evrard, van Wachem, Denner and Paladino2022a,Reference Zanutto, Paladino, Evrard, van Wachem and Dennerb) for evaporating flows and non-ideal mixtures. Two-scalar approaches are presented in Alke et al. (Reference Alke, Bothe, Kroeger and Warnecke2009), Bothe & Fleckenstein (Reference Bothe and Fleckenstein2013) and used in Fleckenstein & Bothe (Reference Fleckenstein and Bothe2015) with a geometric VOF for multicomponent mass transfer with volume effects. A novel implementation is presented in Schulz et al. (Reference Schulz, Wecker, Inguva, Lopatin and Kenig2022), where the mesh is split at the interface based on its geometrical reconstruction. A combination of one- and two-scalar schemes is presented in Farsoiya et al. (Reference Farsoiya, Popinet and Deike2021); the same authors have recently proposed an alternative approach that takes into account volume effects (Farsoiya et al. Reference Farsoiya, Magdelaine, Antkowiak, Popinet and Deike2023). In the present work, a geometric VOF scheme is adopted and the piecewise linear reconstruction (PLIC) of the interface allows for a sharp separation between the disperse and continuous domains. Under these circumstances, a two-scalar method is the preferred choice to prevent any artificial mass transfer to occur during the advection of the interface (Deising et al. Reference Deising, Marschall and Bothe2016).
The rest of this article is organised as follows. The governing equations for two-phase flows with soluble species are introduced in § 2.1, whereas the numerical methodology, which is based on our previous work (Gennari et al. Reference Gennari, Jefferson-Loveday and Pickering2022), is briefly summarised in § 2.2. The numerical framework is validated extensively in 3, whereas the results of bubble dissolution in TC flows are discussed in § 4. It is finally recalled here that the terms continuous (disperse) and liquid (gas) are used interchangeably in the rest of the work.
2. Governing equations and numerical framework
2.1. Governing equations
In this work, the three-dimensional (3-D) Navier–Stokes equations for a two-phase incompressible flow with phase-change are solved in the one-fluid framework (see Tryggvason, Scardovelli & Zaleski Reference Tryggvason, Scardovelli and Zaleski2011 for a rigorous derivation):
$$\begin{gather} \partial_t H+\boldsymbol{\nabla}\boldsymbol{\cdot}(H\boldsymbol{u}) ={-}\frac{\dot{m}}{\rho_c}\delta_\varSigma , \end{gather}$$
$$\begin{gather}\boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{u} = \dot{m}\left(\frac{1}{\rho_d}-\frac{1}{\rho_c}\right)\delta_\varSigma, \end{gather}$$
$$\begin{gather}\partial_t\boldsymbol{u}+\boldsymbol{\nabla}\boldsymbol{\cdot}(\boldsymbol{u}\otimes\boldsymbol{u}) = \frac{1}{\rho}[-\boldsymbol{\nabla} p + \boldsymbol{\nabla}\boldsymbol{\cdot}(2\mu\boldsymbol{\mathsf{S}})] + \boldsymbol{g}+\frac{\sigma \kappa \boldsymbol{n}_\varSigma}{\rho} \delta_\varSigma , \end{gather}$$where (2.1) represents the transport of the Heaviside function, which is used to mark the location of the interface between the continuous and disperse domains,
\begin{equation} H(\boldsymbol{x},t) = \begin{cases} 1, & \text{if } \boldsymbol{x}\in \varOmega_c,\\ 0, & \text{if } \boldsymbol{x} \in \varOmega_d. \end{cases} \end{equation}
Once 
$H(\boldsymbol {x},t)$ is known everywhere, the values of density (
$\rho$) and viscosity (
$\mu$) can be computed as
\begin{equation} \rho=\rho_c H+\rho_d (1-H) \end{equation}and
\begin{equation} \mu=\mu_c H+\mu_d (1-H), \end{equation}
where the subscript 
$c$ 
$(d)$ is used to refer to the continuous (disperse) phase. In the following we use the letter 
$f_c$ to refer to the volume fraction of the continuous phase in a computational cell with volume 
$V$, i.e. 
$f_c = ({1}/{V}) \int _V H \,{\rm d}V$. Equations (2.2)–(2.3) are the balances of mass and momentum, respectively, where the term on the right-hand side of the continuity equation takes into account volume effects when phases with different densities exchange mass. In the system of (2.1)–(2.3), 
$\boldsymbol {u}$ is the velocity field, 
$\dot {m}$ is the mass transfer rate, 
$\delta _\varSigma$ is the interfacial Dirac function, 
$p$ is the pressure, 
$\boldsymbol{\mathsf{S}}$ is the deformation tensor 
$[\boldsymbol {\nabla }\boldsymbol {u} + (\nabla \boldsymbol {u})^{\rm T}]/2$, 
$\boldsymbol {g}$ is the gravitational acceleration, 
$\sigma$ is the surface tension, 
$\kappa$ and 
$\boldsymbol {n}_\varSigma$ are the curvature and the normal vector of the interface.
Mass transfer at the interface of a two-phase system can occur for different physical phenomena, such as evaporation, boiling, chemical reactions and gas solubility. In the present work, the focus is on the solubility of gaseous species in liquid solutions, where the mass transfer is driven by a diffusive process that occurs at the interface (diffusion-driven phase-change) and depends on the species concentration around the interface (
$\varSigma$). Therefore, to close the system of governing equations, the conservation law for soluble species in two-phase flows needs to be included. This takes the form of a system of two transport equations for the molar concentration field (
$c^k$) of each soluble component in the domain and, for the generic 
$k$th species, reads (see Bothe & Fleckenstein Reference Bothe and Fleckenstein2013)
\begin{equation} \left.\begin{array}{c} \partial_t c_c^k + \boldsymbol{u}_c \boldsymbol{\cdot} \boldsymbol{\nabla} c_c^k - \boldsymbol{\nabla} \boldsymbol{\cdot} ( D_c^k\boldsymbol{\nabla} c_c^k ) ={-}\dfrac{\dot{m}^k}{M^k}\delta_\varSigma \quad \text{in } \varOmega_c,\\ \partial_t c_d^k + \boldsymbol{u}_d \boldsymbol{\cdot} \boldsymbol{\nabla} c_d^k - \boldsymbol{\nabla} \boldsymbol{\cdot} ( D_d^k\boldsymbol{\nabla} c_d^k ) =\dfrac{\dot{m}^k}{M^k}\delta_\varSigma \quad \text{in } \varOmega_d, \end{array}\right\} \end{equation}
where the subscripts 
$c,d$ emphasise that the equations of system (2.7) must be integrated in the respective domain only, i.e. 
$\varOmega _{c,d}$; 
$M^k$ and 
$D^k_{c,d}$ are the molar mass and diffusivity in phase (
$c,d$) of the species. The species mass transfer term that appears in the concentration transport equations (2.7), is, by definition, the difference between the species and interface velocities along the normal direction. For the generic 
$k$th component, it reads
\begin{equation} \|\rho^k(\boldsymbol{u}^k-\boldsymbol{u}_\varSigma)\boldsymbol{\cdot}\boldsymbol{n}_\varSigma\| = \|\dot{m}^k\|=0 \quad \text{at } \varSigma, \end{equation}
where the jump notation has been introduced (i.e. 
$\|\rho ^k\|=\rho ^k_c-\rho ^k_d$). Equation (2.8), also known as the Rankine–Hugoniot condition, implies that no mass can be stored at the interface. Here, a generic system of 
$n$ components is considered, where the first 
$n-1$ elements are soluble species (that can be transferred across the interface and appear as dilute components in the liquid phase), and the 
$n$th component is the solvent, which is assumed to be not volatile (i.e. no solvent species exists in the disperse phase). Under these assumptions, the mass transfer rate of a single species can be rearranged into (see Fleckenstein & Bothe Reference Fleckenstein and Bothe2015 for more details)
\begin{equation} \dot{m}^k = \frac{\rho^k}{\rho} \sum_{l=1}^{n-1} \dot{m}^l - D^k \boldsymbol{\nabla}\rho^k \boldsymbol{\cdot} \boldsymbol{n}_\varSigma, \end{equation}
where (2.9) can be computed from either the continuous or disperse side of the interface. A special case arises when the disperse phase is made of a single species only (i.e. no mixtures). In this case, the system contains two components (
$n=2$): the pure gas (
$k=1$) and the solvent (liquid phase, 
$k=2$) and the overall mass transfer (
$\dot {m}$) is entirely given by the transfer rate of the single species which exists in the disperse phase, i.e. 
$\dot {m}=\dot {m}^1$. The mass transfer rate for a pure disperse phase can be computed from (2.9) as
\begin{equation} \dot{m}={-}\frac{D^1_c}{1-y^1_c}\frac{\partial \rho^1_c}{\partial\boldsymbol{n}_\varSigma}, \end{equation}
where 
$y_c^1=\rho ^1/\rho _c$, whereas the subscript c has been added to denote that the mass transfer rate must be computed from the liquid side of the interface (computing 
$\dot {m}$ from (2.9) in 
$\varOmega _d$ gives the identity 
$\dot {m} = \dot {m}$).
One more condition needs to be taken into account at the interface for the chemical partitioning of species densities. In a generic two-phase flow, the species distribution at the interface is discontinuous and, for a system at equilibrium (saturated interface), such jump in the concentration profile can be predicted by Henry's law, which states that the 
$k$th species concentration on the liquid side of 
$\varSigma$ is directly proportional to the partial pressure of the same gaseous species on the liquid. By taking advantage of the perfect gas law, Henry's formula can be written in terms of a jump condition for the species densities at the interface:
\begin{equation} (c^k_c)_\varSigma=\frac{(c^k_d)_\varSigma}{H_e^k}, \end{equation}
where 
$H_e^k$ is the Henry's law coefficient for the 
$k$th species and it is a material property of the system, which generally depends on the temperature and pressure fields near the interface (see Bothe & Fleckenstein Reference Bothe and Fleckenstein2013 for a detailed discussion about the generalised Henry's law). For the applications considered in the present work, 
$H_e^k$ is assumed to be constant for each species and the interface is always treated as saturated.
2.2. Numerical framework
The governing equations presented in § 2.1 are solved with the open-source solver Basilisk (Popinet & collaborators Reference Popinet2013–2024). Basilisk is a finite-volume solver for the solution of partial differential equations on adaptive Cartesian grids and implements a second-order accurate (time and space discretisation) solver for direct numerical simulations (DNS) of two-phase immiscible fluids (Popinet Reference Popinet2009). The interface position is tracked with a geometric VOF method and state-of-the-art numerical techniques are implemented for the computation of the interface curvature, which is particularly relevant to mitigate the numerical effect of spurious currents (Popinet Reference Popinet2009). The Cartesian mesh is organised into a hierarchical tree structure (Popinet Reference Popinet2015) and can be dynamically adapted (i.e. refined and/or coarsened) by means of an adaptive mesh refinement (AMR) technique based on a wavelet estimation of the spatial discretisation error for selected flow fields (van Hooft et al. Reference van Hooft, Popinet, van Heerwaarden, van der Linden, de Roode and van de Wiel2018). The ability of adapting the mesh at each iteration in regions where strong gradients occur makes Basilisk an efficient solver for interfacial flows, where generally a fine mesh is required around the gas–liquid interface and a coarser discretisation can be employed for the remaining part of the domain. In this work, we adopt the phase-change solver presented in Gennari et al. (Reference Gennari, Jefferson-Loveday and Pickering2022) and implemented in Basilisk. In the following, the main ingredients of the numerical algorithms are briefly summarised.
The integration of (2.1) is performed in two steps: first, the advection term is integrated with the PLIC scheme presented in Weymouth & Yue (Reference Weymouth and Yue2010) (based on an operator-split method), then the interface is shifted with a rigid displacement along the normal direction, equivalent to 
$\boldsymbol {h}_\varSigma = -({\dot {m}}/{\rho _c}) ({\Delta t}/{\Delta }) \boldsymbol {n}_\varSigma$. This last term corresponds to the integration of the source term on the right-hand side of (2.1). The VOF scheme is designed to ensure mass conservation for incompressible flows without phase-change and relies on the kinematic constraint 
$\boldsymbol {\nabla } \boldsymbol{\cdot} \boldsymbol {u} = 0$. In case of mass transfer occurring at the interface, the non-divergence-free condition (2.2) introduces a velocity discontinuity that no longer satisfies the conservation of mass. To address this problem, a novel algorithm was proposed in Gennari et al. (Reference Gennari, Jefferson-Loveday and Pickering2022), which consists of a redistribution of the mass transfer term 
$\dot {m}$ from the interfacial cells to a layer of pure gas cells next to the interface. The redistributed term is then used for the numerical discretisation of the continuity equation (2.2), which produces a divergence-free velocity field in both liquid and interfacial cells.
To prevent artificial mass transfer during the integration of the species transport equations (2.7), both advection and diffusion terms must transport the molar concentration in their respective phase only, i.e. no transfer of moles across the interface is allowed at this stage. This is accomplished by advecting the molar concentration with the same geometric fluxes (based on the PLIC reconstruction of the interface) used for the transport of the Heaviside function. For the generic 
$k$th species, the flux reads (López-Herrera et al. Reference López-Herrera, Gañán-Calvo, Popinet and Herrada2015; see also figure 2a)
\begin{equation} F^{{adv},k}_{p,x (i-1/2,j)} = \frac{\Delta V_p}{\Delta t} c^k_{p(i-1/2,j)}, \quad \text{for } p=c,d, \end{equation}
where 
$\Delta V_p$ is the exact (in the sense of the PLIC reconstruction of the interface) amount of volume of phase 
$p$ that crosses the cell edge. The molar concentration on the face is predicted using the upwind scheme of Bell, Colella & Glaz (Reference Bell, Colella and Glaz1989), which performs an extrapolation in time (half-time-step) and in space from the upwind cell centre to the cell boundary. At this point, a correction for the advection of 
$c$ in 
$\varOmega _d$ is required, since the velocity field is no longer divergence-free in the disperse domain near the interface. This is accomplished here using the same approach adopted in Fleckenstein & Bothe (Reference Fleckenstein and Bothe2015), where the global dilation term is subtracted after all the one-dimensional (1-D) advection operations are performed.
Figure 2. (a) Advection of species concentrations confined within the respective phases. The transport fluxes across the cell boundary are based on the PLIC advection of the respective volume of fluids (red and green volumes for the continuous and disperse phases, respectively); 
$u_f$ represents the face-centred velocity field. (b) Unsplit scheme for the computation of the mass transfer term.
The diffusion term is treated with the approach proposed in López-Herrera et al. (Reference López-Herrera, Gañán-Calvo, Popinet and Herrada2015) and Magdelaine-Guillot de Suduiraut (Reference Magdelaine-Guillot de Suduiraut2019), which is equivalent to a standard finite-volume scheme, where the fluxes across the boundaries are computed on all the cell faces and the diffusion coefficient is multiplied by the face fraction (obtained from the PLIC reconstruction) of the respective phase. For the generic 
$k$th species, the flux reads
\begin{equation} F^{diff,k}_{p,x (i-1/2,j)} = \frac{\partial c^k_p}{\partial \boldsymbol{n}} ( D_p^k f_{f,p} ) A, \quad \text{for } p=c,d, \end{equation}
where 
$f_{f,p}$ is the face fraction on the cell boundary of phase 
$p$, i.e. 
$f_{f,p} = A_p/A$, and 
$A$ is the area of the cell face (figure 2a). The gradients along the Cartesian axes are computed with a central finite difference scheme.
Finally, the mass transfer term 
$\dot {m}$ requires the evaluation of the gradient term 
$\boldsymbol {\nabla } c^k \boldsymbol{\cdot} \boldsymbol {n}_\varSigma$ (2.9). This is calculated here from the continuous side, by using the unsplit scheme proposed by Bothe & Fleckenstein (Reference Bothe and Fleckenstein2013) and it reads (see figure 2b)
\begin{equation} {-}\frac{\partial c^k_{c(i,j)}}{\partial \boldsymbol{n}_\varSigma} = f_c\frac{c^k_c(P_1)-c^k_c(P)}{\overline{PP_1}} + (1 - f_c)\frac{c^k_c(P_2)-c^k_c(P)}{\overline{PP_2}}, \end{equation}
where the values of concentration in points 
$P_1$ and 
$P_2$ are obtained from quadratic (bi-quadratic in three dimensions) interpolation, whereas the value at the centroid of the interface 
$c_c^k(P)$ is computed by applying Henry's law (2.11).
3. Validation of the numerical framework
3.1. Single-phase TC flow
In this section, the Basilisk code is validated for single-phase TC flows against available data in the literature. The DNS of (3-D) incompressible flows are performed and wall boundaries are treated with an embedded boundary method, where Dirichlet boundary conditions are enforced with the approach proposed by Schwartz et al. (Reference Schwartz, Barad, Colella and Ligocki2006). The tangential velocity 
$U_{in}=r_{in} \omega _{in}$ is applied at the inner cylinder, whereas the outer one is fixed (i.e. 
$U_{ext}=0$) and periodic boundary conditions are used for the top and bottom ends of the computational domain (see figure 1). The choice of the axial length of the domain (
$L_z$) is particularly relevant when only a section of the apparatus is modelled, since periodic boundaries force the flow to adapt to the available space and constrain the number of Taylor vortices that form within the annulus. Results from linear stability analysis for infinite cylinders (see the appendix by P.H. Roberts in Donnelly et al. Reference Donnelly, Schwarz, Roberts and Chandrasekhar1965) show that the wavelength, i.e. the axial extension of a pair of counter-rotating vortices (see figure 1), is expected to be close to 
$\lambda \approx 2 ( r_{out} - r_{in} )$. However, the results collected in the work of Chouippe et al. (Reference Chouippe, Climent, Legendre and Gabillet2014) from different experimental investigations show that a significant dispersion is observed in the measured wavelengths. The main reason is due to the non-uniqueness feature of the TC flow for which the final observed state of the system depends on the procedure used to reach such state (e.g. acceleration/deceleration rates of the rotor) and not only on the geometrical configuration. Therefore, for the validation of the numerical method, it is important to select an axial length that is a multiple of the observed wavelength (i.e. 
$L_z = n \lambda$), so that a number of 
$n$ vortex pairs is modelled and a sensible comparison can be made against the reference data. In the present work, three configurations are tested, namely 
$\eta = 0.5$, 0.73 and 0.91, at different Reynolds numbers. Details on the main parameters, including the observed wavelength and the critical Reynolds number (
$Re_{cr}$) for the transition from planar Couette flow to TVF, are summarised in table 1 (for a comprehensive summary on the critical values for a range of radius ratios, the reader is referred to Childs Reference Childs2011 and the references therein). The selected choice of configurations allows for a comprehensive validation of the single-phase numerical framework, since the main TC regimes are represented (i.e. TVF, WVF and TTVF). For details on the mesh sensitivity study and characteristics of the selected grids for fully resolved simulations, the reader is referred to Appendix A.
Table 1. Single-phase TC cases.
The cases reported in table 1 are run until an equilibrium configuration is reached and the flow statistics are stationary. This state occurs when the torque exerted by the fluid on the walls is the same for both the inner and outer cylinders (Chouippe et al. Reference Chouippe, Climent, Legendre and Gabillet2014) and an example of the plot of the non-dimensional torque (
$G_w$) for the configuration 
$\eta =0.5$, 
$Re=5000$ is reported in figure 3. The torque is made non-dimensional with the cylinders' axial length and with the liquid density and viscosity:
\begin{equation} G_w^{{in,out}}=\frac{T_w^{{in,out}}}{\rho_c \nu_c^2 L_z}. \end{equation}
The mean torque values for all the tested configurations at their equilibrium points are compared against the experimental formula proposed by Wendt (Reference Wendt1933), where 
$G_w$ scales as 
$Re^{3/2}$:
\begin{equation} G_w^{{Wendt}} = 1.45 \left[\frac{\eta^{3/2}}{(1 - \eta)^{7/4}} \right]Re^{3/2} \end{equation}and the corresponding results are reported in figure 4, where, for all the simulated cases, a good comparison against the experimental data is observed, confirming that the statistically stationary regime is reached for all the tested radius ratios and Reynolds numbers.
Figure 3. Inner and outer cylinder (non-dimensional) torques vs time for the TC configuration with 
$\eta = 0.5$ and 
$Re = 5000$. The absolute value 
$|Gw|$ is plotted here to compare between the two walls. The statistically stationary regime is approximately reached after 50 revolutions.
Figure 4. Comparison of the (non-dimensional) torque exerted on the inner cylinder against the experimental work of Wendt (Reference Wendt1933) (3.2).
The mean azimuthal velocity 
$\langle u_\theta \rangle _{z \theta t}$ and fluctuation 
$\sqrt {\langle {u_\theta '}^2\rangle _{z \theta t}}$ (see Appendix A for their derivation) for the configurations with 
$\eta =0.5$ and 
$\eta =0.91$ are compared against the available numerical data of Dong (Reference Dong2007) and Chouippe et al. (Reference Chouippe, Climent, Legendre and Gabillet2014) and results are reported in figures 5 and 6, respectively. A good comparison is observed for almost all the selected configurations, for both the average azimuthal velocity and the corresponding fluctuation. The profiles of velocity fluctuations show the characteristic shape with two local peaks near the inner and outer walls and an (almost) uniform value in the bulk of the liquid; similar profiles are observed for different turbulent channels configurations (Moser & Moin Reference Moser and Moin1987; Hoyas & Jiménez Reference Hoyas and Jiménez2006). As the Reynolds number increases, the magnitude of the (normalised) fluctuations decreases and the peaks move closer to the respective walls. The configuration with 
$\eta =0.5$, 
$Re=1000$ shows a significant deviation for the azimuthal fluctuation (but not for the main velocity component) from the work of Dong (Reference Dong2007) (figures 5a and 6a). However, the same case compared to the results reported in Chouippe et al. (Reference Chouippe, Climent, Legendre and Gabillet2014) for 
$u'_\theta$ shows an excellent agreement at every distance from the walls. Surprisingly, the radial profile of average azimuthal velocity for the configuration with 
$\eta =0.5$, 
$Re=3000$ (figure 5b) does not match the reference data of Dong (Reference Dong2007) within the bulk of the liquid, where the velocity is under-predicted, but a good agreement is reached in the regions close to the inner and outer walls. However, the same configuration agrees well with both the works of Chouippe et al. (Reference Chouippe, Climent, Legendre and Gabillet2014) and Dong (Reference Dong2007) in terms of velocity fluctuations (figure 6b), although some quantitative discrepancies with the latter reference near the inner wall are observed.
Figure 5. Average radial profiles of the azimuthal velocity component for the configurations with (a) 
$\eta = 0.5$, 
$Re = 1000$, (b) 
$\eta = 0.5$, 
$Re = 3000$, (c) 
$\eta = 0.5$, 
$Re = 5000$ and (d) 
$\eta = 0.91$, 
$Re = 5000$.
Figure 6. Average radial profiles of the azimuthal velocity fluctuation for the configurations with (a) 
$\eta = 0.5$, 
$Re = 1000$, (b) 
$\eta = 0.5$, 
$Re = 3000$ and (c) 
$\eta = 0.5$, 
$Re = 5000$.
A qualitative representation of the flow field and the effect of the Reynolds number for the configurations with 
$\eta =0.5$ and 
$\eta =0.91$ is reported in figure 7, where the contours of axial velocity (
$u_z$) on a cylindrical surface with constant 
$r$ are compared in a (planar) two-dimensional (2-D) plot on the corresponding 
$z$–
$\theta$ plane. Figure 7(a–c) shows the effect of the Reynolds number on the topology of Taylor vortices as the flow regime evolves from WVF to TTVF (see table 1). For 
$Re=1000$ (figure 7a) two organised pairs of counter-rotating vortices develop within the annulus and a thin region of null axial velocity separates each vortex from the adjacent (counter rotating) one. The axial extension of the computational domain was set to twice the expected wavelength (see table 1) and the qualitative results reported here (two pairs of vortices) confirm that the axial length of Taylor cells matches the expected one. The travelling trajectory along the azimuthal direction of each vortex is almost straight, but the onset of a wavy motion is visible from the oscillating boundaries of the vortices, suggesting that the apparatus is in a transitional state from TVF to WVF. As the Reynolds number is increased to 
$Re=3000$ (figure 7b), the flow is fully turbulent and the shape of the vortices is distorted. However, two main regions of counter rotating velocities can still be identified, although Taylor cells are not well-defined as in the case with 
$Re=1000$. Finally, for 
$Re=5000$ (figure 7c) the flow appears chaotic with many flow structures distributed in a random way and Taylor vortices do not form into an organised and clear pattern; these observations are qualitatively confirmed by the results reported in Dong (Reference Dong2007). The effect of the gap size is clearly visible from the comparison between figure 7(c) (
$\eta =0.5$) and figure 7(d) (
$\eta =0.91$), which both run at 
$Re=5000$. For larger radius ratios, the small gap within the cylinders represents a geometric constraint for the formation of Taylor vortices, whose topology appears (even for large and fully turbulent Reynolds numbers) well organised into stable and clearly recognisable pairs of alternating axial velocities.
Figure 7. Contours of axial velocity on the 
$z$ - 
$\theta$ plane for the configurations with (a) 
$\eta = 0.5$, 
$Re = 1000$, (b) 
$\eta = 0.5$, 
$Re = 3000$ and (c) 
$\eta = 0.5$, 
$Re = 5000$ and (d) 
$\eta = 0.91$, 
$Re = 5000$. These plots are obtained from the corresponding cylindrical surface with radius 
$r_{in} + 0.1(r_{out} - r_{in})$ for cases (a–c) and radius 
$r_{in} + 0.25(r_{out} - r_{in})$ for case (d).
The results presented in this section show that the numerical methodology used in the present work to model single-phase TC flows is able to accurately reproduce the features of the main flow regimes for different geometries (radius ratios) and rotating speeds (Reynolds numbers).
3.2. Phase-change solver
In this section the numerical framework presented in § 2.2 is validated for the generic scenario of competing mass transfer of a mixture of soluble species. The concentration of species is non-uniform in both phases and the direction of mass transfer, i.e. from 
$\varOmega _d$ to 
$\varOmega _c$ or vice-versa, can be different for each component, depending on the local concentration at the interface.
3.2.1. Mass transfer in an infinite cylinder
In this test case, a binary gaseous mixture made of two soluble components (species A and B) is confined by a liquid annulus where 
$R_{in}$ and 
$R_{ext}$ are the inner and outer radius, respectively. The liquid phase is therefore confined within the region 
$R_{in} < r < R_{ext}$, whereas the gaseous phase exists for 
$r < R_{in}$. The axial length of the cylinder (
$L_z$) is infinite and the external radius is set to 
$R_{ext} = 1$ mm. The inner radius of the liquid annulus, which represents the interface between the phases, is free to move as some of the species crosses the interface and is initially set to 
$R^{t=0}_{in} = 0.5$ mm. Due to the infinite axial extension, the problem is independent of the axial coordinate and can be represented by a 2-D model; a sketch of the computational domain is shown in figure 8.
Figure 8. Computational domain for an infinite gaseous cylinder (
$\varOmega _d$) confined by a liquid annulus (
$\varOmega _c$).
The properties of the gas–liquid system are reported in table 2 and approximate an air–water system. The case simulated in this section replicates one of the set-ups proposed in Maes & Soulaine (Reference Maes and Soulaine2020), where the gaseous (disperse) phase is initially composed of species B only, i.e. 
$c_d^{B(t=0)}=\rho _d/M^B$. Species A is assumed to be weakly soluble in the liquid solvent, whereas species B is not soluble and the respective Henry's law coefficients are 
$H_e^A=100$ and 
$H_e^B \rightarrow \infty$. By setting Henry's law coefficient to 
$H_e^B \rightarrow \infty$ for species B, the equilibrium value on the liquid side of the interface is 
$( c_c^B )_\varSigma = 0$, regardless of the amount of species within the gaseous domain. Since no species B exists initially in the liquid domain, the mass transfer of B across the interface is prevented (i.e. the solution is saturated with respect to species B), and the species is confined within the gaseous region. The liquid domain is therefore composed of the solvent (not soluble in the disperse phase) and species A, which has a relatively (compared with a typical gas solubility) large Henry's law coefficient and, therefore, is weakly soluble in the liquid solvent; the concentration of species A is kept constant at the external boundary (
$r = R_{ext}$) and set to 
$c_c^A(R_{ext},t) = \rho _d/( M^A H_e^A )$. Diffusivity is the same for both species and is set to 
$D_c^A = D_c^B = 10^{-6}$ m
$^{2}$ s
$^{-1}$ and 
$D_d^A = D_d^B = 10^{-4}$ m
$^{2}$ s
$^{-1}$ in the continuous and disperse phases, respectively.
Table 2. Gas–liquid properties for competing mass transfer in an infinite cylinder.
Due to the symmetry of the problem, the velocity and concentration fields depend only on the radial distance (and time), and the liquid moves along the radial direction only, i.e. 
$\boldsymbol {u}_c = u_c(r,t) \boldsymbol {e}_r$. Under this assumption, the problem can be simplified significantly and the following analytical model is derived (see Maes & Soulaine Reference Maes and Soulaine2020 for the details):
$$\begin{gather} c_c^A(r,t) = \frac{\rho_d}{M^A H_e^A} \left( 1 - \frac{R_{in}^{2(t=0)}}{R_{in}^2(t)} \frac{\ln ( r/R_{ext} )}{\ln ( R_{in}(t)/R_{ext} )} \right), \quad \text{for } r > R_{in}(t), \end{gather}$$
$$\begin{gather}\frac{{\rm d}R_{in}}{{\rm d}t} = \frac{D_c^A R_{in}^{2(t=0)}}{H_e^A R^3_{in}(t) \ln ( R_{in}(t)/R_{ext} )}, \quad \text{for } t>0. \end{gather}$$ A summary of the numerical set-up is given in table 3; the mesh size is set to 
$\varDelta = 1.95 \times 10^{-5}$ mm, whereas the molar masses are the same for both species and equal to 
$M^A = M^B = 1$ kg mol
$^{-1}$. The concentration profile of species A in 
$\varOmega _c$ is initialised with (3.3) at 
$t=0$, coherently with the assumption of solution at equilibrium at every time step. Results are made non-dimensional with the reference length 
$L_{ref} = R_{ext}$, time 
$t_{ref} = \rho _c R^2_{ext}/ \mu _c$ and concentration 
$c_{ref} = \rho _d/M^A$, whereas the reference velocity follows from 
$U_{ref} = L_{ref}/t_{ref}$. The numerical simulation is run for a time of 
$\Delta t=5$ s and the result in terms of interface position (
$R_{in}$) is compared against the analytical solution (3.4) in figure 9, where a good agreement is observed.
Table 3. Numerical set-up for a cylinder of gas expanding in an infinite liquid annulus.
Figure 9. Inner radius of the liquid annulus vs time.
3.2.2. Competing mass transfer in a rising bubble
This benchmark is based on the test case proposed in Fleckenstein & Bothe (Reference Fleckenstein and Bothe2015) and consists of the study of competing mass transfer amongst three soluble species for a bubble rising in a quiescent flow. The properties of the gas–liquid system used for the present test case are reported in table 4. The soluble species that exist in the present model are: CO
$_{2}$, N
$_{2}$ and O
$_{2}$; the respective properties are reported in table 5. The main non-dimensional numbers used for the present analysis are the bubble Reynolds number 
$(Re_b = {\rho _c U_b D_b}/{\mu _c} )$, Galilei 
$( \sqrt {{g D_b^3}/{\nu _c^2}} )$, Bond 
$( {\rho _c g D_b^2}/{\sigma } )$, Schmidt 
$( Sc^k={\nu _c}/{D_c^k} )$ and Péclet 
$( Pe^k = Re_bSc^k )$. In these groups, the index 
$k$ refers to the generic 
$k$th component, whereas 
$U_b$ is the bubble rising velocity.
Table 4. Gas–liquid properties for the competing mass transfer in a rising bubble.
Table 5. Species properties for the competing mass transfer in a rising bubble.
In order to speed up the volume change process and reduce the computational time of the simulation, the diffusivity for all the species in the liquid domain 
$( D_c^k )$ has been increased by a factor of 10 with respect to the real physical property (the same approach was used in the reference case of Fleckenstein & Bothe Reference Fleckenstein and Bothe2015); the corresponding diffusivity in the disperse phase 
$( D_d^k )$ is assumed to be 100 times larger than the continuous one (i.e. 
$D_d^k = D_c^k \times 10^2$). The solubility of CO
$_{2}$ is significantly larger than the solubility of the other species (lower Henry's law coefficient), which means that, for the same concentrations in both phases, the mass transfer from the gaseous region to the liquid (under-saturated solutions) occurs faster for CO
$_{2}$ than N
$_{2}$ and O
$_{2}$; the opposite scenario occurs for super-saturated solutions, where the transfer from the continuous phase to the liquid phase is quicker for N
$_{2}$ and O
$_{2}$ than CO
$_{2}$. In table 5, the Schmidt numbers are computed with the liquid properties reported in table 4 and are similar for all species, since the diffusivity of each component does not change significantly.
The initial diameter of the bubble is set to 
$D_b^{t=0}={0.8}$ mm and the bubble is confined in a large square domain with dimensions 
$L0 \times L0 ={48}\ \textrm {mm} \times {48}\ \textrm {mm}$, where it rises under the effect of the gravitational field 
$g=9.81$ m s
$^{-2}$. Due to the large dimension of the domain compared to the bubble size, end walls effect do not affect the dynamics of the bubble in the present case. The Galilei and Bond numbers are 
$Ga = 79.39$ and 
$Bo = 0.0869$, respectively, and, for these parameters, the bubble is expected to rise vertically, keeping the original spherical shape. Therefore, a 2-D axisymmetric model is used here, where only half of the bubble is considered, and the rising trajectory is the horizontal 
$x$-axis, i.e. 
$\boldsymbol {g} = -g\boldsymbol {e}_x$. An outflow condition is applied to the right boundary to allow the liquid to enter/leave the domain as the bubble volume changes, whereas symmetric conditions are used for the other boundaries; AMR is used to keep the grid at the finest level around the bubble and save computational cells far from the interface. Results are made non-dimensional with the reference length 
$L_{ref}=R_b^{t=0}$, time 
$t_{ref} = \sqrt {L_{ref}/g}$ and the gaseous concentration in 
$\varOmega _d$ when the bubble is composed of CO
$_{2}$ only, i.e. 
$c_{ref} = \rho _d/M^{{CO}_2}$. The bubble is initially composed of CO
$_{2}$ (i.e. 
$c_d^{{CO}_2(t=0)} = 44.59$ mol m
$^{-3}$), whereas the liquid solution is composed of the solvent (not soluble in 
$\varOmega _d$) and species N
$_{2}$ and O
$_{2}$ with concentrations 
$c_c^{{N}_2(t=0)} = 0.51$ mol m
$^{-3}$ and 
$c_c^{{O}_2(t=0)} = 0.27$ mol m
$^{-3}$. The solution is therefore under-saturated for CO
$_{2}$ and super-saturated for the other species.
A mesh sensitivity study is first performed to evaluate the level of grid refinement that is necessary to reach a mesh-independent solution. Four grids are tested (cases A, B, C and D) and the mesh size around the interface, along with the number of cells per diameter of the bubble, is summarised in table 6. The simulations are run for a time interval of 
$\Delta t = 0.25$ s and results in terms of volume change for the bubble are shown in figure 10. The grid convergence analysis shows that a mesh independent solution is reached for case C, which corresponds to approximately 546 cells per diameter at 
$t=0$. For the selected chemical composition of the liquid and gaseous phases, CO
$_{2}$ is transferred from the bubble to the liquid (under-saturation), whereas N
$_{2}$ and O
$_{2}$ flow in the opposite direction (super-saturation). Due to the larger solubility of CO
$_{2}$ compared with the other species and the weak super-saturation ratios for N
$_{2}$ and O
$_{2}$, the competing mass transfer is dominated by CO
$_{2}$ and results in a net flow of mass out of the bubble; the phase volume decreases accordingly. The volume reduces almost linearly in the first part of the simulation (until 
$t^*\approx 10$), where the mass transfer is driven by CO
$_{2}$ and the concentration of N
$_{2}$ and O
$_{2}$ are still marginal. As the chemical composition inside the bubble changes and the mass fractions of N
$_{2}$ and O
$_{2}$ become more relevant, the volume change rate decreases and becomes almost negligible for 
$t^* > 30$. Since the solution does not change significantly between grids C and D, case C is used in the following part of the analysis.
Table 6. Grid convergence study for the competing mass transfer in a rising bubble. The mesh size 
$\varDelta$ refers to the maximum refinement around the interface, whereas the number of cells per diameter is computed assuming a uniform resolution inside the bubble.
Figure 10. Grid convergence for the competing mass transfer in a rising bubble. Plot of bubble volume vs time.
The maximum Péclet number is observed at 
$t^*\approx 5$ for CO
$_{2}$ and is approximately 
$Pe^{{CO}_2} \approx 7800$. The results in terms of grid sensitivity are consistent with the analysis performed in Gennari et al. (Reference Gennari, Jefferson-Loveday and Pickering2022) (see § 4.6) for pure bubbles rising at different Péclet numbers in an under-saturated solution, where a resolution of 
$456\ \textrm {cells}/D_b$ was required to reach a mesh-independent solution at 
$Pe=4650$. Results in terms of chemical composition of the bubble are shown in figure 11 for case C (case A is discussed later in the text). The bubble is initially composed of CO
$_{2}$ only, therefore the mass fractions are 
$m_d^{{CO}_2(t=0)} = 1$ and 
$m_d^{{N}_2(t=0)} = m_d^{{O}_2(t=0)} = 0$. As the phase-change process occurs, CO
$_{2}$ is transferred to the liquid, whereas the other species flow across the interface in opposite directions; the mass fraction of CO
$_{2}$ decreases, whereas the fractions of the other species increase accordingly. Due to the lower solubility of N
$_{2}$ compared with O
$_{2}$ and larger initial concentration in the liquid phase, the mass fraction of N
$_{2}$ grows faster than O
$_{2}$ and reaches the same value of the fraction of CO
$_{2}$ at 
$t^* \approx 34.8$ and becomes the most relevant component of the bubble for 
$t^*>34.8$. The fraction of O
$_{2}$ equals CO
$_{2}$ at 
$t^*\approx 37.7$ and CO
$_{2}$ becomes the most marginal species at the end of the simulation. The sum of the mass fractions is reported in figure 11, which shows that the method is mass conservative since the global mass fraction is always 
$m_d^{{tot}} = m_d^{{CO}_2} + m_d^{{N}_2} + m_d^{{O}_2} = 1$ for 
$t>0$.
Figure 11. Species mass fractions and bubble volume vs time. Results from case A are compared against the work of Fleckenstein & Bothe (Reference Fleckenstein and Bothe2015), where a similar mesh resolution is adopted.
To validate the accuracy of the numerical methodology, results are compared with the work of Fleckenstein & Bothe (Reference Fleckenstein and Bothe2015), where the set-up for this case was taken from. In the reference work, the mesh density corresponds to approximately 102 cells per (initial) diameter, which is similar to the grid refinement used for case A in the present work (see table 6). Case A is therefore used for the comparison against the reference case and results in terms of volume and mass fractions of the bubble are reported in figure 11, where a good agreement is observed for all the plotted quantities.
4. Bubble dissolution in TC flow
4.1. Simulation set-up and governing parameters
In this section, a single (pure) gas bubble is injected at the bottom of a TC device and is let free to exchange mass with the surrounding liquid. The selected apparatus for this study is that with radius ratio of 
$\eta =0.5$ for Reynolds numbers in the range 
$0 \leq Re \leq 5000$, which was validated extensively for the single-phase case in § 3.1. The properties of the gas–liquid system are reported in table 7. The initial bubble diameter is set to 
$D_b^{t=0} = (r_{out} - r_{in})/3 = {5}$ mm and the centre of the bubble is placed in the middle of the gap at 
$z_b^{t=0} = r_{out}/3$ from the bottom of the device (it is recalled here that the axis of the apparatus is aligned to the 
$z$ direction), whereas the solution is assumed initially under-saturated, with no concentration of gas at 
$t=0$ in the continuous phase (i.e. 
$c_c^{t=0} = 0$ mol m
$^{-3}$). The cylinders are oriented vertically and standard gravitational acceleration is assumed here, i.e. 
$\boldsymbol {g} = {-9.81}$ m s
$^{-2}$ 
$\boldsymbol {e}_{\boldsymbol {z}}$. Overall, the problem is defined by 13 dimensional parameters: apparatus radii (
$r_{in}, r_{out}$), rotor angular speed (
$\omega _{in}$), gravitational acceleration (
$\boldsymbol {g}$), densities (
$\rho _c, \rho _d$), viscosities (
$\mu _c, \mu _d$), initial bubble diameter (
$D_b^{t=0}$), surface tension (
$\sigma$), diffusion coefficient of the gaseous species in the liquid phase (
$D_c$) and initial species concentrations 
$(c_c^{t=0}, c_d^{t=0} )$. Given the four units involved (i.e. length, time, mass and amount of substance (mole)), the system can be described by nine independent non-dimensional numbers, reported, along with their definition and values, in table 8.
Table 7. Gas–liquid properties for a dissolving bubble in a TC flow.
Table 8. Independent non-dimensional numbers for a dissolving bubble in a TC flow.
Mass transfer is characterised by the analysis of the (time-dependent) Sherwood number, the definition of which follows as
\begin{equation} Sh=\frac{k_m(t) L_{ref}(t)}{D_c} \end{equation}
and it depends on a reference length, computed here as the equivalent (time-dependent) diameter of a sphere with the same volume as the bubble, i.e. 
$L_{ref} (t) = 2(3V_b/(4 {\rm \pi}))^{1/3}$. The mass transfer coefficient is based on the reference concentration difference 
$\Delta c$ between the continuous and disperse phases:
\begin{equation} k_m ={-}\frac{\displaystyle\int_\varSigma \dot{m} \,{\rm d}S}{A_\varSigma M \Delta c}, \end{equation}
where the reference area 
$A_\varSigma$ is the effective area of the interface. Other useful non-dimensional parameters can be derived from those reported in table 8, such as the bubble Morton number 
$Mo = Bo^3/Ga^4 = 3.22 \times 10^{-11}$ and Péclet number 
$Pe = Re_b Sc$, where 
$Re_b$ is the rising bubble Reynolds number. Finally, we introduce here the Froude number that will be used later to compare the effects of the inertial features of TC flows and gravity:
\begin{equation} Fr = \frac{u^{TC}}{\sqrt{gD_b}}, \end{equation}
where 
$u^{TC}$ is a characteristic velocity of TC flow.
Simulations are first started from rest (null liquid velocity) and the bubble is kept fixed until a (statistically) stationary regime is reached (see § 3.1) and Taylor vortices are completely formed. During this initialisation stage, the transport of species and volume change are not computed, i.e. the bubble volume remains constant. After the TC regime is established, the bubble is set free to move within the device and the full phase-change solver is run. As the volume of the bubble decreases, more liquid needs to be introduced within the apparatus for the conservation of mass. However, the considered domain is a closed system in the sense that the boundaries of the fluid domain consist of two solid walls (inner and outer cylinders) and two periodic boundaries (top and bottom), which do not allow for any net flow of liquid towards the apparatus. This issue is solved by making a small circular hole where the reference pressure value is set and a homogeneous Neumann boundary condition is applied for the velocity field. The hole has a diameter 
$r_{out}/12$ and is placed halfway along the axial length (i.e. at 
$L_z/2$) on the external cylinder, thus allowing the liquid to enter the domain as the bubble dissolves (see figure 12).
Figure 12. Opening of the outer cylinder for the passage of liquid (section taken at 
$z=L_z/2$). This modification is necessary to ensure the continuity of mass when the volume of the gas fraction decreases.
Simulations are run in a non-dimensional form, where the selected reference quantities for the four units involved are 
$L_{ref} = D_b/10$, 
$\rho _{ref} = \rho _c$, 
$U_{ref} = U_{in}$ and 
$c_{ref} = \rho _d/M$. The simulations are run for a (physical) time interval 
$\Delta t = 0.12$ s. To facilitate the comparison between different TC Reynolds numbers, the time scales and mass transfer coefficients are presented in the following in the corresponding dimensional form.
4.2. Mesh sensitivity
A different physical process (i.e. the mass transfer at the interface) requires a new mesh sensitivity study to find a suitable grid for mesh-independent solutions. The selected configuration for this analysis is that with a steady rotor (i.e. 
$Re=0$), which consists of a bubble rising in a quiescent flow bounded by cylindrical walls. The advantage of this configuration is that the finest mesh resolution is only needed around the bubble (with an AMR technique) and it is therefore significantly cheaper to run compared with the cases with Taylor vortices. The selected fields for grid refinement (see the brief introduction to AMR in § 2.2 and the references therein) are the volume fraction, species concentration and velocity field, with a threshold tolerance of 0.01 (made non-dimensional with 
$c_d$ and 
$U_{in}$ for concentration and velocity, respectively). The resulting mesh has the maximum level of refinement near the interface and in the wake region (as well as around the cylindrical walls), thus providing a suitable grid to capture both mass transfer in thin concentration boundary layers and the dynamics of highly deformed bubbles. At this point, it is important to remind here that the requirements in terms of grid density for the mass transfer depend on the Péclet number and this can obviously be affected by the rotor speed. However, for the considered bubble size (
$D_b^{t=0} = 5$ mm), the bubble Reynolds number (
$Re_b$) is mainly determined by the rising velocity and, therefore, the 
$Pe$ number is weakly dependent on the rotor speed (see the results in the following). Three different mesh refinements are compared here and the list of cases for the grid sensitivity study is reported in table 9. Results in terms of volume dissolution rates for the three considered meshes are reported in figure 13. Mesh M1 over-predicts the volume ratio as a result of the under-resolved concentration boundary layer at the gas–liquid interface, whereas meshes M2 and M3 are indistinguishable until 
$t\approx 0.05$ s and produce similar results for 
$t<0.08$ s. As the bubble volume is further reduced, mesh M2 deviates from mesh M3 because not enough cells are distributed around the interface. This is a common issue for dissolving bubbles, since no mesh can be fine enough to capture the mass transfer until complete dissolution. However, given the mesh-independent solution obtained for a volume reduction of up to 
$80\,\%$ (i.e. 
$V_b(t)/V_b^{t=0} = 0.2$) and the cheaper computational cost compared with case M3, mesh M2 is selected for all the other cases presented in this section.
Table 9. List of cases for the grid convergence analysis of a dissolving bubble in a TC device with no rotation.
Figure 13. Grid convergence for a dissolving bubble in a TC device with no rotation. Plot of bubble volume ratio vs time.
4.3. Effect of inner cylinder rotating speed
The effect of the Reynolds number is investigated by comparing the cases with 
$Re=0$ (no rotation) and 
$Re=1000$, 3000 and 5000, where the TC flow regime moves from WVF to TTVF (see table 1). The complete list of cases presented in the rest of this section is summarised in table 10, along with four cases for the study of wake effects for two bubbles. Cases A–D represent a realistic configuration, where the motion of the bubble is determined by two major components: the gravitational acceleration and the transport induced by the carrier liquid (TC flow). Although the effects on the distribution of the dissolved species within the device are clearly dependent on the rotor speed (as is shown later), for the selected bubble dimension (
$D_b^{t=0} = {5}$ mm) the motion is mainly dominated by the gravitational force. Therefore, for a better understanding on the role of Taylor vortices on the mass transfer of bubbles, gravity has been neglected in cases E–G and the bubble motion is made completely dependent on the carrier flow. Results are first presented for the cases with gravity (cases A–D, § 4.3.1) and subsequently the removal of the buoyancy force is discussed in cases E–G (§ 4.3.2). A study into the flow scales that control mass transfer is proposed in § 4.3.3, whereas wake effects are investigated in § 4.3.4.
Table 10. List of cases for the study of dissolving bubbles in a TC device at different rotating speeds and gravitational accelerations.
4.3.1. Single bubble with gravity
Results for cases A–D in terms of volume ratio against time are reported in figure 14. As was anticipated before, the velocity magnitude of the bubble is basically determined by the rising component and the volume dissolution rate for these cases is not affected significantly by the rotation of the inner cylinder (minor differences are observed at the start and at the end of the simulation, where cases C and D dissolve slightly faster than cases A and B, coherently with the larger rotating speed). The plot of the volume ratio shows a linear trend until 
$V_b(t)/V_b^{t=0} \approx 0.4$ and, after that, the slope progressively decreases as the bubble dissolves; a similar behaviour was observed for the mass transfer of a rising bubble in a quiescent flow (figure 11).
Figure 14. Volume ratio vs time for a dissolving bubble in a TC device at different rotating speeds. For the selected configuration, gravity is dominant and the TC flow plays a marginal role in the dissolution rate.
The (time-dependent) Sherwood number is monitored during the simulation and results are plotted in figure 15. The plots of the Sherwood number show a similar profile until 
$t \approx {0.06}$ s, where the size of the bubble is larger and the buoyancy effects are more relevant. However, for 
$t > {0.06}$ s, two different patterns that characterise cases A and B and cases C and D, respectively, are clearly observable. In the higher rotating speed cases (
$Re=3000$ and 5000), the Sherwood number 
$Sh$ is enhanced by the turbulent TC flow structures that develop within the apparatus, whereas almost no difference is observed between the steady case (
$Re=0$) and the wavy vortex regime (
$Re=1000$). However, such differences occur when the bubble volume is already reduced significantly (
$V_b/V_b^{t=0} < 0.3$) and no relevant effects in terms of dissolution rates can be observed afterwards. Cases A and B show a local peak around 
$t\approx {0.08}$ s that is larger than the values of 
$Sh$ for cases C and D; as will be shown later, this effect is due to the corresponding rising speed of the bubble.
Figure 15. Sherwood number vs time for a dissolving bubble in a TC device at different rotating speeds. The Sherwood number 
$Sh$ is based on the diameter of the equivalent sphere (4.1).
When the rotating speed of the inner cylinder is increased, the magnitude of the main (azimuthal) velocity component of the carrier fluid grows and the motion of the bubble is affected accordingly. Figure 16 compares the trajectory of the bubble centre on different planes for cases A–D (it is recalled here that the axis of the cylinders is aligned with the 
$z$ direction). When no rotation is applied, the bubble rises following an almost perfect rectilinear trajectory (figure 16a,b). As the Reynolds number is increased, the liquid velocity (combined with gravity) induces an irregular motion of the bubble, which results in a net anticlockwise displacement on the 
$xy$ plane (coherently with the rotation of the rotor) and in a more developed trajectory in the azimuthal direction for the most turbulent case (figure 16c). The rising trajectory of a bubble is the result of the interaction among several factors, such as the vorticity shed into the liquid phase, the shape deformation and the topology of the carrier flow. In the present cases, a very complex interaction of the aforementioned parameters is observed, where the bubble experiences shear rates on both the azimuthal (
$r\theta$) and 
$rz$ planes as well as chaotic fluctuations for the turbulent cases. Studies that have investigated the rising trajectory of bubbles in simple (planar) shear flows (e.g. Ervin & Tryggvason Reference Ervin and Tryggvason1997; Hidman et al. Reference Hidman, Stróm, Sasic and Sardina2022) have shown that a change in the sign of the lift force (i.e. the component acting perpendicular to the main rising direction) occurs when the shear rate increases. The bubble considered in the present work (with 
$Ga = 1050.7$, 
$Bo = 3.4$) belongs to a region in the 
$Ga$–
$Bo$ plane where such a change is observed. The case with 
$Re=1000$ does not deviate significantly from a rectilinear trajectory for most of the simulated time (figure 16a,b), leading to the conclusion that the TC flow pattern is not strong enough to perturb the buoyancy-dominated dynamics. On the other hand, when the volume of the bubble decreases and buoyancy is less predominant, a clear effect of the carrier flow is observed. When the rotating speed is increased (
$Re=3000$ and 
$5000$), a deviation from the straight rising path is induced by the combination of azimuthal and TC vortical flow patterns. Interestingly, the lift direction is opposite (figure 16b), with the case 
$Re = 3000$ initially attracted towards the inner cylinder and subsequently towards the outer one, whereas the 
$Re=5000$ case starts to rise vertically in the 
$yz$ plane before moving towards the rotating wall. A similar mechanism to that observed for lift forces that change sign with increasing shear rates and interfacial deformation could be at work in these cases. However, it is important to recall here that the present configuration (3-D and fully turbulent flow with phase-change) differs significantly from planar shear flows investigated in the literature. The first instances of the bubble motion are also affected by the relative position between the bubble centre and the Taylor vortices when the bubble is released into the flow. As was shown in § 3.1, TC flow changes from a well-organised and steady pattern of alternating vortices for 
$Re=1000$ to a fully turbulent and time-dependent configuration for 
$Re=3000$ and 5000. The initial position of the bubble is always the same for all the cases considered in this work (we recall here that the initial axial location of the bubble centre is 
$z_b^{t=0} = r_{out}/3$). Based on the flow visualisations shown in figure 7, for 
$Re=3000$ the bubble is initially located within a vortex (closer to the lower part of the vortical cell), whereas for 
$Re=5000$, the centre lies at the boundary between two adjacent vortices. In the 
$Re=3000$ case, the bubble initially experiences a negative net radial velocity, which results in a displacement towards the inner wall. In the most turbulent case, the flow around the bubble at 
$t = 0$ s has a positive radial component, although its initial displacement occurs mainly along the azimuthal direction. This effect can be explained with the strong unsteadiness that characterises the chaotic flow at 
$Re = 5000$, where the vortices are not steady but move considerably inside the reactor.
Figure 16. Bubbles rising trajectories projected on the (a) 
$xz$, (b) 
$yz$ and (c) 
$xy$ planes at different rotating speeds. Bubbles are initialised at 
$x=0$, 
$y = -(3/2) r_{in}$ and 
$z = (2/3) r_{in}$.
In order to compare the effects of the TC flow features and buoyancy on the dynamics of the rising bubble, we define two different Froude numbers. The first definition takes into account only the main azimuthal component of the liquid medium and reads
\begin{equation} Fr^{\theta} = \frac{\langle u_\theta\rangle _{z \theta t}(r_b)}{\sqrt{g D_b}}, \end{equation}
where 
$r_b$ is the radial position of the bubble centre. The velocity is averaged along the axial and azimuthal direction as well as in time and is evaluated from the single-phase cases (§ 3.1). The effect of the vertical component of Taylor vortices (which can enhance or reduce the rising speed of the bubbles) is quantified by the following Froude number:
\begin{equation} Fr^{TV} = \frac{|u^{TV}(r_b)|}{\sqrt{g D_b}}, \end{equation}
where 
$u^{TV}$ is the characteristic (axial) velocity component of Taylor vortices from the undisturbed flow (Chouippe et al. Reference Chouippe, Climent, Legendre and Gabillet2014). Taylor cells are assumed to be squared (of the size of the reactor gap) and the velocity profile is varied linearly between the two walls. The profiles of 
$Fr^{\theta }$ and 
$Fr^{TV}$ against time for the configurations with 
$Re=1000$, 3000 and 5000 are reported in figures 17(a) and 17(b), respectively. The laminar case (
$Re = 1000$) shows almost no influence of TC features (for both the main azimuthal and axial Taylor vortices components), consistently with the observed dynamics of the rising bubble that resembles the case without rotation very closely (e.g. rectilinear trajectory and Sherwood profile). The azimuthal component of TC flows becomes stronger as 
$Re$ increases and 
$Fr^{\theta }$ follows accordingly for 
$Re=3000$ and 5000, with 
$Fr^{\theta } \approx 1$ for 
$Re=5000$ at the end of the simulation. The axial component of Taylor vortices (figure 17b) is more effective for the intermediate case 
$Re=3000$ rather than for the most turbulent one. In the first part of the simulation (
$t<0.04$ s), the 
$Re=5000$ case moves exclusively along the azimuthal direction, with almost no deviation on the 
$yz$ plane (figure 16b). Therefore, the bubble stays at the centre of one of the Taylor cells and experiences no significant axial velocity field. On the other hand, the 
$Re=3000$ case deviates almost immediately from the centre of the reactor gap and moves towards a region with non-negligible axial flow. For 
$t > 0.06$ s the case with 
$Re=3000$ (5000) approaches the outer (inner) walls and the corresponding 
$Fr^{TV}$ increases accordingly. The Froude number for 
$Re=3000$ is still larger than the corresponding one for 
$Re=5000$ because the former approaches the wall region faster. It is also important to recall here that the ratio 
$u^{TV}/U_{in}$ is not constant, but decreases for increasing Reynolds numbers (Chouippe et al. Reference Chouippe, Climent, Legendre and Gabillet2014). Both Froude numbers ((4.4)–(4.5)) are below one for all the selected cases, leading to the conclusion that the presence of TC flow features introduces perturbations into the system (e.g. trajectory), but the dynamics of the dissolving bubbles is mainly driven by buoyancy.
Figure 17. Froude number based on (a) the main TC azimuthal flow and (b) the axial Taylor vortex component vs time for a dissolving bubble in a TC device at different rotating speeds.
Although the volume dissolution rates are basically the same for cases A–D, such different trajectories provide some useful information for the operation of TC devices as chemical reactors. Indeed, when the gas extracted from the disperse phase is needed to perform a chemical reaction within the liquid phase, the more the distribution of the dissolved species is spread in a wide area the more likely is that the reagents react and produce the desired product. The case with 
$Re=5000$ results in a more extended trajectory compared with the other cases, which helps distribute the gas in a wider region within the reaction vessel and, eventually, promote reactions. The effects of the trajectory and Taylor vortices on the distribution of species are shown in figure 18, where the concentration contours of the gas released into the liquid on a generic 
$rz$ plane that intersects the bubble and the corresponding isosurfaces are compared for cases A–D. The figure clearly shows the increasing complexity of the wake region and species distribution as the rotor speed increases. When the rotor is steady (figure 18a), a symmetric isosurface develops around the bubble and inside the wake. As the rotor is accelerated, the topology of the isosurface becomes more distorted and, in the fully turbulent case at 
$Re=5000$ (figure 18d), the distribution of species becomes well-mixed within a wide region below the bubble. As explained earlier, this is the most desirable scenario for the enhancement of the yield of a chemical reaction when the dissolved gas is one of the reactant species. Therefore, it can be concluded that, although no major differences are observed in these cases for the dissolution rates, the promotion of turbulent (chaotic) Taylor vortices is a desirable feature for the enhancement of species mixing within the reactor and, eventually, the production of chemical compounds.
Figure 18. Contours of dissolved gas concentration on a 
$rz$ plane (left) and corresponding isosurfaces with 
$c_c = 0.1 \rho _d/M$ (right) in a TC device at (a) 
$Re = 0$, (b) 
$Re=1000$, (c) 
$Re=3000$ and (d) 
$Re = 5000$. The outer cylinder has been removed to improve the clarity of the figure. Snapshots taken at 
$t = 0.1$ s.
Many attempts have been made in the literature to provide formulae for the prediction of Sherwood numbers in rising bubbly flows and, although no formula can be generic enough to be independent of the specific flow configuration, most of the available correlations relate Sherwood with Reynolds or Péclet numbers in a proportionality law. To the best of the authors’ knowledge, no specific relationships have been investigated for the mass transfer of a single bubble in a TC flow at different rotating speeds (and TC flow regimes). Here, the correlation between Sherwood and Reynolds numbers is first investigated for cases A–D and the results are reported in figure 19, where the reference length used for 
$Re_b$ is the equivalent diameter of a sphere (as is done for 
$Sh$). In all the tested configurations, the plots of the Reynolds numbers exhibit a similar trend until 
$t \approx {0.07}$ s, where a maximum peak is observed. In the first part of the simulation, the buoyancy force makes the bubble less sensitive to the carrier flow, which explains why the plots have a similar shape. Interestingly, the magnitude of the maximum 
$Re$ is larger for the no-rotation (figure 19a) and 
$Re=1000$ (figure 19b) cases than for the high-speed configurations (
$Re=3000$ and 5000 in figure 19c,d respectively), meaning that the presence of turbulent Taylor vortices induces a strong downwards (liquid) motion that limits the upwards (rising) bubble velocity component as induced by gravity. This effect is significantly stronger as the strength of Taylor vortices increases and explains why the maximum observed peak of Reynolds number is larger in cases A and B than the fully turbulent cases C and D. For 
$t>{0.07}$ s, cases A and B have a similar trend with a strong fluctuating profile and an almost constant mean value, whereas cases C and D have weaker oscillations but an average decreasing value of 
$Re$ over time. Interestingly, the 
$Re = 3000$ case has a larger 
$Sh$ compared with the most-turbulent case for 
$0.08\ \textrm {s} < t < 0.095\ \textrm {s}$, meaning that the local TC pattern (i.e. the upwards and downwards velocity regions) has a stronger effect than the magnitude of the rotating speed, coherently with the results reported in figure 17. The plots of Sherwood numbers in figure 19 clearly show that 
$Sh$ and 
$Re$ are intrinsically related, since both profiles appear similar and the peaks occur approximately at the same time (with a small delay in the Sherwood plot) for all the tested configurations. Given this correlation, it is not surprising that cases A and B show a larger 
$Sh$ number than cases C and D at 
$t\approx {0.07}$ s, as was observed (but not explained) in figure 15.
Figure 19. Plots of 
$Sh$ and 
$Re$ numbers vs time for a dissolving bubble in a TC device at (a) 
$Re=0$, (b) 
$Re=1000$, (c) 
$Re=3000$ and (d) 
$Re=5000$. The similarity of the profiles suggests a functional relationship between 
$Sh$ and 
$Re$, as found for rising bubbles in (unbounded) quiescent flows.
Following the qualitative results presented in figure 19, a conceptually equivalent proportionality law between 
$Sh$ and 
$Re$ to those proposed in the literature for a rising bubble is expected to be valid also for cases A–D. Here the corresponding Sherwood profiles are compared against the theoretical formulae proposed by Oellrich et al. (Reference Oellrich, Schmidt-Traub and Brauer1973) for small bubbles,
\begin{equation} Sh = 2 + 0.651 \frac{Pe^{1.72}}{1 + Pe^{1.22}}, \quad \text{for } Re_b \rightarrow 0, Sc \rightarrow \infty, \end{equation}and for large bubbles,
\begin{equation} Sh = 2 + \frac{0.232 Pe^{1.72}}{1+0.205 Pe^{1.22}}, \quad \text{for } Re_b \rightarrow \infty, Sc \rightarrow 0. \end{equation}
Equations (4.6)–(4.7) provide two boundary curves for 
$Sh$ and are generally used to predict the mass transfer of a single rising bubble in a steady-state regime, i.e. when 
$Pe$ is time-independent (Deising, Bothe & Marschall Reference Deising, Bothe and Marschall2018). Oellrich et al. (Reference Oellrich, Schmidt-Traub and Brauer1973) showed that the Sherwood number of spherical bubbles rising at constant speed is a function of both 
$Pe$ and 
$Sc$ and approaches equation (4.6) (equation (4.7)) for small (large) Péclet numbers. The Schmidt number affects how quickly a rising bubble migrates from (4.6) to 4.7 as 
$Pe$ increases: the larger the Schmidt value, the later such transition occurs. It is finally noted that (4.7) approaches the well-known potential flow solution 
$Sh=( 2/\sqrt {{\rm \pi} } ) \sqrt {Pe}$, for 
$Pe \rightarrow \infty$ (Levich Reference Levich1962). For the considered application, the Péclet number (
$=Re_{b}Sc$) changes over time and formulae (4.6)–(4.7) are compared against the numerical results by replacing 
$Pe$ with 
$Pe(t)$ in figure 20. Since correlation formulae for 
$Sh$ are generally based on the surface of the equivalent sphere (
$A_{sphere}$), a correction factor (
$Sr$) is needed for the numerical results (which are based on the effective surface 
$A_\varSigma$) to compare against the theoretical equations:
\begin{equation} Sr = \frac{A_\varSigma}{A_{sphere}}. \end{equation}
Here 
$Sr$, which is always 
$\geq$1, is also known as a shape factor and provides a parameter for the estimation of the bubble deformation. As the bubble dissolves, the surface tension becomes more relevant (larger curvature) and the bubble approaches the spherical shape, i.e. 
$Sr \rightarrow 1$. The results reported in figure 20 show that the qualitative trend of the corrected Sherwood number (i.e. 
$Sh \times Sr$) is reproduced correctly by the theoretical formulae of Oellrich et al. (Reference Oellrich, Schmidt-Traub and Brauer1973), where the solution is closer to (4.7) in the first part of the simulation (where 
$Pe$ is larger due to the buoyancy-induced rising speed and larger size of the bubble) and progressively approaches equation (4.6) as the bubble dissolves (and 
$Pe$ decreases), coherently with the range of validity of these formulae. The trend of a decreasing Sherwood when the Reynolds number reduces (e.g. in the last part of the simulation, for 
$t > {0.08}$ s) is also reproduced correctly. Similar conclusions were obtained in Maes & Soulaine (Reference Maes and Soulaine2020) for a dissolving bubble rising in a quiescent flow and the present results confirm that volume change effects can be qualitatively taken into account by replacing the steady-state non-dimensional numbers with the corresponding time-dependent values in the appropriate correlation formulae.
Figure 20. Comparison of the corrected Sherwood number against the theoretical formulae proposed by Oellrich et al. (Reference Oellrich, Schmidt-Traub and Brauer1973) for (a) 
$Re=0$, (b) 
$Re=1000$, (c) 
$Re=3000$ and (d) 
$Re=5000$.
As shown in figure 20, (4.6)–(4.7) can be used as qualitative references for the expected Sherwood number of a rising bubble in a Taylor–Couette reactor. However, a quantitative accurate match between the present results and these correlations cannot be obtained, as the theoretical formulae are derived assuming a spherical shape of the bubbles and a rectilinear rising trajectory. For the analysed configurations, the combined effect of gravity, TC flow and phase-change induce strong deformations (
$Sr > 1$) in the bubble shape, which are compared in figure 21 for cases A–D, along with the corresponding shape factors. Bubbles are initialised as perfect spheres (i.e. 
$Sr^{t=0} = 1$) and, as soon as the buoyancy force makes the bubble rise, the interface assumes the typical dimple shape that can be observed at 
$t \approx {0.02}$ s. The shape factor increases accordingly until 
$t \approx {0.06}$ s for cases A and B (figure 21a,b) and 
$t \approx {0.055}$ s for cases C and D (figure 21c,d), where a local maximum peak is reached. The corresponding deformations are different between the first two cases (ellipsoidal shape) and the fully turbulent cases (reverse dimple); the relative shape factors also differ and are stronger for cases A and B (
$Sr \approx 1.65$) than for configurations C and D (
$Sr \approx 1.56$). After this peak, two different behaviours can be observed: for the no-rotation and 
$Re=1000$ cases, a second maximum peak is reached slightly after 
$t={0.08}$ s of approximately 
$Sr \approx 1.75$, where the bubbles approach a (less-pronounced) dimple shape, whereas for cases 
$Re=3000$ and 5000 the profiles do not have such a significant peak and irregular shapes can be observed. As the volume of the bubble decreases, the surface tension force becomes dominant and all the bubbles move towards a spherical shape (
$Sr \rightarrow 1$). The time evolution of the shapes represented in figure 21 suggests that the bubble interface experiences a very complex dynamics due to a combination of wobbling effects (initial Morton number 
$Mo = 3.2 \times 10^{-11}$; see Clift, Grace & Weber Reference Clift, Grace and Weber2005) and volume dissolution driven by mass transfer. Such irregular interfacial deformations result in the primary cause of the fluctuations that characterise the Sherwood plots for the steady case (
$Re = 0$), where the carrier flow is at rest and does not exhibit any other time-dependent feature. However, as the bubble reduces its size, the perturbations induced by the TC flow (particularly the toroidal vortices) result in a larger effect and change the bubble dynamics significantly (e.g. rising trajectory, bubble shape, 
$Sh$ and 
$Re$ fluctuations).
Figure 21. Shape factor and bubble shapes vs time for a dissolving bubble in a TC device at (a) 
$Re=0$, (b) 
$Re=1000$, (c) 
$Re=3000$ and (d) 
$Re=5000$.
4.3.2. Single bubble without gravity
The motion induced by the buoyancy force is the most relevant component for the configurations analysed so far, i.e. cases A–D. In this section, the focus is on cases E–G (see table 10), where the initial bubble size is kept the same (i.e. 
$D_b^{t=0} = {0.005}$ m) and gravity is neglected. The effect of the rotor speed is first investigated by comparing the bubble volume dissolution rates in figure 22. The bubble dissolves now significantly faster as the inner cylinder is accelerated, in contrast to the cases with gravity (see figure 14) where the dissolution rates were independent of the rotor speed. This is the expected behaviour, since the bubble velocity is now entirely given by the carrier liquid, whose main velocity component 
$(u_\theta )$ increases with the rotating speed of the apparatus.
Figure 22. Volume ratio vs time for a dissolving bubble in a TC device at different rotating speeds. Gravity is not taken into account.
The effect on the Sherwood number is shown in figure 23. As expected, 
$Sh$ increases as the rotor is accelerated and, after a transient regime where 
$Sh$ decreases whilst a concentration boundary layer develops around the bubble interface, the profiles approach a quasi-steady-state solution. Case G exhibits a constant value over time, whereas cases E and F have a slightly decreasing trend. Some qualitative differences between the low-Reynolds-number case (
$Re=1000$) and the fully turbulent cases (
$Re=3000$ and 
$5000$) can be observed in the plots of figure 23. The presence of unstable and chaotic Taylor vortices induce some fluctuations in the Sherwood profiles for the turbulent cases, whereas the well-organised and steady flow structures that develop in the WVF regime do not introduce analogous perturbations in case E.
Figure 23. Sherwood number vs time for a dissolving bubble in a TC device at different rotating speeds. Gravity is not taken into account.
The volume ratio and Sherwood number are integral parameters that are mainly affected in these cases by the main component of the TC flow and do not provide insights into the effects of the different TC regimes that characterise the apparatus at different Reynolds numbers. To look into the effects of Taylor vortices on the distribution of the dissolved species in the liquid phase, the contours of species concentration for cases E–G are compared in figure 24. The concentration for case E (figure 24a,b) appears uniform around the interface of the bubble and quite similar to the symmetric distribution that characterises a suspended bubble in a stagnant flow, meaning that the effect of Taylor vortices is marginal at 
$Re=1000$. On the other hand, in cases F (figure 24c,d) and G (figure 24e, f), the effect of the turbulent Taylor cells is clearly visible in the spatial distributions of species concentration, which now assume irregular and non-symmetric shapes around the bubble. The position of the bubble centre in the vertical plane can be tracked by looking at the wake left by the dissolution of species (figure 24b,d, f), and it can be observed that the bubble stays at a constant axial position for 
$Re=1000$, whereas in the turbulent cases (
$Re=3000$ and 
$5000$) it moves upwards, transported by the upward velocity induced by the vortices. These results confirm that, in case E, Taylor cells play a marginal role and the bubble behaves as a particle transported by the azimuthal velocity component, whereas for the TTVF regime (cases F and G) Taylor vortices actively contribute to the dynamics of the bubble and distribute the concentration of the dissolved species in a wider region around the interface, which is a desirable scenario for a good mixing of species. It is finally observed that the concentration patterns shown in figure 24 have a significantly different structure compared with the case of a rising bubble. Indeed, for rising bubbles, the concentration boundary layer is thinner on top (where advection counteracts the effect of diffusion) and becomes thicker towards the rear of the bubble. For the case of a bubble transported by a TC flow without gravity, the convective transport induced by the azimuthal velocity component has the same magnitude on both the top and bottom sides of the bubble and its effect is uniform around the interface (figure 24a,c,e), in contrast to the convective component induced by Taylor vortices, which acts on the radial–axial plane and depends on the bubble position and flow configuration.
Figure 24. Contours of species concentration and bubble interface in a TC device without gravity at (a,b) 
$Re=1000$, (c,d) 
$Re=3000$ and (e, f) 
$Re = 5000$. Top view (left) and side view (right). The outer cylinder has been removed to improve the clarity of the figure. Snapshots taken at 
$t = 0.1$ s.
Figure 24 also shows the shape of the bubbles, which appears almost spherical (
$Sr \approx 1$) for all the tested configurations. This happens because the shear rate induced by the TC flow is not strong enough to overcome the surface tension and induce significant deformations of the interface, in contrast to cases A–D where gravity was responsible for strong deviations from the spherical shape (see figure 21).
4.3.3. Mass transfer models
In order to gain insights into the underlying physics of the problem and discern the relevant flow scales that control mass transfer, the surface-renewal theory (Danckwerts Reference Danckwerts1951) is applied to the cases under investigation. The fundamental interphase mass transfer mechanism of the surface-renewal theory follows from the penetration model of Higbie (Reference Higbie1935) and assumes that the species-absorbing fluid next to the interface is continuously refreshed with new elements from the bulk liquid. The corresponding mass transfer coefficient is
\begin{equation} k_m \sim \sqrt{\frac{D_c}{\varTheta}}, \end{equation}
where 
$\varTheta$ is a characteristic residence time of a fluid element adjacent to the interface. The characteristic time 
$\varTheta$ is not known a priori and some assumptions regarding the scales controlling mass transfer are required. In the following, we present two approaches for the prediction of 
$\varTheta$.
The first approach is based on the assumption that mass transfer is driven by the macroscopic flow pattern, i.e. the combination of buoyancy and TC flow that transports the interface. At this point, a distinction between the cases with (cases A–D) and without (cases E–G) gravity is required. As shown in the analysis of the Froude number (figure 17) for the cases with gravity, the TC flow introduces small perturbations to the dynamics of the bubble and mass transfer is mainly affected by the rising speed induced by buoyancy. Under this circumstance, the relative velocity between the phases can be assumed equal to the bubble velocity 
$U_b$ and the residence time follows as
\begin{equation} \varTheta \propto \frac{D_b}{U_b}. \end{equation}
By replacing equation (4.10) into (4.9) and using the definition of bubble Reynolds 
$( Re_b = \rho _c U_b D_b / \mu _c )$ and Sherwood numbers, it follows that
\begin{equation} \frac{Sh}{\sqrt{Sc}} \propto \sqrt{Re_b}. \end{equation}
In (4.11), 
$Sh$ is a function of the solely bubble Reynolds and Schmidt numbers and corresponds to the well-known functional relationship 
$Sh \propto \sqrt {Pe}$. This is indeed the case for the configurations with gravity considered in the present work (figure 20) and it further confirms that mass transfer is controlled by buoyancy in those cases. In the following, the focus is on the cases without gravity in order to discern the relevant scales involved in the mass transfer process for configurations entirely driven by TC flows.
In cases E–G (i.e. without gravity) the bubble is subject to a shear rate in the azimuthal direction, which depends on the radial distance from the inner wall and increases with the TC Reynolds number. In contrast to the gravity-driven cases, the liquid (shear) flow moves with the bubble (i.e. the whole fluid domain is rotating). A relative motion still exists due to the varying velocity field induced by the shear rate, which results in a flow direction (relative to the bubble centre) towards increasing 
$\theta$ around the bubble side exposed to the inner cylinder (it is recalled here that the rotor rotates towards increasing 
$\theta$); the opposite scenario occurs for the side that faces the outer wall. Such relative motion can be observed in figure 24(a) (anticlockwise rotation), where the species distribution tends to move towards increasing 
$\theta$ faster than the centre of the bubble on the side that faces the inner wall, whereas the opposite trend is observed for the other side. However, the average shear-rate within the TC device is not particularly strong for the cases considered in this work, except for the two regions near the walls (see figure 5). Given the initial size of the bubbles modelled in the present work, the surface-renewal mechanism related to the macroscopic (shear) flow is not expected to be at work in the cases under consideration and is not discussed further.
The second approach is based on the small-eddy model of Lamont & Scott (Reference Lamont and Scott1970), where the smallest turbulent eddies are expected to control the exchange of mass at the interface. In this scenario, the characteristic turbulent length (
$l_K$) and velocity (
$u_K$) scales are computed as
\begin{equation} l_K = \left( \frac{\nu_c^3}{\epsilon} \right)^{1/4} \quad \text{and} \quad u_K = ( \nu_c \epsilon )^{1/4}, \end{equation}
where 
$\epsilon$ is the rate of turbulent dissipation. The turbulent time scale 
$t_K$ follows from the corresponding length and velocity quantities and the residence time is assumed to be (Theofanous, Houze & Brumfield Reference Theofanous, Houze and Brumfield1976; Herlina & Wissink Reference Herlina and Wissink2016)
\begin{equation} \varTheta \propto t_K = \frac{l_K}{u_K} = \sqrt{\frac{\nu_c}{\epsilon}}. \end{equation}Finally, the mass transfer coefficient can be formulated as (4.9):
\begin{equation} k_m \propto \sqrt{D_c} \left( \frac{\epsilon}{\nu_c} \right)^{1/4}. \end{equation}
The average turbulent dissipation rate can be obtained as a function of the TC Reynolds number and geometry when the flow statistics are stationary, i.e. when the inner/outer torques balance out (
$T^{in} = T^{out} = T$; see figure 3). Since the mechanical power applied to the internal cylinder must be dissipated by the fluid viscosity, the average dissipation rate follows as (Tokgoz et al. Reference Tokgoz, Elsinga, Delfos and Westerweel2012)
\begin{equation} \bar{\epsilon} = \frac{T \omega_{in}}{\rho_c V}, \end{equation}
where 
$V$ is the volume of liquid contained inside the reactor. By replacing the torque 
$T$ with the corresponding non-dimensional one (
$G$) and applying Wendt's formula to predict its value (3.2), 
$\bar {\epsilon }$ can be re-formulated as
\begin{equation} \bar{\epsilon} = \frac{G \nu_c^2 \omega_{in}}{{\rm \pi} ( r^2_{out} - r^2_{in} )}. \end{equation}
Finally, the average 
$\bar {\epsilon }$ is substituted in (4.14) and the prediction of the mass transfer coefficient of the small-eddy model follows as
\begin{equation} k_m \propto \sqrt{D_c} \left( \frac{G \nu_c \omega_{in}}{{\rm \pi} ( r^2_{out} - r^2_{in} )} \right)^{1/4}. \end{equation} The results of the small-eddy model are reported in figure 25. The analytical prediction of (4.14) shows an increasing trend of mass transfer coefficient for increasing Reynolds (coherently with the dissolution rates reported in figure 22) and shows a good agreement with the computed values of 
$k_m$ for cases 
$Re=3000$ and 5000. The 
$Re=1000$ case is reported here for reference, but it is not surprising that it is significantly off compared with the analytical model, since this configuration is laminar and no turbulent eddies can be at work in this case. The good agreement offered by the small-eddy model suggests that, for fully turbulent cases, mass transfer is controlled by the dissipative turbulent structures, as recently found for bubbles dissolving in homogeneous isotropic turbulence (Farsoiya et al. Reference Farsoiya, Magdelaine, Antkowiak, Popinet and Deike2023). Coherently with this result, the 
$Re = 5000$ case is independent of the bubble size and approaches a steady-state solution, whereas the case with 
$Re = 3000$ exhibits a quasi-steady solution with a slight decreasing trend over time.
Figure 25. Comparison of the small-eddy model (4.17) against the quasi-steady-state mass transfer coefficients of cases E–G. The proportionality coefficient is 0.51. Here 
$k_m$ values are averaged over time for 
$0.08\ \textrm {s} < t < 0.1\ \textrm {s}$.
It is finally recalled here that the initial position of the bubbles is always the same for cases A–G, i.e. the centre at 
$t=0$ is placed halfway between the inner and outer walls. For cases without gravity, where bubbles are entirely transported by the carrier liquid flow field, there might be a dependency on the initial position. This is expected to be particularly important for small bubbles that can be entirely trapped within the velocity boundary layer near the cylindrical walls. These regions show steep velocity gradients and fluctuations (figures 5 and 6) and the resulting mass transfer rate can be affected.
4.3.4. Wake effect
In this section, the interaction between two (identical) bubbles in a TC flow at different Reynolds number is investigated in terms of volume dissolution rates. The set-up of the apparatus is the same as the one presented in § 4.1 (i.e. 
$\eta =0.5$ and 
$Re=0$, 1000, 3000 and 5000), but two bubbles (referred to as 
$b1$ and 
$b2$) with diameter 
$D_{b1}^{t=0} = D_{b2}^{t=0} = (r_{out} - r_{in})/3 = {5}$ mm are initially placed at 
$z_{b1}^{t=0} = r_{out}/3$ and 
$z_{b2}^{t=0} = 7r_{out}/12$, respectively, and the same 
$(x,y)$ coordinates (i.e. the minimum distance between the interfaces is equal to 
$D_{b1,b2}/2$). Bubble 
$b1$ is placed at the same initial axial location as the single-bubble cases presented in §§ 4.3.1–4.3.2, so that a straightforward quantification of the wake effect induced by bubble 
$b2$ can be achieved by simply monitoring the evolution of volume over time. A summary of the cases for the wake effect study is reported in table 10 (cases H–K).
A qualitative representation is shown in figure 26, where the 3-D shape of the bubbles is plotted, along with contours of species concentration in the continuous phase on a 
$rz$ plane. In the case with no rotation of the inner cylinder (figure 26a), the solution is axisymmetric and the bubbles approach a similar shape, whilst rising along a vertical trajectory. In this scenario, the top bubble 
$(b2)$ rises in a clean environment (i.e. no concentration of species is distributed in the continuous phase around the upstream side of the interface). On the other hand, the bottom bubble 
$(b1)$ is affected by the wake of the top one, which modifies both the velocity and concentration fields around the interface. In particular, 
$b1$ rises in a contaminated environment, where the species released from 
$b2$ as it dissolves increases (locally) the saturation ratio of the liquid solvent. As a consequence, the difference in species concentration 
$(\Delta c)$ between the (saturated) interface and the surrounding liquid at the top side of bubble 
$b1$ (which drives the transport of molecules across the interface) decreases and a lower dissolution rate is expected, according to (2.9).
Figure 26. Two dissolving bubbles and contours of species concentration on a 
$rz$ plane at (a) 
$Re=0$, (b) 
$Re=1000$, (c) 
$Re=3000$ and (d) 
$Re=5000$. The outer cylinder has been removed to improve the clarity of the figure. Snapshots taken at 
$t = 0.057$ s.
As the rotating speed of the inner cylinder increases, the development of counter-rotating toroidal vortices induces non-null velocity components along the axial and radial directions that break the symmetry of the problem (figure 26b–d). The deviation from the symmetrical solution becomes larger as the Reynolds number of the apparatus increases from 
$Re=1000$ (figure 26b, where the bubbles follow slightly different trajectories but keep a similar shape) to fully turbulent (figure 26c,d), where the rising path and shapes of 
$b1$ and 
$b2$ are strongly decoupled. Therefore, for increasing 
$Re$, it is expected that the dissolution rate of the downstream bubble becomes less affected by the wake of the upstream one, as both bubbles follow different trajectories and rise in a clean environment.
The above qualitative observations are confirmed, from a quantitative point of view, in figure 27, where the time-evolving volumes of bubbles 
$b1$ and 
$b2$ are plotted and compared against the single-bubble case. For all the selected Reynolds numbers, the top bubble behaves in the same way as the corresponding single-bubble simulation, confirming that 
$b2$ rises in a clean environment and is not affected by the presence of the second bubble (nor by a different initial position within the reactor). The plot of the downstream bubble exhibits two different regimes. For 
$t < 0.04$ s, the volume changes with the same dissolution rate as 
$b2$ (and the single-bubble case), whereas for 
$t > 0.04$ s the bubble dissolves with a slower rate. In the first part of the simulation, the wake of the top bubble is not fully developed yet and the species released into the liquid (due to its dissolution) does not affect the concentration field around the downstream bubble. Therefore, in the first part of the simulation, both 
$b1$ and 
$b2$ rise in a clean environment and exchange moles at the same rate. For 
$t > 0.04$ s, the wake of 
$b2$ is extended sufficiently to interact with the mass transfer process of the downstream bubble. In particular, the local concentration gradient at the upstream side of the interface decreases, due to the increase in the bulk concentration in 
$b2$'s wake region, and the diffusive transfer of moles becomes slower accordingly. The effect of the apparatus Reynolds number on the dissolution rate of 
$b1$ can be easily understood by comparing figure 27(a,b) against figure 27(c,d). For no rotation of the inner cylinder or low rotating speeds, the wake effect is stronger as the downstream bubble clearly dissolves more slowly than the upstream one; however, for the fully turbulent cases, such difference is much less relevant. This is due to the chaotic motion induced by the strong vortices at 
$Re=3000$ and 5000 that decouple the trajectories of the two bubbles. Figure 26(c,d) shows that the part of the interface of bubble 
$b1$ that is affected by the upstream wake is minimal (compared with the other two cases), making the mass transfer dynamics almost identical to the case of a bubble rising in a clean environment. The decoupling of the trajectories is mainly due to the presence of strong TC vortices in the turbulent cases, since for the configurations with 
$Re = 0$ and 1000 the wake interaction has a significant effect on the dissolution rate of the bottom bubble. Bubbles are initially exposed to concurrent upwards and downwards velocity fields and experience (with different ratios) the influence of two adjacent TC cells. However, a different local flow topology between the two bubbles is not sufficient to decouple the trajectories (see, for example, the case with 
$Re = 1000$), which happens only when the flow induced by the vortices is strong enough to change the sign of the lift force acting on each bubble.
Figure 27. Volume ratio vs time for two dissolving bubbles in a TC device at (a) 
$Re=0$, (b) 
$Re=1000$, (c) 
$Re=3000$ and (d) 
$Re=5000$. The top bubble is not affected by the bottom bubble and is equivalent to the single-bubble case. The bottom bubble dissolves slower and the wake effect becomes less relevant as the rotating speed increases.
5. Conclusions
In this work we adopted our recent numerical framework (Gennari et al. Reference Gennari, Jefferson-Loveday and Pickering2022) to investigate bubble dissolution in a wide range of TC flows. The methodology is first validated for single-phase TC flows (with radius ratios varying between 
$0.5< \eta < 0.91$) at Reynolds numbers in the range 
$338 < Re < 5000$, where the main regimes (TVF, WVF and TTVF) are all reproduced and good agreement is observed against previous investigations.
Bubble dissolution in TC flows is first studied for a single bubble with a 5 mm diameter in a reactor with a radius ratio of 
$\eta =0.5$ and a gap size of 15 mm. For this configuration, the buoyancy force is predominant over the velocity induced by the rotating wall and the global dissolution rate is almost unaltered in the range 
$0< Re<5000$. However, the concentration of species released from the bubble is affected significantly by the TC regime, as a fully turbulent and chaotic flow distributes the dissolved species in a wider region and enhances the mixing within the reactor. A clear correlation between 
$Sh$ and 
$Re$ numbers is observed for all the modelled TC regimes where the bubble is rising. The theoretical predictions proposed by Oellrich et al. (Reference Oellrich, Schmidt-Traub and Brauer1973) for the 
$Sh$ number of spherical bubbles rising along a straight trajectory are modified by replacing the constant Péclet number with the corresponding time-dependent version and by introducing a correction factor to take into account shape deformations. The results show that large bubbles tend to agree with the predictions for 
$Re_b \rightarrow \infty$, whereas small bubbles are close to the expected behaviour for 
$Re_b \rightarrow 0$, even in the case of large interfacial deformations. The modelling of a bubble in absence of gravity provides useful information to quantify the effect of the different TC regimes for the cases in which the buoyancy force is marginal (e.g. small bubbles). In this specific case, volume dissolution occurs significantly faster for increasing rotating speeds, and all the simulated TC Reynolds numbers approach a quasi-steady-state solution.
The mass transfer mechanism is investigated by applying the surface-renewal theory of Danckwerts (Reference Danckwerts1951). The characteristic time and length scales of the macroscopic flow field control the mass transfer process for the cases where buoyancy is predominant. This theory produces the well-known functional relationship 
$Sh \propto \sqrt {Pe}$, which is consistent with the formulae of Oellrich et al. (Reference Oellrich, Schmidt-Traub and Brauer1973) for 
$Pe \rightarrow \infty$. For the cases without gravity, the small-eddy theory of Lamont & Scott (Reference Lamont and Scott1970) is combined with the surface-renewal approach and a simple analytical model for the prediction of the mass transfer coefficient is proposed. Our results show that the smallest turbulent scales control the exchange of mass between the phases in fully turbulent TC flows.
The wake effects are studied by placing two bubbles (aligned vertically) in the reactor. It is shown that the top bubble is unaffected by the presence of the second bubble and dissolves as in the single-bubble case. However, the bottom bubble rises into a contaminated flow and for null (
$Re=0$) or low rotating speeds (
$Re=1000$) the dissolution rate decreases significantly (at 
$t=0.1$ s, the bottom bubble has a volume 41 % and 52 % bigger than the top bubble for 
$Re=0$ and 
$1000$, respectively). On the other hand, for turbulent TC flows, the trajectories of the bubbles are decoupled and similar global mass transfer rates are observed.
Acknowledgements
We thank M. Poliakoff, D. Giddings and M. Magnini for helpful discussions.
Funding
We thank EPSRC grant EP/P013341/1 for financial support. The calculations were performed using the University of Nottingham High Performance Computing Facility and Sulis at HPC Midlands+, which was funded by EPSRC grant EP/T022108/1.
Declaration of interests
The authors report no conflict of interest.
Appendix A. Mesh sensitivity and grid selection for fully resolved simulations of TC flows
A mesh sensitivity study is carried out for one of the most-demanding cases in terms of mesh resolution (i.e. 
$\eta =0.5$ and 
$Re=5000$), where the flow regime is fully turbulent and strong velocity fluctuations are expected near the walls. The octree grid structure of Basilisk is used for the discretisation of the domain and two cylindrical regions with thickness 
$\Delta h_{in} = \Delta h_{out} = 0.05( r_{out} - r_{in} )$ are used to set different mesh refinements near the walls (see figure 28). Therefore, three different subdomains can be identified within the annulus, i.e. the inner, outer and bulk regions. Three meshes are tested for the selected configuration and the corresponding parameters are reported in table 11. Mesh M1 has a uniform resolution within the gap, whereas meshes M2 and M3 take advantage of the two refinement regions to increase the grid density near the cylindrical walls (M2 and M3 have the same resolution near the walls, but a different mesh density in the bulk region). Numerical modelling of TC flows requires that enough grid points are distributed within the gap between the cylinders, in order to capture the complex flow features that develop as the rotating speed is increased. Meshes M1 and M2 have a similar number of radial points (i.e. 
$N_r = 61$ and 
$N_r=73$ respectively), where 
$N_r$ is computed as 
$N_r = N_r^{b} + N_r^{in} + N_r^{out}$. However, the cost in terms of total number of cells for this marginal increment of resolution along the radial direction is significantly large (see table 11). This is a limitation of the octree Cartesian grid structure, where mesh stretching is not allowed, i.e. the aspect ratio of each cell is fixed to one. Results from the selected meshes are compared for the average azimuthal velocity 
$\langle u_\theta \rangle _{z \theta t}$ (where the operator 
$\langle \,\rangle _{z \theta t}$ refers to the average in time and along the axial 
$z$ and azimuthal 
$\theta$ directions) and for the corresponding fluctuating component,
\begin{equation} u'_\theta = u_{\theta} - \langle u_\theta\rangle _t, \end{equation}which can be averaged in time in the following way:
\begin{equation} \langle {u_\theta'}^2\rangle _t = \langle u_\theta^2\rangle _t - \langle u_\theta\rangle ^2_t. \end{equation}
The time interval used for the computation of the average and fluctuating quantities corresponds to five revolutions, i.e. 
$\Delta t = 5 t_{rev}$, where 
$t_{rev} = 2 {\rm \pi}r_{in}/U_{in}$. Results for 
$\langle u_\theta \rangle _{z \theta t}$ and 
$\sqrt {\langle {u_\theta '}^2\rangle _{z \theta t}}$ are plotted in figure 29(a,b), respectively, and compared against the numerical study of Chouippe et al. (Reference Chouippe, Climent, Legendre and Gabillet2014). The results reported in figure 29 show that the average radial profiles of the plotted quantities are not significantly affected by an increase in the mesh resolution. Mesh M1 tends to slightly over-predict the velocity fluctuations near the inner wall and the coarser resolution around the cylinders, combined with the embedded boundary method, results in a underestimation of the tangential velocity at the inner rotor; meshes M2 and M3 provide almost the same results. The grids are compared in terms of wall unit resolutions in table 12, where the average viscous length scales 
$\delta ^{*,in}$ and 
$\delta ^{*,out}$ at the inner and outer cylinders, respectively, are computed as
\begin{equation} \delta^*_{in,out} = \frac{\nu_c}{u^*_{in,out}}, \end{equation}
where the friction velocity 
$u^*$ is obtained from the shear stress 
$\tau _w$:
\begin{equation} u^*_{in,out} = \sqrt{\frac{|\tau_w^{in,out}|}{\rho_c}}. \end{equation}
The shear stress in (A4) is the average value on the cylinders and follows from the integral torque 
$T_w$:
\begin{equation} \tau_w^{in,out} = \frac{T_w^{in,out}}{2 {\rm \pi}r^2_{in,out} L_z}. \end{equation}
The values of 
$\varDelta _{r^+}^{in,out}$ reported in table 12 are computed with the average wall shear stress (A5) and, due to the Cartesian structure of the mesh, the non-dimensional quantities 
$\varDelta _{z^+}^{in,out}$ and 
$r_{in,out}\varDelta _{\theta ^+}$ are the same as 
$\varDelta _{r^+}^{in,out}$. Meshes M2 and M3 have the same refinement near the walls and both have at least four cells within the viscous sublayer region, i.e. 
$r^+ < 5$. Given the results reported in figure 29 and the requirements in terms of mesh resolution for DNS (i.e. 
$\varDelta _{r^+} < 1$), mesh M2 is selected as the reference grid for the modelling of TC flows; the grids used for the other configurations have similar characteristics and their details are reported in table 13. All the meshes have the first cell centre within the non-dimensional distance 
$\varDelta _{r^+}^{in,out} < 1$ from the walls and have at least three cells within the regions 
$r^+_{in,out} < 5$. Exceptions are the configurations with 
$\eta =0.73$, 
$Re=1014$ and 
$\eta =0.91$, 
$Re=5000$, where 
$\varDelta _{r^+}$ is slightly above one at the wall. In the last case (
$\eta =0.91$), this is due to the small gap within the cylinders, where the maximum number of cells is limited by the Cartesian topology of the grid and a further level of refinement would generate too many cells along the axial and azimuthal directions that cannot be handled with the available computational resources.
Figure 28. Mesh refinements within two cylindrical regions (dashed lines) around the inner and outer walls.
Table 11. Mesh sensitivity study for the configuration 
$\eta =0.5$ and 
$Re=5000$. Here 
$N_z$, 
$N_r$ and 
$N_\theta$ are the number of cells along the axial, radial and azimuthal directions, respectively. The superscripts 
$N^{b}$, 
$N^{in}$ and 
$N^{out}$ refer to the bulk, inner and outer regions within the domain (see figure 28).
Figure 29. Mesh sensitivity study for the configuration with 
$\eta = 0.5$ and 
$Re = 5000$. The radial profiles of (a) the average azimuthal velocity and (b) fluctuations are compared against the work of Chouippe et al. (Reference Chouippe, Climent, Legendre and Gabillet2014).
Table 12. Mesh characteristics in terms of wall units and number of cells in the viscous sublayer for the configuration 
$\eta =0.5$, 
$Re=5000$.
Table 13. Selected mesh characteristics for the single-phase TC cases.
Open access
