2.1 Non-dimensionalization

To be able to explicitly explore the effect of chemical reaction rate and the diffusion of the species ( $A,B,C$ ) on the dynamics, we employ for non-dimensionalization the characteristic length, time, velocity, solute concentration, viscosity and pressure, as

(2.3a-c ) $$\begin{eqnarray}\displaystyle \boldsymbol{x}=\frac{\widetilde{\boldsymbol{x}}}{\sqrt{Qt_{f}}},\quad t=\frac{\widetilde{t}}{t_{f}},\quad \boldsymbol{u}=\frac{\widetilde{\boldsymbol{u}}}{\sqrt{Q/t_{f}}}, & & \displaystyle\end{eqnarray}$$

(2.3d-f ) $$\begin{eqnarray}\displaystyle (a,b,c)=\frac{(\widetilde{a},\widetilde{b},\widetilde{c})}{a_{0}},\quad \unicode[STIX]{x1D707}=\frac{\widetilde{\unicode[STIX]{x1D707}}}{\unicode[STIX]{x1D707}^{\ast }},\quad p=\frac{\widetilde{p}}{Q\unicode[STIX]{x1D707}^{\ast }/\unicode[STIX]{x1D705}}. & & \displaystyle\end{eqnarray}$$

Here, $t_{f}$ is the final duration of the fluid injection. The characteristic viscosity $\unicode[STIX]{x1D707}^{\ast }$ is defined below. The resultant initial value problem written in dimensionless form satisfies coupled nonlinear partial differential equations:

(2.4a-c ) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{u}=0,\quad \boldsymbol{u}=-\frac{1}{\unicode[STIX]{x1D707}(c)}\unicode[STIX]{x1D735}p,\quad \unicode[STIX]{x1D707}(c)=\exp (R_{c}c), & & \displaystyle\end{eqnarray}$$

and

(2.4d ) $$\begin{eqnarray}\displaystyle \frac{\unicode[STIX]{x2202}a}{\unicode[STIX]{x2202}t}+(\boldsymbol{u}\boldsymbol{\cdot }\unicode[STIX]{x1D735})a=\frac{1}{Pe}\unicode[STIX]{x1D6FB}^{2}a-Da\,ab, & & \displaystyle\end{eqnarray}$$

(2.4e ) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\unicode[STIX]{x2202}b}{\unicode[STIX]{x2202}t}+(\boldsymbol{u}\boldsymbol{\cdot }\unicode[STIX]{x1D735})b=\frac{1}{Pe}\unicode[STIX]{x1D6FB}^{2}b-Da\,ab, & \displaystyle\end{eqnarray}$$

(2.4f ) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\unicode[STIX]{x2202}c}{\unicode[STIX]{x2202}t}+(\boldsymbol{u}\boldsymbol{\cdot }\unicode[STIX]{x1D735})c=\frac{1}{Pe}\unicode[STIX]{x1D6FB}^{2}c+Da\,ab, & \displaystyle\end{eqnarray}$$

associated with the initial conditions

(2.5a ) $$\begin{eqnarray}\displaystyle \boldsymbol{u}=(x[2\unicode[STIX]{x03C0}(x^{2}+y^{2})]^{-1},y[2\unicode[STIX]{x03C0}(x^{2}+y^{2})]^{-1}), & & \displaystyle\end{eqnarray}$$

and

(2.5b ) $$\begin{eqnarray}\displaystyle a(x,y,t=0)=\left\{\begin{array}{@{}ll@{}}1, & 0<x^{2}+y^{2}\leqslant d^{2}/4\\ 0, & \text{otherwise},\end{array}\right. & & \displaystyle\end{eqnarray}$$

(2.5c,d ) $$\begin{eqnarray}\displaystyle b(x,y,t=0)=1-a(x,y,t=0),\quad c(x,y,t=0)=0\quad \forall x,y, & & \displaystyle\end{eqnarray}$$

wherein $d$ is the dimensionless diameter of a circular region. The three dimensionless parameters appearing in equations (2.4a–c )–(2.4f ) are

(2.6a-c ) $$\begin{eqnarray}\displaystyle R_{c}=\ln \left(\frac{\unicode[STIX]{x1D707}_{C}}{\unicode[STIX]{x1D707}^{\ast }}\right),\quad Da=\frac{t_{f}}{1/(ka_{0})},\quad Pe=\frac{Q}{D}, & & \displaystyle\end{eqnarray}$$

where $\unicode[STIX]{x1D707}^{\ast }=\unicode[STIX]{x1D707}(a+b=1,c=0)$ , $\unicode[STIX]{x1D707}_{C}=\unicode[STIX]{x1D707}(a+b=0,c=1)$ . Although, in general, the diffusion coefficients of the three species $A,B,C$ are different and not necessarily constant, we assume these to be identical and equal to a constant $D$ . This simplification allows us to focus on the understanding of the influence of a simple chemical reaction on the fingering instability in an otherwise stable miscible displacement. The definition of the Péclet number ( $Pe$ ) is similar to that used in experiments by Nagatsu et al. (Reference Nagatsu, Matsuda, Kato and Tada2007), which is a consequence of the scaling used, allowing a direct comparison with experiments. While $Pe$ determines the relative importance of the advection to the diffusion of the solute, the Damköhler number, $Da$ , is defined as the ratio of the advective time scale to the reaction time scale. A higher $Da$ corresponds to a smaller reaction time (a faster reaction rate) and hence the product is formed fast. Thus for every fixed $Pe$ , varying the two dimensionless parameters $Da$ and $R_{c}$ , we explore the effect of chemical reaction on fingering dynamics. Equations (2.4a–c )–(2.4f ) are solved using a hybridization of a compact finite difference method and the pseudo-spectral method for the spatial derivative, while a fully explicit third-order Runge–Kutta method is used for the time integration (Chen et al. Reference Chen, Huang, Gadêlha and Miranda2008, and refs. therein). Similar numerical schemes have been quite effective with the non-reactive counterparts, where good quantitative agreement with the theoretical results are available in the literature. To name a few, Chen & Meiburg (Reference Chen and Meiburg1998) reported a good agreement with the growth rate on a miscible interface, while for immiscible fluids, the number of fingers triggered (Chen, Huang & Miranda Reference Chen, Huang and Miranda2011, Reference Chen, Huang and Miranda2014) was found to be consistent with the theoretical results. We briefly discuss our method for the reactive miscible fluids below.

2.2 Numerical method

Numerical computations are performed in a square box $\unicode[STIX]{x1D6FA}=[-L,L]\times [-L,L]$ . In the absence of any viscosity gradient, the radial flow is a potential flow with velocity vector, $\boldsymbol{u}_{pot}=(x[2\unicode[STIX]{x03C0}(x^{2}+y^{2})]^{-1},y[2\unicode[STIX]{x03C0}(x^{2}+y^{2})]^{-1})$ . However, this potential flow is perturbed as the chemical reaction induces a viscosity variation at the miscible interface between the two reactants. In order to capture this deviation of the velocity from the potential flow velocity $\boldsymbol{u}_{pot}$ , we decompose the velocity vector at any time into two parts,

(2.7) $$\begin{eqnarray}\displaystyle \boldsymbol{u}=\boldsymbol{u}_{pot}+\boldsymbol{u}_{rot}. & & \displaystyle\end{eqnarray}$$

The singularity in the potential flow at the origin is regularized by introducing a Gaussian source of core $\unicode[STIX]{x1D70E}$ as in Chen et al. (Reference Chen, Huang, Gadêlha and Miranda2008):

(2.8) $$\begin{eqnarray}\displaystyle \boldsymbol{u}_{pot}=\left(\frac{(1-\text{e}^{-(x^{2}+y^{2})/\unicode[STIX]{x1D70E}^{2}})x}{2\unicode[STIX]{x03C0}(x^{2}+y^{2})},\frac{(1-\text{e}^{-(x^{2}+y^{2})/\unicode[STIX]{x1D70E}^{2}})y}{2\unicode[STIX]{x03C0}(x^{2}+y^{2})}\right), & & \displaystyle\end{eqnarray}$$

here $\unicode[STIX]{x1D70E}\leqslant d/2$ , the radius of the circle initially occupied by reactant $A$ . We take $d=0.15$ in all the simulations. Introducing a streamfunction $\unicode[STIX]{x1D713}$ such that $\boldsymbol{u}_{rot}=(\unicode[STIX]{x2202}\unicode[STIX]{x1D713}/\unicode[STIX]{x2202}y,-\unicode[STIX]{x2202}\unicode[STIX]{x1D713}/\unicode[STIX]{x2202}x)$ , we write the Darcy’s law, equation (2.4a–c ), by eliminating the pressure, as

(2.9a,b ) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D6FB}^{2}\unicode[STIX]{x1D713}=-\unicode[STIX]{x1D714},\quad \unicode[STIX]{x1D714}=R_{c}[v(\unicode[STIX]{x2202}c/\unicode[STIX]{x2202}x)-u(\unicode[STIX]{x2202}c/\unicode[STIX]{x2202}y)]. & & \displaystyle\end{eqnarray}$$

We discretize the advection and diffusion terms in equations (2.4d )–(2.4f ) using fourth-order and sixth-order finite difference formulae (Lele Reference Lele1992) to obtain semi-discretized equations

(2.10a-c ) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x2202}_{t}a=f(a,b,c),\quad \unicode[STIX]{x2202}_{t}b=g(a,b,c),\quad \unicode[STIX]{x2202}_{t}c=h(a,b,c), & & \displaystyle\end{eqnarray}$$

where $f,g,h$ are the discretized versions of the advection–diffusion–reaction operators corresponding to their continuous form in equations (2.4d )–(2.4f ). No-flux boundary conditions are used for the solute concentrations and $\unicode[STIX]{x1D713}=0$ across all four boundaries. The equations (2.10) are associated with the initial conditions (2.5a )–(2.5c,d ). For numerical convergence, the step-like initial condition is regularized as follows. Equation (2.5b ) tells us that the reactant $A$ is initially contained within a circle of radius $d/2$ , which is equivalent to a dimensionless diffusive time $t=d^{2}/8$ . Thus, the axisymmetric base state solution of Tan & Homsy (Reference Tan and Homsy1987) at $t=d^{2}/8$ is used as an initial condition for $a$ . We introduce a random disturbance of amplitude $\unicode[STIX]{x1D709}\sim O(10^{-2})$ around the points where $a=0.5$ , so that the concentration at those points is given as $a=0.5+\unicode[STIX]{x1D709}\sin (2\unicode[STIX]{x03C0}N),$ where $N$ is a random number between $0$ and $1$ . The same set of random numbers is used for all the simulations. The initial condition for $a$ subtracted from unity gives the initial condition for $b$ .

The system of equations (2.10) is solved using a fully explicit third-order Runge–Kutta method with adaptive time stepping that satisfies the Courant–Friedrichs–Lewy condition. Sixth-order compact finite difference approximations are used in equation (2.9) to compute $\unicode[STIX]{x1D714}$ . A fast and efficient solution of the Poisson equation in (2.9) is achieved using a hybridization of the compact finite difference and the pseudo-spectral methods. We employ a Galerkin-type discretization in the $x$ -direction and a sixth-order compact finite difference in the $y$ -direction for $\unicode[STIX]{x1D713}$ . The simulations are stopped as soon as the concentration near the boundary becomes more than $0.1$ for $a$ or less than $0.9$ for $b$ . In this study, we take $\unicode[STIX]{x1D6FA}=[-0.75,0.75]\times [-0.75,0.75]$ with $512\times 512$ grid points.

3.1 Quantitative analysis

For the quantitative analysis, we make the transformations from $(x,y)$ coordinates to $(r,\unicode[STIX]{x1D703})$ coordinates. The conversions are made from $r=\unicode[STIX]{x1D70E}$ to $r=L$ in order to avoid the effect of a Gaussian source, if any, near the origin. So the domain in the $(r,\unicode[STIX]{x1D703})$ coordinates is $[\unicode[STIX]{x1D70E},L]\times [0,2\unicode[STIX]{x03C0})$ . We carry the quantitative analysis in the $(r,\unicode[STIX]{x1D703})$ plane in the next two subsections §§ 3.1.1 and 3.1.2.

3.1.1 Onset of instability

Interfacial length is the length of the product ring for a stable displacement and the length of the curve along the fingering pattern for an unstable displacement. In order to quantify the interfacial length, we define it in the $(r,\unicode[STIX]{x1D703})$ plane as,

(3.1) $$\begin{eqnarray}\displaystyle I(t)=\int _{\unicode[STIX]{x1D70E}}^{L}\int _{0}^{2\unicode[STIX]{x03C0}}\sqrt{\left(\frac{\unicode[STIX]{x2202}c}{\unicode[STIX]{x2202}r}\right)^{2}+\frac{1}{r^{2}}\left(\frac{\unicode[STIX]{x2202}c}{\unicode[STIX]{x2202}\unicode[STIX]{x1D703}}\right)^{2}}r\,\text{d}\unicode[STIX]{x1D703}\,\text{d}r. & & \displaystyle\end{eqnarray}$$

The numerical integration is carried out using the trapezoidal rule. In the $(r,\unicode[STIX]{x1D703})$ plane, $\text{d}r$ is fixed by the spatial discretization in the $(x,y)$ plane as $\text{d}r=\sqrt{\text{d}x^{2}+\text{d}y^{2}}$ but $\text{d}\unicode[STIX]{x1D703}=\tan ^{-1}(\text{d}y/\text{d}x)$ results in a poor discretization. So, we fix $\text{d}\unicode[STIX]{x1D703}=10^{-4}$ and carry out the conversions from the $(x,y)$ to the $(r,\unicode[STIX]{x1D703})$ plane. It is verified that decreasing $\text{d}\unicode[STIX]{x1D703}$ further does not qualitatively or quantitatively change our results. $I(t)$ normalized with $d\unicode[STIX]{x1D703}$ is shown in figure 5.

Figure 5.

Interfacial length as a function of time for $Da=100,Pe=3000$ and various $R_{c}$ . Inset shows zoomed plot around time $t\in (0.04,0.44)$ in order to highlight the onset time for $R_{c}=3.5$ and $R_{c}=-3.5$ . Clearly, the instability is triggered late for $R_{c}=-3.5$ .

$I(t)$ can be used as a tool to gain an insight into the onset of instability. The instant, when $I(t)$ for $R_{c}\neq 0$ departs from that for $R_{c}=0$ , is the time when instability triggers and is termed the onset time. It is observable from figure 5 that the difference between the onset time for $R_{c}<0$ and the corresponding positive $R_{c}$ decreases as we increase $|R_{c}|$ and ultimately the onset times become equal. This is because, before fingering occurs, the product is formed in a ring and the difference in onset time is attributed to the fluid velocity at the unstable interface of the ring of product. For $R_{c}<0$ , fingers are formed at the outer interface of the ring, where the fluid velocity of the unperturbed flow is smaller than the velocity at the inner interface of the product ring. Therefore the advective force, which is responsible for the fingering instability, is weaker in the former case. As $|R_{c}|$ increases, the onset time of instability decreases and it is independent of width of the ring. Hence for a given $Da$ , onset time for $R_{c}<0$ is always less than, or equal to, that for the corresponding positive $R_{c}$ . This is in contradiction to the linear stability analysis in a rectilinear geometry (Hejazi et al. Reference Hejazi, Trevelyan, Azaiez and De Wit2010) and can be attributed to the spatially varying velocity of the potential flow. Thus the geometry used for the nonlinear simulations does play a role in the dynamics and can lead to different results. Also, $I(t)$ for many non-zero $R_{c}$ coincides with that for $R_{c}=0$ for all times. This restates the existence of stable flows for a range of $R_{c}$ despite an unfavourable viscosity contrast.

3.1.2 Effect of barrier

It is evident from figure 2 that the non-monotonicity introduced in the viscosity profile by the chemical reaction results in a hydrodynamically stable interface where a more viscous fluid displaces a less viscous one. The stable interface hinders the fingers growing at the other interface. This is observable from the density plots in figure 4 where the fingers either only grow inwards for $R_{c}>0$ or outwards for $R_{c}<0$ . In order to quantify the effect of the barrier, we define the angular average concentration as

(3.2) $$\begin{eqnarray}\displaystyle \langle c(r,t)\rangle =(2\unicode[STIX]{x03C0})^{-1}\int _{0}^{2\unicode[STIX]{x03C0}}c(r,\unicode[STIX]{x1D703},t)\,\text{d}\unicode[STIX]{x1D703}. & & \displaystyle\end{eqnarray}$$

This quantification is possible only due to the transformation from the $(x,y)$ to $(r,\unicode[STIX]{x1D703})$ plane. Figure 6 (a) shows $\langle c(r,t)\rangle$ at $t=0.9$ for various $R_{c}$ and $Da=100,Pe=3000$ . The bumps towards the left for $R_{c}=7$ and towards the right for $R_{c}=-7$ show that fingers predominantly grow inwards and outwards, respectively for positive and negative $R_{c}$ . Further, we use $\langle c(r,t)\rangle$ to find the length of the fingers. We define the fingers growing at the outer interface of the product ring as outward fingers while the fingers at the inner interface are defined as inward fingers. To analyse the effect of the barrier induced by the chemical reaction, we define the outward and the inward finger lengths as

(3.3a,b ) $$\begin{eqnarray}\displaystyle L^{+}=|r_{0,out}-r_{out}|\quad \text{and}\quad L^{-}=|r_{0,in}-r_{in}|, & & \displaystyle\end{eqnarray}$$

respectively. Here, $r_{0,in}$ and $r_{0,out}$ are the radii when $\langle c(r,t)\rangle \approx 10^{-4}$ for $R_{c}=0$ as shown in figure 6(a) and $r_{in}$ and $r_{out}$ are the corresponding values for $R_{c}\neq 0$ . Figure 6(b) shows $L^{+}$ is larger (smaller) than $L^{-}$ for negative (positive) $R_{c}$ , thus capturing the effect of the barrier introduced by the non-monotonicity, on the finger length.

Figure 6.

(a) Angular average concentration, $\langle c(r,t)\rangle$ at $t=0.9$ for $Da=100,Pe=3000$ and various $R_{c}$ . $r_{0,in}$ , $r_{0,out},r_{in}$ and $r_{out}$ used to compute $L^{+}$ and $L^{-}$ are also shown. (b) Outward and inward finger lengths respectively for $R_{c}=-7$ and 7 for $Da=20$ .

Figure 7.

Outward and inward finger lengths respectively for $R_{c}=-7$ and 7 for $Da=100$ , 20 and 10.

We capture the effect of the chemical reaction on the finger length in figure 7. It is clear that a reaction producing a more viscous product is more unstable as the length of the fingers for positive $R_{c}$ is always more than that for corresponding $R_{c}<0$ , for fixed $Da$ . This is due to the radially varying potential component of velocity which delays the onset and consequently reduces the finger length for $R_{c}<0$ . Further, as $Da$ increases, the two finger lengths are found to approach each other, with $L^{-}$ ultimately surpassing $L^{+}$ . This is attributed to the barrier in the direction of bulk fluid motion for positive $R_{c}$ which hinders the growth of the fingers.

3.2 Critical $R_{c}$ for the fingering instability

The above discussions reveal that, for each $Da$ , a minimum viscosity contrast is required for the fingering instability to occur. To determine the critical parameters for the fingering instability, we use the interfacial length corresponding to $R_{c}=0$ as the reference. A set of parameters ( $Da,R_{c}$ ) is assigned an unstable displacement if the corresponding interfacial length exceeds the reference value by more than $5\,\%$ ; otherwise it is assigned as a stable displacement. Thus obtained stable and unstable cases are summarized in figure 8 for $Pe=3000$ , which depicts that the $Da-R_{c}$ plane consists of a stable region around $R_{c}=0$ , sandwiched between two unstable regions. This stable region exists for all $Da$ values considered and narrows as $Da$ increases. Therefore, for every finite $Da$ , we capture critical values of $R_{c}$ denoted as  $R_{c}^{-}$  and  $R_{c}^{+}$ , such that no VF occurs when $R_{c}^{-}\leqslant R_{c}\leqslant R_{c}^{+}$ and VF is observed outside this interval. Symbols in figure 8 denote the critical parameters $R_{c}^{-}$ (squares), $R_{c}^{+}$ (circles), while the lines represent a least squares fit to the function $|R_{c}^{\pm }|\propto Da^{-\unicode[STIX]{x1D6FD}_{\pm }}$ . A naive observation to be made here is that a higher viscosity gradient is required to trigger the instability for the reaction producing a less viscous product as $|R_{c}^{-}|>|R_{c}^{+}|$ for all $Da$ . The numerical simulations are carried out for a wide range of $Pe$ and the qualitative results presented in figure 8 remain unaffected for other values of $Pe$ as well, as is evident from the $Da-R_{c}$ parameter space for $Pe=1000$ in figure 9. It is verified that the stable region reduces as $Pe$ increases, but never completely disappears.

Figure 8.

$Da-R_{c}$ parameter space showing stable as well as unstable regions for $Pe=3000$ . ▫ denote $R_{c}^{-}$ and ○ denote $R_{c}^{+}$ , while the solid lines represent the power laws $R_{c}^{+}=8.8318Da^{-0.3365}$ and $R_{c}^{-}=-11.2329Da^{-0.3833}$ . The truncated density plots of product concentration in $(r,\unicode[STIX]{x1D703})$ coordinates representing each zone are also shown in the inset.

Figure 9.

Widening of the stable region on decreasing $Pe$ . $Da-R_{c}$ parameter space for $Pe=1000$ with the solid lines representing the power laws $|R_{c}^{\pm }|=\unicode[STIX]{x1D6FC}_{\pm }Da^{-\unicode[STIX]{x1D6FD}_{\pm }}$ . Critical $R_{c}$ is shown by symbols with ▫ and ○ respectively denoting $R_{c}^{-}$ and $R_{c}^{+}$ .

Now, we give a relation between the critical parameters so as to mathematically determine the bounds of the stable region for any given moderate $Da$ . Bischofberger, Ramachandran & Nagel (Reference Bischofberger, Ramachandran and Nagel2014) experimentally determined the existence of a critical viscosity contrast ( $R_{crit}$ ) for the fingering instability to occur in the miscible VF. Following the stability analysis of the fingering instability, we can write the critical viscosity ratio as $\unicode[STIX]{x1D707}_{crit}\approx \text{e}^{R_{crit}}$ . We extend this idea for the current case of reactive fingering instability in miscible fluids. We write the critical viscosity $\unicode[STIX]{x1D707}_{crit}^{\pm }$ as $\text{e}^{|R_{c}^{\pm }|f(Da)}$ , where $f(Da)$ is a function of $Da$ , to be determined. Setting $\unicode[STIX]{x1D707}_{crit}^{\pm }=\unicode[STIX]{x1D707}_{crit}$ and using a best fit of the numerically obtained critical parameter, we obtain

(3.4a,b ) $$\begin{eqnarray}\displaystyle R_{c}^{+}=\unicode[STIX]{x1D6FC}_{+}Da_{c}^{-\unicode[STIX]{x1D6FD}_{+}},\quad R_{c}^{-}=-\unicode[STIX]{x1D6FC}_{-}Da_{c}^{-\unicode[STIX]{x1D6FD}_{-}}, & & \displaystyle\end{eqnarray}$$

where $\unicode[STIX]{x1D6FC}_{\pm }$ and $\unicode[STIX]{x1D6FD}_{\pm }$ are the fitting parameters which depend on $Pe$ . These parameters for $Pe=3000,1000$ are given in figures 8 and 9, respectively. Clearly the stable region widens on decreasing $Pe$ . This clearly is a consequence of the effect of the stabilizing diffusive force on miscible VF. Thus, for a given $Pe$ , we can divide the  $Da-R_{c}$ parameter space into three regions which can be used for suitably choosing the parameters to control VF.