TY - JOUR
AU - Daripa, Prabir
AB - Abstract We perform a linear stability analysis of three-layer radial porous media and Hele-Shaw flows with variable viscosity in the middle layer. A nonlinear change of variables results in an eigenvalue problem that has time-dependent coefficients and eigenvalue-dependent boundary conditions. We study this eigenvalue problem and find upper bounds on the spectrum. We also give a characterization of the eigenvalues and prescribe a measure for which the eigenfunctions form an orthonormal basis of the corresponding |$L^2$| space. This allows for any initial perturbation of the interfaces and viscosity profile to be easily expanded in terms of the eigenfunctions by using the inner product of the |$L^2$| space, thus providing an efficient method for simulating the growth of the perturbations via the linear theory. The limit as the viscosity gradient goes to zero is compared with previous results on multi-layer radial flows. We then numerically compute the eigenvalues and obtain, among other results, optimal profiles within certain classes of functions. 1. Introduction Saffman–Taylor instability occurs when a less viscous fluid drives a more viscous fluid in a porous medium. This phenomenon occurs in many applications including carbon sequestration, filtration, hydrology and petrology. One important application is oil recovery in which it is common practice to use water to displace oil. Oil recovery was the application driving the seminal work on this type of instability in Saffman & Taylor (1958). In order to simplify their experiments, Saffman and Taylor studied the instability in the context of Hele-Shaw flows—flows between two parallel plates with a small gap between them. Hele-Shaw flows are a good model of porous media flows because they are also governed by Darcy’s Law. It is now well known that a positive interfacial viscosity jump in the direction of rectilinear flow produces an unstable flow with interfacial tension stabilizing short waves. In applications in which it is advantageous to suppress the instability, one strategy is to use intermediate fluids that are more viscous than the displacing fluid but less viscous than the displaced fluid in order to limit the size of the viscosity jumps. This strategy is effective when the interfacial tension at the interfaces is comparable to that of a two-layer flow. This has been demonstrated in Daripa (2008) in the case of rectilinear flow using upper bounds on the growth rates, which were also used to give stabilization criteria. Another approach to limiting the instability is using a variable viscosity fluid as the intermediate fluid. This can be achieved, for example, in chemical enhanced oil recovery (EOR) in which a polymer is used to increase the viscosity of the displacing fluid and variable polymer concentration leads to variable viscosity. A viscosity gradient allows for even smaller viscosity jumps at the interfaces than a constant viscosity fluid but the layer itself becomes unstable. There are several studies on the stability of multi-layer variable viscosity porous media and Hele-Shaw flows in a rectilinear geometry. Gorell & Homsy (1983) first theoretically studied the stability of three-layer flows in such a geometry in which the middle fluid has variable viscosity. However, they studied the restricted case in which the trailing interface is a miscible interface with no interfacial tension. Daripa & Pasa (2006) dropped this restriction and studied the case of three-layer variable viscosity flow with both interfaces immiscible. They used the variational form of the problem to derive upper bounds on the growth rate of instabilities and demonstrated that the growth rate can be decreased in comparison to the two-layer layer flow that does not have the intermediate fluid. Saffman–Taylor instability is also studied in a radial flow geometry in which there is a point source in the center of the flow and the fluid moves outward radially with circular interfaces between fluids. In the context of oil recovery, the radial geometry models the flow near the injection and production well whereas the rectilinear geometry models the flow far from the wells. Modeling near the well is more difficult because the well acts as an approximate singularity (source or sink depending on whether the well is a production well or injection well). This paper advances understanding of such near-well flows and has direct bearing on oil recovery technologies. Moreover, due to the rich and interesting dynamics of the interfaces in radial flows, which are still not well understood from a physical standpoint, this topic has been and will remain of interest in the near future. Radial flow is one of several cases investigated in Muskat (1934, 1946) for the two-layer problem but in the case of zero interfacial tension. Paterson (1981) later performed a linear stability analysis for two-layer radial flow with interfacial tension. There are relatively few works on the stability of ‘multi-layer’ radial flows. Cardoso & Woods (1995) studied the stability of three-layer radial flows in the limiting case in which the inner interface is completely stable and looked at the break up of the middle layer into drops. Beeson-Jones & Woods (2015) analyzed general three-layer flows, and Gin & Daripa (2015) performed a linear stability analysis for an arbitrary number of fluid layers in the radial geometry. In each of these works, the viscosity of the fluids increases in the direction of the flow, but the viscosity in each layer is constant. To date, there are no known studies of multi-layer radial flows with ‘variable’ viscosity. The development of the theory for radial flow lags behind that of rectilinear flow because of the challenges due to the time-dependence of the problem. In particular, the curvature of the interfaces, the length of the middle layer(s) of fluid and the spatially dependent viscosity profile are all time-dependent. Therefore, the analysis of rectilinear flow in, for example, Daripa & Pasa (2006) does not hold and a novel approach is required. Previous stability studies of radial Hele-Shaw flows (Cardoso & Woods, 1995; Paterson, 1981; Kim et al., 2009; Dallaston & McCue, 2013; Anjos et al., 2015; Beeson-Jones & Woods, 2015, 2017) allow for a time-dependent growth rate of instabilities, typically making use of a quasi-steady-state approximation. Despite this approximation, the stability results have often shown agreement with both physical and numerical experiments. In what follows, we keep with this tradition and use a quasi-steady-state approximation to derive an eigenvalue problem that has time-dependent coefficients and thus time-dependent eigenvalues. Additionally, the eigenvalue problem has eigenvalue-dependent boundary conditions, which make the study of the eigenvalue problem mathematically interesting and challenging. The eigenvalues depend on many different parameters including the viscosity profiles of the fluid layers, the interfacial tension at each interface, the curvature of the interfaces and the fluid injection rate. Therefore, a numerical exploration of the vast parameter space is necessary. It is also important for applications in chemical EOR to investigate which viscosity profiles minimize the instability for a given set of parameters, as is done in Daripa & Ding (2012) in the rectilinear case. There is also considerable interest in studying viscous fingering in cases where certain parameters of the flow or geometry are time-dependent. Li et al. (2009) numerically and experimentally studied radial Hele-Shaw flows with a time-dependent injection rate. Dias & Miranda (2010) used a weakly nonlinear approach to analytically study both a time-dependent injection rate and a time-dependent gap width of the Hele-Shaw cell. Zheng et al. (2015) showed through both linear stability analysis and experiments how to control viscous fingering for a radial Hele-Shaw flow using a multitude of time-dependent strategies with a particular emphasis on a time-dependent gap width. The results of Zheng et al. (2015) were recently extended with a more detailed stability analysis in Vaquero-Stainer et al. (2019). Morrow et al. (2019) used numerical simulations to compare constant versus time-dependent injection rates in standard, tapered and rotating Hele-Shaw cells. All of these studies were for the case of a single fluid interface. There are several physical mechanisms through which a variable viscosity profile can arise in the context of porous media flows. The simplest model to study incompressible porous media flow is to use the Hele-Shaw model which is based on the incompressibility condition and the Darcy’s law in each of the phases. This model does not allow any mixing between the oil and water phases macroscopically and maintains sharp interfaces between the phases. Another model, called the Buckley–Leverett model, builds on this Hele-Shaw model by adding a saturation equation, a nonlinear hyperbolic conservation law with non-convex flux function which allows mixing between the phases (macroscopically) due to rarefaction waves behind the leading saturation front sweeping the oil ahead (Daripa et al., 1988). The rarefaction waves create a viscosity profile behind the front with viscosity gradually increasing towards the moving front. This region with a viscosity profile is usually finite in length, and its length grows with time. The study of the stability of such composite solutions to conservation laws is relevant for porous media but less well-developed. An appropriate model for this problem is the one under study which models the effect of rarefaction waves with a viscosity profile between two interfaces and the shock front with a material interface having appropriate viscosity jump at the interface. Such a study was partially carried out in Daripa (2008) in the rectilinear geometry but the present study in the radial geometry is much harder as mentioned above and is also more relevant for porous media flow. A variable viscosity profile can also arise in a Hele-Shaw flow when there are particles in the displacing fluid. Due to shear-induced migration, the particles accumulate near the leading interface. The non-uniform concentration profile leads to a viscosity profile with the viscosity increasing in the direction of the flow. It has been shown experimentally that an otherwise stable radial Hele-Shaw flow can experience viscous fingering due to the viscosity gradient (Tang et al., 2000; Kim et al., 2017; Luo et al., 2018; Xu & Lee, 2019). The concentration profile, and hence the effective viscosity profile, are often measured as a function of the radial distance from the injection point in the Hele-Shaw cell (Luo et al., 2018; Xu & Lee, 2019). The paper is laid out as follows. In Section 2 we perform a linear stability analysis of a point source driven three-layer ‘radial’ Hele-Shaw flow in which fluid between the two interfaces has a variable viscosity profile. We use a time-dependent coordinate transformation to freeze the basic motion of the two interfaces and derive the associated eigenvalue problem in this new coordinate system. The growth rate of disturbances in the transformed coordinate system is related to the physical growth of disturbances of the interfaces in Section 3. Section 4 gives the restriction of the problem to the case of constant viscosity. In Section 5, we follow the approach of Daripa & Pasa (2006) to derive upper bounds on the growth rate using the variational form of the problem. The nature of the eigenvalues and the completeness of the eigenfunctions are investigated in Section 6. Numerical evaluation of the eigenvalues and their dependence on certain physical parameters are given in Section 7, and the validity of the quasi-steady-state approximation is explored. We conclude in Section 8. 2. Preliminaries We consider a radial Hele-Shaw flow consisting of three regions of incompressible, immiscible fluid. By averaging across the gap, we may consider a two-dimensional flow domain in polar coordinates, |$\varOmega := (r,\theta ) = \mathbb{R}^{2}$|⁠. The least viscous fluid with constant viscosity |$\mu _i$| is injected into the center of the cell at an injection rate, |$Q$|⁠. The most viscous fluid, with constant viscosity |$\mu _o$|⁠, is the outermost fluid. The middle layer fluid has a smooth, axisymmetric viscosity profile |$\mu (r)$| where |$\mu _i < \mu (r) < \mu _o$|⁠. The fluid flow is governed by the following equations $$\begin{align} &\boldsymbol{\nabla\cdot}{\textbf{u}} = 0,\qquad \boldsymbol{\nabla}\;p=-\mu\;{\textbf{u}}, \qquad \frac{\partial \mu}{\partial t} + {\textbf{u}} \cdot \nabla \mu = 0, \qquad \textrm{for} r \neq 0. \end{align}$$(2.1) The first equation (2.1) is the continuity equation for incompressible flow, the second equation (2.1) is Darcy’s Law and the third equation (2.1) is an advection equation for viscosity. Note that Darcy’s Law for Hele-Shaw flows contains a permeability term |$K = b^2/12$|⁠, but for the sake of simplicity here we have scaled viscosity by this term. Therefore, in what follows |$\mu $| denotes the modified viscosity. We start with the fluids separated by circular interfaces with radii |$R_1(0)$| and |$R_2(0)$|⁠, where |$R_1(t)$| and |$R_2(t)$| are the positions of the interfaces at time |$t$|⁠. This setup is shown in Fig. 1. Fig. 1. Open in new tabDownload slide The basic solution for three-layer flow. Fig. 1. Open in new tabDownload slide The basic solution for three-layer flow. The equations admit a simple basic solution in which all of the fluid moves outward radially with velocity |${\textbf{u}}:= \left (u_r,u_{\theta }\right ) = \left (Q/(2 \pi r),0\right )$|⁠. The interfaces remain circular and their radii are given by |$R_1(t) = \sqrt{Qt/\pi + R_1(0)^2}$| and |$R_2(t) = \sqrt{Qt/\pi + R_2(0)^2}$|⁠. The pressure, |$p_b = p_b(r)$|⁠, may be obtained by integrating equation (2.1). We define the quantity |$R_0(t) = \sqrt{Qt/\pi }$| and define the following coordinate transformation: $$\begin{align}& \zeta = \frac{r^2 - R_0^2(t)}{R_2^2(t) - R_0^2(t)} = \frac{r^2 - R_0^2(t)}{R_2^2(0)}, \qquad \alpha = \theta, \qquad \tau = t. \end{align}$$(2.2) Using the chain rule, the relationship between the partial derivatives in the original coordinate system and the new coordinate system is as follows. $$\begin{align*}& \frac{\partial}{\partial r} = \frac{2r}{R_2^2(0)}\frac{\partial}{\partial \zeta}, \qquad \frac{\partial}{\partial \theta} = \frac{\partial}{\partial \alpha}, \qquad \frac{\partial}{\partial t} = \frac{\partial}{\partial \tau} - \frac{Q}{\pi R_2^2(0)} \frac{\partial}{\partial \zeta}. \end{align*}$$ The basic solution after the coordinate change is |$(u_{\zeta },u_{\alpha }) = (0,0)$| with the interfaces stationary at |$\zeta = \zeta _1:= R_1^2(0)/R_2^2(0)$| and |$\zeta = 1$|⁠. |$\mu = \mu (\zeta )$| is now independent of time. We perturb this basic solution |$ (u_{\zeta } = 0, u_{\alpha } = 0, p_b, \mu)$| by |$ (\tilde{u}_{\zeta },\tilde{u}_{\alpha },\tilde{p}, \tilde{\mu } )$| where the disturbances are assumed to be small. The linearized equations which govern these disturbances are $$\begin{align} \begin{cases} & \frac{\partial \tilde{u}_{\zeta}}{\partial \zeta} + \frac{1}{\zeta} \frac{\partial \tilde{u}_{\alpha}}{\partial \alpha} = 0 \\ &\frac{\partial \tilde{p}}{\partial \zeta} = -\frac{R_2^4(0)}{4(\zeta R_2^2(0) + R_0^2(\tau))} \mu \tilde{u}_{\zeta} - \frac{QR_2^2(0)}{4 \pi (\zeta R_2^2(0) + R_0^2(\tau))} \tilde{\mu} \\ &\frac{\partial \tilde{p}}{\partial \alpha} = - \frac{\zeta R_2^2(0) + R_0^2(\tau)}{\zeta} \mu \tilde{u}_{\alpha} \\ \frac{\partial \tilde{\mu}}{\partial \tau} + \tilde{u}_{\zeta} \frac{\partial \mu}{\partial \zeta} = 0. \end{cases} \end{align}$$(2.3) We use separation of variables and decompose the disturbances into Fourier modes in the |$\alpha $| coordinate so the disturbances are of the form $$\begin{align} & \big(\tilde{u}_{\zeta},\tilde{u}_{\alpha},\tilde{p}, \tilde{\mu}\big) = \big(f(\zeta,\tau), \delta(\zeta,\tau), \psi(\zeta,\tau),\phi(\zeta,\tau)\big) e^{in\alpha}. \end{align}$$(2.4) Using this ansatz in equation (2.3) yields the following relations: $$\begin{align} &\frac{\partial \phi(\zeta,\tau)}{\partial \tau} = - \frac{d \mu}{d \zeta} f(\zeta,\tau), \end{align}$$(2.5) $$\begin{align} &\frac{\partial}{\partial \zeta} \left\{ \left( \zeta R_2^2(0) + R_0^2(\tau) \right) \mu \frac{\partial f(\zeta,\tau)}{\partial \zeta} \right\} -\frac{n^2R_2^4(0)}{4(\zeta R_2^2(0) + R_0^2(\tau))} \mu f(\zeta,\tau) \\ = &\frac{Qn^2R_2^2(0)}{4 \pi (\zeta R_2^2(0) + R_0^2(\tau))} \phi(\zeta,\tau).\notag \end{align}$$(2.6) In the innermost and outermost layers, the viscosity is constant and therefore |$\phi (\zeta ,\tau ) \equiv 0$|⁠. In these regions, the solution of (2.6) is of the form $$\begin{align} & f(\zeta,\tau) = \widetilde{C_1} \left(\zeta R_2^2(0) + R_0^2(\tau)\right)^{\frac{n}{2}} + \widetilde{C_2} \left(\zeta R_2^2(0) + R_0^2(\tau)\right)^{-\frac{n}{2}}. \end{align}$$(2.7) 2.1 Interface Conditions Recall that the inner interface in the |$\zeta $|-coordinate system is located at |$\zeta = \zeta _1:= R_1^2(0)/R_2^2(0)$|⁠. Let the disturbance of this interface be of the form |$\eta _i = C_n(\tau ) e^{in\alpha }$|⁠. The linearized kinematic condition at the inner interface is given by $$\begin{align} & C_n^{\prime}(\tau) = f(\zeta_1,\tau). \end{align}$$(2.8) The outer interface is located at |$\zeta = 1$|⁠. If the disturbance of the interface is of the form |$\eta _o = D_n(\tau ) e^{in\alpha }$|⁠, then the linearized kinematic condition is $$\begin{align} & D_n^{\prime}(\tau) = f(1,\tau). \end{align}$$(2.9) The linearized dynamic interface condition at the inner interface is given by $$\begin{align*}& \frac{2 R_1^2(\tau)}{R_2^2(0)} \big( \tilde{p}^+(\zeta_1) - \tilde{p}^-(\zeta_1) \big) - \eta_i \frac{Q}{2 \pi} \big(\mu(\zeta_1) - \mu_i \big) = T_1 \left( \frac{\eta_i + \frac{\partial^2 \eta_i}{\partial \alpha^2}}{R_1(\tau)}\right), \end{align*}$$ where |$T_1$| is the interfacial tension (see Gin (2015) for the derivation). Using the ansatz (2.4) and the system (2.3), $$\begin{align} & \frac{2 R_1^2(\tau)}{R_2^2(0)} \big(\mu_i (f^-)^{\prime}(\zeta_1,\tau) - \mu(\zeta_1) (f^+)^{\prime}(\zeta_1,\tau) \big) = \left\{ \frac{Qn^2}{2 \pi R_1^2(\tau)} \big(\mu(\zeta_1) - \mu_i \big) - T_1 \frac{n^4 - n^2}{R_1^3(\tau)}\right\} C_n(\tau). \end{align}$$(2.10) When |$\zeta < \zeta _1$|⁠, |$f$| is of the form given by (2.7). When |$\tau = 0$|⁠, in order to avoid a singularity when |$\zeta \to 0$|⁠, |$f$| must be of the form |$f(\zeta ,\tau ) = \widetilde{C_1} (\zeta R_2^2(0)+ R_0^2(\tau ))^{\frac{n}{2}}.$| We assume this also to be true for |$\tau> 0$|⁠. Then $$\begin{align} & (f^-)^{\prime}(\zeta_1,\tau) = \frac{n R_2^2(0)}{2 R_1^2(\tau)} f(\zeta_1,\tau). \end{align}$$(2.11) Using (2.11) in (2.10), $$\begin{align} & n \mu_i f(\zeta_1,\tau) - \frac{2 R_1^2(\tau)}{R_2^2(0)} \mu(\zeta_1) (f^+)^{\prime}(\zeta_1,\tau) = \left\{ \frac{Qn^2}{2 \pi R_1^2(\tau)} \big(\mu(\zeta_1) - \mu_i \big) - T_1 \frac{n^4 - n^2}{R_1^3(\tau)}\right\} C_n(\tau). \end{align}$$(2.12) Combining this with the kinematic interface condition (2.8), $$\begin{align} & C_n^{\prime}(\tau) = \frac{2 R_1^2(\tau)}{n R_2^2(0)} \frac{\mu(\zeta_1)}{\mu_i} (f^+)^{\prime}(\zeta_1,\tau) + \frac{F_1}{\mu_i} C_n(\tau), \end{align}$$(2.13) where |$F_1$| is given by $$\begin{align} & F_1 = \frac{Qn}{2 \pi R_1^2(\tau)} \big(\mu(\zeta_1) - \mu_i \big) - T_1 \frac{n^3 - n}{R_1^3(\tau)}. \end{align}$$(2.14) A similar procedure for the outer interface yields the interface condition $$\begin{align} & D_n^{\prime}(\tau) = - \frac{2 R_2^2(\tau)}{n R_2^2(0)} \frac{\mu(1)}{\mu_o} (f^-)^{\prime}(1,\tau) + \frac{F_2}{\mu_o} D_n(\tau), \end{align}$$(2.15) where |$F_2$| is given by $$\begin{align} & F_2 = \frac{Qn}{2 \pi R_2^2(\tau)} \big(\mu_o - \mu(1) \big) - T_2 \frac{n^3 - n}{R_2^3(\tau)}, \end{align}$$(2.16) and |$T_2$| is the interfacial tension at the outer interface. 2.2 Eigenvalue problem To this point we have the following system of equations where the field equations hold in the domain |$(\zeta _1, 1)$|⁠: $$\begin{align} \begin{cases} & \frac{\partial \phi(\zeta,\tau)}{\partial \tau} = - \frac{d \mu}{d \zeta} f(\zeta,\tau), \\ & \frac{\partial}{\partial \zeta} \left\{ \left( \zeta R_2^2(0) + R_0^2(\tau) \right) \mu \frac{\partial f(\zeta,\tau)}{\partial \zeta} \right\} -\frac{n^2R_2^4(0)}{4(\zeta R_2^2(0) + R_0^2(\tau))} \mu f(\zeta,\tau) = \frac{Qn^2R_2^2(0)\phi(\zeta,\tau)}{4 \pi (\zeta R_2^2(0) + R_0^2(\tau))} \\ & C_n^{\prime}(\tau) = \frac{2 R_1^2(\tau)}{n R_2^2(0)} \frac{\mu(\zeta_1)}{\mu_i} f^{\prime}(\zeta_1) + \frac{F_1}{\mu_i} C_n(\tau) \\ & D_n^{\prime}(\tau) = - \frac{2 R_2^2(\tau)}{n R_2^2(0)} \frac{\mu(1)}{\mu_o} f^{\prime}(1) + \frac{F_2}{\mu_o} D_n(\tau), \end{cases} \end{align}$$(2.17) where we have dropped the superscripts “+” and “-”. Using a quasi-steady-state approximation (QSSA) in which |$\tau $| (and hence |$R(\tau )$|⁠) are frozen, the functions |$\phi (x,\tau )$|⁠, |$C_n(\tau )$| and |$D_n(\tau )$| experience short-time exponential growth satisfying $$\begin{align} \begin{cases} &\frac{\partial \phi(\zeta,\tau)}{\partial \tau} = \sigma(\tau) \phi(\zeta,\tau), \\ &C_n^{\prime}(\tau) = \sigma(\tau) C_n(\tau), \\ &D_n^{\prime}(\tau) = \sigma(\tau) D_n(\tau), \end{cases} \end{align}$$(2.18) for some growth rate |$\sigma (\tau )$|⁠. Plugging (2.18) into (2.17) and using (2.8) and (2.9), |$(f,\sigma )$| is a solution to the following eigenvalue problem in the domain |$(\zeta _1, 1)$|⁠: $$\begin{align} \begin{cases} & \Big(\left(\zeta R_2^2(0) + R_0^2(\tau) \right) \mu f^{\prime}(\zeta) \Big)^{\prime} - \frac{n^2 R_2^4(0)}{4(\zeta R_2^2(0) + R_0^2(\tau))} \mu f(\zeta) = - \frac{Q n^2 R_2^2(0)}{4 \pi (\zeta R_2^2(0) + R_0^2(\tau))} \frac{1}{\sigma(\tau)} \frac{d \mu}{d \zeta} f(\zeta), \\ & \frac{2 R_1^2(\tau)}{n R_2^2(0)} \mu(\zeta_1) f^{\prime}(\zeta_1) = \left(\mu_i - \frac{F_1}{\sigma(\tau)} \right) f(\zeta_1), \\ & -\frac{2 R_2^2(\tau)}{n R_2^2(0)} \mu(1) f^{\prime}(1) = \left(\mu_o - \frac{F_2}{\sigma(\tau)} \right) f(1). \end{cases} \end{align}$$(2.19) The eigenvalues of system (2.19) are the time-dependent growth rates of the disturbances of the system. The QSSA allows for the considerable analysis and computation of growth rates that follows. We provide an in-depth investigation into the validity of the QSSA in Section 7.4. 3. Relating the growth of interfacial disturbances in the |$\zeta $|-coordinates with the physical coordinates We now relate the growth of the interfacial disturbances in the |$\zeta $|-coordinates to the same in the physical coordinate system. We start with the inner interface. Recall that in the transformed coordinates, the inner interface was disturbed by |$C_n(\tau ) e^{in\alpha }$|⁠. Therefore, it is located at |$\zeta = \zeta _1 + C_n(\tau ) e^{in\alpha }$|⁠. Thus, the position of the interface in the physical coordinates is $$\begin{align*}& r = \sqrt{\zeta R_2^2(0) + R_0^2(t)} = \sqrt{(\zeta_1 + C_n(\tau) e^{in\alpha}) R_2^2(0) + R_0^2(t)}. \end{align*}$$ Expanding about |$\zeta = \zeta _1$|⁠, $$\begin{align*}& r = R_1(\tau) + \frac{R_2^2(0)}{2 R_1(\tau)} C_n(\tau) e^{in\alpha} + \mathcal{O}(C_n^2(\tau)). \end{align*}$$ If we write the disturbance in the physical coordinates as |$A_n(t) e^{in\theta }$| (that is, the interface is located at |$r = R_1(t) + A_n(t) e^{in\theta }$|⁠), then, within linear approximation, $$\begin{align} & A_n(t) = \frac{R_2^2(0)}{2 R_1(\tau)} C_n(\tau). \end{align}$$(3.1) This implies that $$\begin{align} & \frac{A_n^{\prime}(t)}{A_n(t)} = \frac{C_n^{\prime}(\tau)}{C_n(\tau)} -\frac{Q}{2\pi R_1^2(t)}. \end{align}$$(3.2) Following the same process, the growth rate of the disturbance of the outer interface is $$\begin{align} & \frac{B_n^{\prime}(t)}{B_n(t)} = \frac{D_n^{\prime}(\tau)}{D_n(\tau)} -\frac{Q}{2\pi R_2^2(t)}, \end{align}$$(3.3) where the outer interface is located at |$r = R_2(t) + B_n(t) e^{in\theta }$|⁠. 4. Constant Viscosity Fluids We now consider the case in which all of the fluids have constant viscosity, first for two-layer flows and then for three-layer. We do this to demonstrate that the variable viscosity formulation can recover previous results in the constant viscosity limit. Through this process, we also present some new results on three-layer constant viscosity flows. When there are only two fluids (i.e. one interface located at |$r = R(t)$|⁠), the above analysis holds with the coordinate transformation $$\begin{align*}& \zeta = \frac{r^2 - R_0^2(t)}{R^2(0)}. \end{align*}$$ In the new coordinates, the basic solution has the interface fixed at |$\zeta = 1$|⁠. Let |$\mu _i$| denote the viscosity of the inner fluid and |$\mu _o$| denote the viscosity of the outer fluid. Analogous to equations (2.9) and (2.10), the kinematic interface condition is $$\begin{align} & D_n^{\prime}(\tau) = f(1,\tau), \end{align}$$(4.1) and the dynamic interface condition is $$\begin{align} & \frac{2 R^2(\tau)}{R^2(0)} \big(-\mu_o (f^+)^{\prime}(1,\tau) + \mu_i (f^-)^{\prime}(1,\tau) \big) = \left\{ \frac{Qn^2}{2 \pi R^2(\tau)} \big(\mu_o - \mu_i \big) - T \frac{n^4 - n^2}{R^3(\tau)}\right\} D_n(\tau), \end{align}$$(4.2) where |$T$| is the interfacial tension and |$D_n(\tau )$| is the amplitude of the disturbance of the interface with wave number |$n$|⁠. Also, as stated in the derivation of the interface conditions above, $$\begin{align*}& f(\zeta,\tau) = \widetilde{C_1} \Big(\zeta R^2(0) + R_0^2(\tau)\Big)^{\frac{n}{2}}, \qquad \zeta < 1, \end{align*}$$ and $$\begin{align*}& f(\zeta,\tau) = \widetilde{C_2} \Big(\zeta R^2(0) + R_0^2(\tau)\Big)^{-\frac{n}{2}}, \qquad \zeta> 1. \end{align*}$$ Using these in equations (4.1) and (4.2) gives the two-layer growth rate $$\begin{align} & \sigma:= \frac{D_n^{\prime}(\tau)}{D_n(\tau)} = \frac{Qn}{2 \pi R^2(\tau)} \frac{\mu_o - \mu_i}{\mu_o + \mu_i} - \frac{T}{\mu_i + \mu_o} \frac{n^3 - n}{R^3(\tau)}. \end{align}$$(4.3) This is an expression for the growth rate of the disturbance of the interface in the |$\zeta $|-coordinate system. This problem can be solved in the original |$r$|-coordinate system, and the result is a classic one (Paterson, 1981). Paterson’s result, using the notation in this paper, is that a disturbance with wave number |$n$| and amplitude |$B_n(t)$| has the growth rate $$\begin{align} & \frac{B_n^{\prime}(t)}{B_n(t)} = \frac{Qn}{2 \pi R^2(t)} \frac{\mu_o - \mu_i}{\mu_o + \mu_i} - \frac{T}{\mu_o + \mu_i} \frac{n^3 - n}{R^3(t)} - \frac{Q}{2 \pi R^2(t)}. \end{align}$$(4.4) The relationship between equations (4.3) and (4.4) is consistent with the comparison of the growth rates in the two different coordinate systems given by equation (3.3). We now turn to three-layer flows in which the fluid in the middle layer also has constant viscosity, |$\mu _1$|⁠. This situation has been investigated in Beeson-Jones & Woods (2015) and Gin & Daripa (2018), and it has been found that the magnitudes of interfacial disturbances |$A_n(t)$| and |$B_n(t)$| are governed by the following system of ODE’s $$\begin{align} & \frac{d}{dt} \left(\begin{array}{c} A_n(t) \\ B_n(t) \end{array}\right) = {\textbf{M}}_1^r(t) \left(\begin{array}{c} A_n(t) \\ B_n(t) \end{array}\right), \end{align}$$(4.5) where |${\textbf{M}}_1^r(t)$| is the |$2 \times 2$| matrix with entries given by $$\begin{align} \begin{cases} \Big({\textbf{M}}_1^r(t)\Big)_{11} &= \frac{\left\{(\mu_o + \mu_1) - (\mu_o - \mu_1) \left(\frac{R_1}{R_2}\right)^{2} \right\} F_1} {(\mu_1 - \mu_i) (\mu_o - \mu_1) \left(\frac{R_1}{R_2}\right)^{2} + (\mu_1 + \mu_i) (\mu_o + \mu_1)} - \frac{Q}{2 \pi R_1^2}, \\ \Big({\textbf{M}}_1^r(t)\Big)_{12} &= \frac{2 \mu_1 \left(\frac{R_1}{R_2}\right)^{n-1} F_2}{(\mu_1 - \mu_i) (\mu_o - \mu_1) \left(\frac{R_1}{R_2}\right)^{2} + (\mu_1 + \mu_i) (\mu_o + \mu_1)}, \\ \Big({\textbf{M}}_1^r(t)\Big)_{21} &= \frac{2 \mu_1 \left(\frac{R_1}{R_2}\right)^{n+1} F_1}{(\mu_1 - \mu_i) (\mu_o - \mu_1) \left(\frac{R_1}{R_2}\right)^{2} + (\mu_1 + \mu_i) (\mu_o + \mu_1)}, \\ \Big({\textbf{M}}_1^r(t)\Big)_{22} &= \frac{\left\{(\mu_1 + \mu_i) + (\mu_1 - \mu_i) \left(\frac{R_1}{R_2}\right)^{2} \right\} F_2} {(\mu_1 - \mu_i) (\mu_o - \mu_1) \left(\frac{R_1}{R_2}\right)^{2} + (\mu_1 + \mu_i) (\mu_o + \mu_1)} - \frac{Q}{2 \pi R_2^2}. \end{cases} \end{align}$$(4.6) We recall equation (3.1), which compares the interfacial disturbance of the inner interface in the |$r$|-coordinates and the |$\zeta $|-coordinates and also consider the corresponding equation for the outer interface: $$\begin{align} & A_n(t) = \frac{R_2^2(0)}{2 R_1(\tau)} C_n(\tau), \qquad B_n(t) = \frac{R_2^2(0)}{2 R_2(\tau)} D_n(\tau). \end{align}$$(4.7) Equations (4.5) and (4.7) give the matrix equation $$\begin{align} & \frac{d}{dt} \left(\begin{array}{c} A_n(t) \\ B_n(t) \end{array}\right) = {\textbf{M}}_1^r(t) \frac{R_2^2(0)}{2} {\textbf{R}}^{-1} \left(\begin{array}{c} C_n(\tau) \\ D_n(\tau) \end{array}\right), \end{align}$$(4.8) where $$\begin{align}& \left( {\textbf{R}} = \begin{array}{cc} R_1 & 0 \\ 0 & R_2 \end{array}\right). \end{align}$$(4.9) Taking derivatives of (4.7) and rewriting the resulting equations in matrix form, we obtain $$\begin{align} & \frac{d}{dt} \left(\begin{array}{c} A_n(t) \\ B_n(t) \end{array}\right) = \frac{R_2^2(0)}{2} {\textbf{R}}^{-1} \frac{d}{d \tau} \left(\begin{array}{c} C_n(\tau) \\ D_n(\tau) \end{array}\right) - \frac{R_2^2(0)}{2} {\textbf{R}}^{-1} {\textbf{Q}} \left(\begin{array}{c} C_n(\tau) \\ D_n(\tau) \end{array}\right), \end{align}$$(4.10) where $$\begin{align}& {\textbf{Q}} = \left(\begin{array}{cc} \frac{Q}{2 \pi R_1^2} & 0 \\ 0 & \frac{Q}{2 \pi R_2^2} \end{array}\right). \end{align}$$(4.11) Combining (4.8) and (4.10), $$\begin{align} & \frac{d}{d\tau} \left(\begin{array}{c} C_n(\tau) \\ D_n(\tau) \end{array}\right) = {\textbf{M}}_1^{\zeta}(\tau) \left(\begin{array}{c} C_n(\tau) \\ D_n(\tau) \end{array}\right), \end{align}$$(4.12) where |${\textbf{M}}_1^{\zeta } = {\textbf{R}} {\textbf{M}}_1^r {\textbf{R}}^{-1} + {\textbf{Q}}$|⁠. The entries of |${\textbf{M}}_1^{\zeta }(\tau )$| are $$\begin{align} \begin{cases} \Big({\textbf{M}}_1^{\zeta}(\tau)\Big)_{11} &= \frac{\left\{(\mu_o + \mu_1) - (\mu_o - \mu_1) \left(\frac{R_1}{R_2}\right)^{2} \right\} F_1} {(\mu_1 - \mu_i) (\mu_o - \mu_1) \left(\frac{R_1}{R_2}\right)^{2} + (\mu_1 + \mu_i) (\mu_o + \mu_1)}, \\ \Big({\textbf{M}}_1^{\zeta}(\tau)\Big)_{12} &= \frac{2 \mu_1 \left(\frac{R_1}{R_2}\right)^{n} F_2}{(\mu_1 - \mu_i) (\mu_o - \mu_1) \left(\frac{R_1}{R_2}\right)^{2} + (\mu_1 + \mu_i) (\mu_o + \mu_1)}, \\ \Big({\textbf{M}}_1^{\zeta}(\tau)\Big)_{21} &= \frac{2 \mu_1 \left(\frac{R_1}{R_2}\right)^{n} F_1}{(\mu_1 - \mu_i) (\mu_o - \mu_1) \left(\frac{R_1}{R_2}\right)^{2} + (\mu_1 + \mu_i) (\mu_o + \mu_1)}, \\ \Big({\textbf{M}}_1^{\zeta}(\tau)\Big)_{22} &= \frac{\left\{(\mu_1 + \mu_i) + (\mu_1 - \mu_i) \left(\frac{R_1}{R_2}\right)^{2} \right\} F_2} {(\mu_1 - \mu_i) (\mu_o - \mu_1) \left(\frac{R_1}{R_2}\right)^{2} + (\mu_1 + \mu_i) (\mu_o + \mu_1)}. \end{cases} \end{align}$$(4.13) There are several important facts about the relationship between the matrices |${\textbf{M}}_1^{\zeta }$| and |${\textbf{M}}_1^r$|⁠: It was demonstrated in Gin & Daripa (2015) that |${\textbf{M}}_1^r$| can have complex eigenvalues. However, it is shown in the next section (Section 5) that the problem in the |$\zeta $|-coordinates has real growth rates. This analysis holds even for constant viscosity. Therefore, |${\textbf{M}}_1^{\zeta }$| has real eigenvalues. |${\textbf{M}}_1^{\zeta }$| can also be expressed as |${\textbf{M}}_1^{\zeta } = {\textbf{R}} \left ({\textbf{M}}_1^r + {\textbf{Q}} \right ) {\textbf{R}}^{-1}$|⁠. Since a similarity transformation does not change eigenvalues, the eigenvalues of |${\textbf{M}}_1^{\zeta }$| are the eigenvalues of |${\textbf{M}}_1^r + {\textbf{Q}}$| where |${\textbf{Q}}$| is a diagonal matrix. Thus it is this diagonal matrix |${\textbf{Q}}$| which, when added to |${\textbf{M}}_1^r$|⁠, converts the complex eigenvalues to real and leaves real eigenvalues real. Both |${\textbf{M}}_1^r$| and |${\textbf{M}}_1^{\zeta }$| have real eigenvalues when |$F_1$| and |$F_2$| defined respectively in (2.15) and (2.17) have the same sign. Define the matrices: $$\begin{align} &{\textbf{E}} = \left(\begin{array}{cc} R_1 \sqrt{|F_1|} & 0 \\ 0 & R_2 \sqrt{|F_2|} \end{array}\right), \qquad {\textbf{F}} = \left(\begin{array}{cc} \sqrt{|F_1|} & 0 \\ 0 & \sqrt{|F_2|} \end{array}\right) \end{align}$$(4.14) then |${\textbf{E}} {\textbf{M}}_1^r {\textbf{E}}^{-1}$| and |${\textbf{F}} {\textbf{M}}_1^{\zeta } {\textbf{F}}^{-1}$| are real symmetric matrices. Therefore, |${\textbf{M}}_1^r$| and |${\textbf{M}}_1^{\zeta }$| are similar to real symmetric (i.e. self-adjoint) matrices and have real eigenvalues. However, the argument breaks down when |$F_1$| and |$F_2$| have opposite signs because |${\textbf{E}} {\textbf{M}}_1^r {\textbf{E}}^{-1}$| and |${\textbf{F}} {\textbf{M}}_1^{\zeta } {\textbf{F}}^{-1}$| are not symmetric. This shows that the complex eigenvalues of |${\textbf{M}}_1^r$| can only occur when |$F_1 F_2 \leq 0$|⁠. 5. Upper Bounds We follow a process similar to Daripa & Pasa (2006) in order to derive an upper bound on the growth rate. First we take an inner product of (2.19)1 with |$f$|⁠. Using integration by parts along with the boundary conditions (2.19)2 and (2.19)3 and solving for |$\sigma $| yields $$\begin{align} & \sigma = \frac{2 \pi n R_2^2(0) F_1 |f(\zeta_1)|^2 + 2 \pi n R_2^2(0) F_2 |f(1)|^2 + Qn^2R_2^2(0) I_1}{2 \pi n R_2^2(0)\mu_i |f(\zeta_1)|^2 + 2 \pi n R_2^2(0) \mu_o|f(1)|^2 + 4 \pi I_2+ \pi n^2 R_2^4(0)I_3}, \end{align}$$(5.1) where $$\begin{align} I_1 &= \int_{\zeta_1}^1 \frac{\mu^{\prime}(\zeta)}{\zeta R_2^2(0) + R_0^2(\tau)} |f(\zeta)|^2 d\zeta, \end{align}$$(5.2) $$\begin{align} I_2 &= \int_{\zeta_1}^{1} \left(\zeta R_2^2(0) + R_0^2(\tau) \right) \mu(\zeta) |f^{\prime}(\zeta)|^2 d\zeta, \end{align}$$(5.3) $$\begin{align} I_3 &= \int_{\zeta_1}^1 \frac{\mu(\zeta)}{\zeta R_2^2(0) + R_0^2(\tau)} |f(\zeta)|^2 d\zeta. \end{align}$$(5.4) Note that all terms in (5.1) are real. Therefore, ‘|$\sigma $| is real for all wave numbers.’ This is a product of the change of variables from the |$r$|-coordinates to the |$\zeta $|-coordinates. It is shown in Gin & Daripa (2015) that the growth rate can be complex for constant viscosity flows in the |$r$|-coordinates. When |$\sigma> 0$|⁠, we may ignore the positive term containing |$I_2$| in the denominator and get $$\begin{align*}& \sigma < \frac{2 \pi n R_2^2(0) F_1 |f(\zeta_1)|^2 + 2 \pi n R_2^2(0) F_2 |f(1)|^2 + Qn^2R_2^2(0) I_1}{2 \pi n R_2^2(0)\mu_i |f(\zeta_1)|^2 + 2 \pi n R_2^2(0) \mu_o |f(1)|^2 + \pi n^2 R_2^4(0)I_3}. \end{align*}$$ We use the following inequality $$\begin{align} & \frac{\sum\limits_{i=1}^{N} A_ix_i}{\sum\limits_{i=1}^{N} B_ix_i} \leq \max_i \left\{\frac{A_i}{B_i}\right\}, \end{align}$$(5.5) which holds for any |$N$| if |$A_i> 0$|⁠, |$B_i> 0$| and |$x_i> 0$| for all |$i = 1,...,N$|⁠. By using this inequality with |$N = 3$|⁠, $$\begin{align*}& \sigma < \max \left\{\frac{F_1}{\mu_i}, \frac{F_2}{\mu_o}, \frac{Q}{\pi R_2^2(0)} \frac{I_1}{I_3} \right\}. \end{align*}$$ But $$\begin{align*} \frac{I_1}{I_3} < \frac{{\sup_{\zeta \in (\zeta_1,1)}} \mu^{\prime}(\zeta)} {{\displaystyle\inf_{\zeta \in (\zeta_1,1)}} \mu(\zeta)} < \frac{{\displaystyle\sup_{\zeta \in (\zeta_1,1)}} \mu^{\prime}(\zeta)}{\mu_i}. \end{align*}$$ Therefore, $$\begin{align} & \sigma < \max \left\{\frac{F_1}{\mu_i}, \frac{F_2}{\mu_o}, \frac{Q}{\pi R_2^2(0)} \frac{1}{\mu_i} \sup_{\zeta \in (\zeta_1,1)} \mu^{\prime}(\zeta) \right\}. \end{align}$$(5.6) Using the definitions of |$F_1$| and |$F_2$| given by (2.14) and (2.16), $$\begin{align} \sigma < \max & \left\{ \frac{Qn}{2 \pi R_1^2(\tau)}\left(\frac{\mu(\zeta_1) - \mu_i}{\mu_i}\right) - \frac{T_1}{\mu_i} \frac{n^3-n}{R_1^3(\tau)}, \right.\nonumber \\ & \left. \frac{Qn}{2 \pi R_2^2(\tau)}\left(\frac{\mu_o - \mu(1)}{\mu_o}\right) - \frac{T_2}{\mu_o} \frac{n^3-n}{R_2^3(\tau)}, \frac{Q}{\pi R_2^2(0)} \frac{1}{\mu_i} \sup_{\zeta \in (\zeta_1,1)} \mu^{\prime}(\zeta) \right\}, \end{align}$$(5.7) which is the modal upper bound for a wave with wave number |$n$|⁠. We can find an absolute upper bound for all wave numbers by taking the maximum of the first two terms over all values of |$n$|⁠. The absolute upper bound is $$\begin{align} \sigma < \max & \left\{\frac{2 T_1}{\mu_i R_1^3(\tau)} \left(\frac{Q R_1(\tau)}{6 \pi T_1} (\mu(\zeta_1) - \mu_i) + \frac{1}{3}\right)^{\frac{3}{2}}, \right. \nonumber\\ & \left. \frac{2 T_2}{\mu_o R_2^3(\tau)} \left(\frac{Q R_2(\tau)}{6 \pi T_2} (\mu_o - \mu(1)) + \frac{1}{3}\right)^{\frac{3}{2}}, \frac{Q}{\pi R_2^2(0)} \frac{1}{\mu_i} \sup_{\zeta \in (\zeta_1,1)} \mu^{\prime}(\zeta) \right\}. \end{align}$$(5.8) 6. Characterization of the Eigenvalues and Eigenfunctions Using |$\lambda = 1/\sigma $|⁠, the eigenvalue problem (2.19) can be written as $$\begin{align} \begin{cases} & \Big(\left(\zeta R_2^2(0) + R_0^2(\tau) \right) \mu f^{\prime}(\zeta) \Big)^{\prime} - \left( \frac{n^2 R_2^4(0)}{4(\zeta R_2^2(0) + R_0^2(\tau))} \mu -\frac{Q n^2 R_2^2(0)}{4 \pi (\zeta R_2^2(0) + R_0^2(\tau))} \mu^{\prime} \lambda \right) f(\zeta) = 0, \\ & \left(\mu_i - \lambda F_1 \right) f(\zeta_1) - \frac{2 R_1^2(\tau)}{Rn_2^2(0)} \mu(\zeta_1) f^{\prime}(\zeta_1) = 0, \\ & \left(\mu_o - \lambda F_2 \right) f(1) + \frac{2 R_2^2(\tau)}{Rn_2^2(0)} \mu(1) f^{\prime}(1) = 0. \end{cases} \end{align}$$(6.1) This is a Sturm–Liouville eigenvalue problem with eigenvalue-dependent boundary conditions. Because the boundary conditions depend on the eigenvalues and are not the natural boundary conditions, the classic Sturm–Liouville theory cannot immediately be applied. In what follows, we examine conditions for which the properties of regular Sturm–Liouville problems hold and determine a function space for which the operator is self-adjoint. Note that |$F_1$| and |$F_2$| are positive for small values of |$n$| and negative for large values of |$n$| (see equations (2.14) and (2.16)). In particular, both |$F_1$| and |$F_2$| are positive when $$\begin{align}& n < \min \left\{\sqrt{\frac{Q R_1(\tau)}{2 \pi T_1} (\mu(\zeta_1) - \mu_i) + 1}, \sqrt{\frac{Q R_2(\tau)}{2 \pi T_2} (\mu_o - \mu(1)) + 1} \right\}. \end{align}$$(6.2) From the upper bound (5.6), we can see that as long as the viscosity gradient |$\mu ^{\prime}(\zeta )$| is not too large, the maximum value of |$\sigma $| will occur when |$F_1$| and |$F_2$| are positive. For this range of wave numbers and for monotonically increasing viscosity profiles, we have the following characterization of the eigenvalues and eigenfunctions. Theorem 6.1 Let |$F_1$|⁠, |$F_2$|⁠, |$Q$|⁠, |$n$|⁠, |$\mu _i$|⁠, |$\mu _o> 0$|⁠. Let |$\mu (\zeta )$| be a positive, strictly increasing function in |$C^1([\zeta _1,1])$|⁠. Then the eigenvalue problem (6.1) has a countably infinite number of real eigenvalues that can be ordered $$\begin{align*}& 0 < \lambda_0 < \lambda_1 < \lambda_2 <... \end{align*}$$ with the property that for the corresponding eigenfunctions, |$\left \{f_i\right \}_{i=0}^{\infty }$|⁠, |$f_i$| has exactly |$i$| zeros in the interval |$(\zeta _1,1)$|⁠. Additionally, the eigenfunctions are continuous with a continuous derivative. Proof. The fact that there are a countably infinite number of real eigenvalues that can be ordered and corresponding eigenfunctions with the prescribed number of zeros is proven in (Ince, 1956, p. 232–233) in Theorem I and Theorem II using $$\begin{align*}& \begin{cases} &a = \zeta_1, \qquad b = 1, \qquad K(x,\lambda) = \left(x R_2^2(0) + R_0^2(\tau) \right)\mu(x), \\ &G(x,\lambda) = \frac{n^2 R_2^4(0)}{4(x R_2^2(0) + R_0^2(\tau))} \mu(x) -\frac{Q n^2 R_2^2(0)}{4 \pi (x R_2^2(0) + R_0^2(\tau))} \mu^{\prime}(x) \lambda, \\ &\alpha = \frac{2 R_1^2(\tau)}{Rn_2^2(0)} \mu(\zeta_1), \qquad \alpha^{\prime} = \mu_i - \lambda F_1, \\ &\beta = \frac{2 R_2^2(\tau)}{Rn_2^2(0)} \mu(1), \qquad \beta^{\prime} = \mu_o - \lambda F_2. \end{cases} \end{align*}$$ The regularity of the eigenfunctions comes from the existence theorem of (Ince, 1956, p. 73). We saw from equation (5.1) that |$\sigma $| is real for all |$n$|⁠, and a closer look at each term in (5.1) shows that if |$F_1,F_2> 0$| and |$\mu (\zeta ), \mu ^{\prime}(\zeta )>0$|⁠, then all terms are positive and |$\sigma> 0$|⁠. 6.1 Self-Adjointness and Expansion Theorem We rewrite equation (6.1) as $$\begin{align} \begin{cases} & -\Big(\left(\zeta R_2^2(0) + R_0^2(\tau) \right) \mu f^{\prime}(\zeta) \Big)^{\prime} + \left( \frac{n^2 R_2^4(0)}{4(\zeta R_2^2(0) + R_0^2(\tau))} \right) \mu f(\zeta) = \frac{Q n^2 R_2^2(0)}{4 \pi (\zeta R_2^2(0) + R_0^2(\tau))} \mu^{\prime} \lambda f(\zeta), \\ & - \left(-\frac{\mu_i}{F_1} f(\zeta_1) + \frac{2 R_1^2(\tau)}{Rn_2^2(0) F_1} \mu(\zeta_1) f^{\prime}(\zeta_1)\right) = \lambda f(\zeta_1), \\ & -\left(-\frac{\mu_o}{F_2} f(1) - \frac{2 R_2^2(\tau)}{Rn_2^2(0)F_2} \mu(1) f^{\prime}(1)\right) = \lambda f(1). \end{cases} \end{align}$$(6.3) This is of the form $$\begin{align} \begin{cases} & Tf:= \frac{1}{r} \left\{-\left(p f^{\prime}\right)^{\prime} + qf \right\} = \lambda f, \qquad\qquad\qquad\qquad \qquad \zeta_1 < \zeta < 1, \\ & - \left(\beta_{11} f(\zeta_1) - \beta_{12} f^{\prime}(\zeta_1) \right) = \lambda \left( \alpha_{11} f(\zeta_1) - \alpha_{12} f^{\prime}(\zeta_1) \right), \\ &- \left(\beta_{21} f(1) - \beta_{22} f^{\prime}(1) \right) = \lambda \left( \alpha_{21} f(1) - \alpha_{22} f^{\prime}(1) \right), \end{cases} \end{align}$$(6.4) where $$\begin{align} \begin{cases} p(\zeta) = \left(\zeta R_2^2(0) + R_0^2(\tau) \right) \mu(\zeta), & \\ q(\zeta) = \frac{n^2 R_2^4(0) \mu(\zeta)}{4(\zeta R_2^2(0) + R_0^2(\tau))}, & \\ r(\zeta) = \frac{Q n^2 R_2^2(0) \mu^{\prime}(\zeta)}{4 \pi (\zeta R_2^2(0) + R_0^2(\tau))}, & \\ \beta_{11} = -\frac{\mu_i}{F_1}, & \beta_{12} = -\frac{2 R_1^2(\tau)}{Rn_2^2(0) F_1} \mu(\zeta_1), \\ \alpha_{11} = 1, & \alpha_{12} = 0, \\ \beta_{21} = -\frac{\mu_o}{F_2}, & \beta_{22} = \frac{2 R_2^2(\tau)}{Rn_2^2(0) F_2} \mu(1), \\ \alpha_{21} = 1, & \alpha_{22} = 0. \end{cases} \end{align}$$(6.5) Given the same assumptions as in Theorem 6.1, we have the following theorem from Walter (1973). Theorem 6.2 Let |$F_1$|⁠, |$F_2$|⁠, |$Q$|⁠, |$n$|⁠, |$\mu _i$|⁠, |$\mu _o> 0$|⁠. Let |$\mu (\zeta )$| be a positive, strictly increasing function in |$C^1([\zeta _1,1])$|⁠. Let |$p(\zeta )$|⁠, |$q(\zeta )$| and |$r(\zeta )$| be defined by (6.5). Let $$\begin{align*}& L^2_r(\zeta_1,1) = \left\{f(\zeta) \Big\vert \int_{\zeta_1}^{1} |f(\zeta)|^2 r(\zeta) d \zeta < \infty \right\}, \end{align*}$$ and define the operator |$T$| on |$L^2_r(\zeta _1,1)$| by $$\begin{align*}& Tf:= \frac{1}{r} \left\{-\left(p f^{\prime}\right)^{\prime} + q f\right\}. \end{align*}$$ Define the measure: $$\begin{align} \nu(M):= \begin{cases} \frac{nR_2^2(0)F_1}{2}, \qquad & \textrm{for}\ M = \{\zeta_1\} \\ \int_{M} r(\zeta) d\zeta, \qquad & \textrm{for}\ M \subset (\zeta_1,1) \\ \frac{nR_2^2(0)F_2}{2}, \qquad & \textrm{for}\ M = \{1\}. \end{cases} \end{align}$$(6.6) We consider the Hilbert space |$H:= L^2([\zeta _1,1];\nu )$|⁠. Consider the operator |$A$| with domain $$\begin{align}& D(A) = \{f \in H \vert f,f^{\prime} \textrm{ absolutely continuous in} (\zeta_1,1), T \in L^2_r(\zeta_1,1) \}, \end{align}$$(6.7) and defined by $$\begin{align} (Af)(\zeta) = \begin{cases} \displaystyle{\lim_{\zeta \to \zeta_1}} \left(\frac{\mu_i}{F_1} f(\zeta) - \frac{2 R_1^2(\tau)}{Rn_2^2(0) F_1} \mu(\zeta_1) f^{\prime}(\zeta)\right), \qquad & \textrm{if} \zeta = \{\zeta_1\} \\ (Tf)(\zeta), \qquad & \textrm{if} \zeta \in (\zeta_1,1) \\ \displaystyle{\lim_{\zeta \to 1}} \left(\frac{\mu_o}{F_2} f(\zeta) + \frac{2 R_2^2(\tau)}{Rn_2^2(0) F_2} \mu(1) f^{\prime}(\zeta)\right), \qquad & \textrm{if} \zeta = \{1\}. \end{cases} \end{align}$$(6.8) Then |$(f,\lambda )$| satisfies (6.3) if and only if |$Af = \lambda f$|⁠. |$A$| is a self-adjoint operator on |$H$| and for any |$u \in H$|⁠, $$\begin{align} & u = \sum_{k=0}^{\infty} f_k \int_{\zeta_1}^{1} u(\zeta) f_k(\zeta) d \nu, \end{align}$$(6.9) where the |$f_k$| are the eigenfunctions of |$A$|⁠. 7. Numerical Results We now investigate the growth rate of disturbances by numerically computing the eigenvalues of the eigenvalue problem (6.1). This eigenvalue problem has time-dependent coefficients and boundary conditions which depend on the eigenvalues. Thus the dispersion relation for this problem depends on time. The eigenvalues are computed using a pseudo-spectral Chebyshev method. The eigenvalues |$\lambda $| are then inverted to find the growth rates |$\sigma $|⁠. Recall that for a given wave number |$n$|⁠, there are infinitely many eigenvalues. In the results that follow, |$\sigma $| refers to the maximum over all eigenvalues. Because we are considering a circular domain, which is |$2 \pi $|-periodic, the wave number |$n$| only takes on integer values. However, in order to show smooth curves in the plots and make clear how |$\sigma $| varies with |$n$|⁠, we plot the dispersion relation for all values of |$n$|⁠. |$\sigma _{max}$| refers to the maximum over all eigenvalues and over all wave numbers. For consistency, we often use the same parameter values throughout our results. Unless otherwise stated, |$\mu _i = 2$|⁠, |$\mu _o = 10$|⁠, |$T_1 = T_2 = 1$| and |$Q = 10$|⁠. Therefore, the inner and outer layer fluids have constant viscosity 2 and 10 respectively for all our studies here. The viscosity profile of the middle layer fluid, however, is a free variable which can be taken as constant or variable in our studies below. In the rest of the paper, we will characterize the flow by the viscosity of the middle layer. 7.1 Constant versus Variable Viscosity We begin by comparing the growth rate of disturbances for a constant viscosity profile with that for a variable viscosity profile. In Fig. 2a, the dispersion relations at a fixed time (⁠|$\tau = 0$|⁠) are plotted for four different viscosity profiles shown in Fig. 2b. The constant viscosity case is given by the solid (black) line and the viscosity of the middle layer fluid is |$\mu = 6$|⁠. The stability of three-layer constant viscosity flows has been studied extensively in Gin & Daripa (2015, 2018). Note that there is a maximum growth rate and that short waves are stable due to interfacial tension. For comparison, three linear viscosity profiles are considered. The dotted (red) line corresponds to a linear viscosity profile with |$\mu (R_1) = 5.9$| and |$\mu (R_2) = 6.1$|⁠, the dashed (blue) line corresponds to |$\mu (R_1) = 5$| and |$\mu (R_2) = 7$|⁠, and the dash-dot (green) line corresponds to |$\mu (R_1) = 4$| and |$\mu (R_2) = 8$|⁠. There are several important features to notice. First, the dispersion relation for each of the variable viscosity profiles has a local maximum for a wave number that is similar to the maximum for the case of a constant viscosity profile. For profiles with smaller viscosity jumps at the interfaces, the local maximum is smaller. Therefore, this local maximum can be attributed to the instability of the interfaces due to the positive viscosity jump. The second thing to notice is that short waves are unstable for variable viscosity profiles, even when the viscosity profile is nearly constant (see the dotted line). As the gradient of the viscosity profiles increase, the growth rate of short waves also increases. Therefore, the short wave behavior is dominated by the instability of the middle layer fluid itself due to the viscosity gradient. The final observation which can be drawn from the dispersion relations is that the maximum growth rate can be smaller for a variable viscosity profile than it is for a constant viscosity profile with constant viscosity equal to the average of the values of viscosity at the two interfaces in the middle layer of the variable viscosity profile. Also shown in Fig. 2a are the modal upper bounds given by equation (5.7) associated with each of the four viscosity profiles. The upper bounds for each viscosity profile are the same color as the corresponding dispersion relation. The qualitative shapes of the modal upper bounds match the dispersion relations, but the upper bounds are not sharp. This is due to ignoring the |$I_2$| term in equation (5.1) and using the inequality (5.5) in the derivation of the upper bounds. There is a possibility that these upper bounds can be improved. Fig. 2. Open in new tabDownload slide A comparison of dispersion relations for four linear viscosity profiles. (a)The dispersion relations (⁠|$\sigma $| versus |$n$|⁠) at time |$\tau = 0$| as curves and the associated upper bounds given by equation (5.7) as markers of the same color. (b) The associated viscosity profiles. The parameter values are |$Q = 10$|⁠, |$\mu _i = 2$|⁠, |$\mu _o = 10$|⁠, |$T_1 = T_2 = 1$|⁠, |$R_1(0) = 20$| and |$R_2(0) = 30$|⁠. Fig. 2. Open in new tabDownload slide A comparison of dispersion relations for four linear viscosity profiles. (a)The dispersion relations (⁠|$\sigma $| versus |$n$|⁠) at time |$\tau = 0$| as curves and the associated upper bounds given by equation (5.7) as markers of the same color. (b) The associated viscosity profiles. The parameter values are |$Q = 10$|⁠, |$\mu _i = 2$|⁠, |$\mu _o = 10$|⁠, |$T_1 = T_2 = 1$|⁠, |$R_1(0) = 20$| and |$R_2(0) = 30$|⁠. 7.2 Optimal Profile In Section 7.1, it is shown that some particular variable viscosity profiles are less unstable (i.e. have a smaller maximum growth rate) than a particular constant viscosity profile. This leads to some more general questions: Are there variable viscosity profiles that are less unstable than all constant viscosity profiles? What is the optimal viscosity profile? The question of the optimal viscosity profile is a difficult one so we start by using some simplifying assumptions. First, recall that the viscosity profile and the growth rate are both time-dependent. In this section we only consider the growth rate at time |$\tau = 0$|⁠. This is reasonable because it is advantageous to control the instability at early times. Therefore, for the present discussion the term “optimal” refers to the viscosity profile that minimizes the maximum growth rate |$\sigma _{max}$| at time |$\tau = 0$|⁠. The second simplification is that we first consider only viscosity profiles in the middle layer that are linear at time |$\tau = 0$|⁠. Note that linear profiles are uniquely determined by the values |$\mu (R_1)$| and |$\mu (R_2)$|⁠. Other types of viscosity profiles will be considered later in this section. Fig. 2a shows the value of |$\sigma _{max}$| for each linear viscosity profile such that the viscosity of the middle layer is between |$\mu _i$| and |$\mu _o$|⁠. The optimal viscosity profile in this case is |$\mu (R_1) = 3.41$| and |$\mu (R_2) = 5.09$|⁠. Note that all possible constant viscosity profiles have been considered as a subset of the set of linear viscosity profiles. Therefore, the fact that the optimal profile is not constant shows that variable viscosity profiles can be used to reduce the instability of a flow. The dispersion curve for the optimal viscosity profile shown in Fig. 2b approaches the value indicated by the dotted horizontal line as |$n \to \infty $|⁠. The two local maxima in this plot have the same value as this limit. This is because the optimal viscosity profile is the one which balances the instabilities of the interfaces with the instability of the middle layer. Next we investigate the optimal viscosity profile under several different values of interfacial tension. Plots of |$\sigma _{max}$| versus the different linear profiles is given in Fig. 4. Figure 4a has the smallest value of interfacial tension with |$T_1 = T_2 = 0.25$|⁠. The optimal viscosity profile has endpoint viscosities of |$\mu (R_1) = 3.20$| and |$\mu (R_2) = 5.65$|⁠. Figure 4b uses |$T_1 = T_2 = 1$| and is a repeat of Fig. 3a. As noted above, the optimal profile is |$\mu (R_1) = 3.41$| and |$\mu (R_2) = 5.09$|⁠. Figure 4c has the largest values of interfacial tension with |$T_1 = T_2 = 4$|⁠. The optimal linear viscosity profile is |$\mu (R_1) = 3.47$| and |$\mu (R_2) = 4.53$|⁠. The trend is that larger values of interfacial tension correspond to optimal viscosity profiles with a smaller viscosity gradient. This is because, as mentioned previously, the optimal viscosity profile is the one which balances the instabilities of the interfaces with the instability of the middle layer. A larger value of interfacial tension decreases the instability of the interfaces so the gradient of the middle layer must also decrease in order to match the interfacial instability. Fig. 3. Open in new tabDownload slide (a) The value of |$\sigma _{max}$| for each linear viscosity profile which is defined by the values |$\mu (R_1)$| and |$\mu (R_2)$|⁠. The optimal profile is |$\mu (R_1) = 3.41$| and |$\mu (R_2) = 5.09$| and its dispersion relation is given in (b). The parameter values are |$Q = 10$|⁠, |$\mu _i = 2$|⁠, |$\mu _o = 10$|⁠, |$T_1 = T_2 = 1$|⁠, |$R_1(0) = 20$| and |$R_2(0) = 30$|⁠. Fig. 3. Open in new tabDownload slide (a) The value of |$\sigma _{max}$| for each linear viscosity profile which is defined by the values |$\mu (R_1)$| and |$\mu (R_2)$|⁠. The optimal profile is |$\mu (R_1) = 3.41$| and |$\mu (R_2) = 5.09$| and its dispersion relation is given in (b). The parameter values are |$Q = 10$|⁠, |$\mu _i = 2$|⁠, |$\mu _o = 10$|⁠, |$T_1 = T_2 = 1$|⁠, |$R_1(0) = 20$| and |$R_2(0) = 30$|⁠. Fig. 4. Open in new tabDownload slide Plots of |$\sigma _{max}$| for all linear viscosity profiles for three different values of interfacial tension: (a) |$T_1 = T_2 = 0.25$|⁠, (b) |$T_1 = T_2 = 1$| and (c) |$T_1 = T_2 = 4$|⁠. Other parameter values are |$Q = 10$|⁠, |$\mu _i = 2$|⁠, |$\mu _o = 10$|⁠, |$R_1(0) = 20$| and |$R_2(0) = 30$|⁠. Fig. 4. Open in new tabDownload slide Plots of |$\sigma _{max}$| for all linear viscosity profiles for three different values of interfacial tension: (a) |$T_1 = T_2 = 0.25$|⁠, (b) |$T_1 = T_2 = 1$| and (c) |$T_1 = T_2 = 4$|⁠. Other parameter values are |$Q = 10$|⁠, |$\mu _i = 2$|⁠, |$\mu _o = 10$|⁠, |$R_1(0) = 20$| and |$R_2(0) = 30$|⁠. In the case of chemical EOR by polymer flooding, the viscosity profile of the middle layer fluid is determined by the concentration of polymer. The use of large quantities of polymer can be expensive so it is useful to explore which viscosity profile minimizes the instability (i.e. |$\sigma _{\textrm{max}}$|⁠) given a fixed amount of total polymer. Assuming that there is a linear relationship between the concentration of polymer and the viscosity, this can be viewed as minimizing the instability for a fixed value of average viscosity. The results of this type of optimization are given in Fig. 5. The value of |$\sigma _{max}$| for all linear viscosity profiles is plotted in Fig. 5a using the same parameter values as Fig. 3a. For each possible average viscosity between |$\mu _i$| and |$\mu _o$|⁠, the profile which minimizes |$\sigma _{max}$| was found and is marked by an ‘|$x$|’ in Fig. 5a. The viscosity profile can be identified by its slope |$a = (\mu (R_2) - \mu (R_1))/(R_2-R_1)$|⁠. The slopes of the optimal profiles are plotted versus the average viscosity of the middle layer in Fig. 5b. Fig. 5. Open in new tabDownload slide Plots of the optimal linear viscosity profile for a fixed value of average viscosity of the middle layer fluid. (a) The value of |$\sigma _{max}$| versus the different linear viscosity profiles with |$x$|’s to denote the optimal profiles. (b) A plot of the slope of the optimal profile versus the average viscosity. The parameter values are |$Q = 10$|⁠, |$\mu _i = 2$|⁠, |$\mu _o = 10$|⁠, |$T_1 = T_2 = 1$|⁠, |$R_1(0) = 20$| and |$R_2(0) = 30$|⁠. Fig. 5. Open in new tabDownload slide Plots of the optimal linear viscosity profile for a fixed value of average viscosity of the middle layer fluid. (a) The value of |$\sigma _{max}$| versus the different linear viscosity profiles with |$x$|’s to denote the optimal profiles. (b) A plot of the slope of the optimal profile versus the average viscosity. The parameter values are |$Q = 10$|⁠, |$\mu _i = 2$|⁠, |$\mu _o = 10$|⁠, |$T_1 = T_2 = 1$|⁠, |$R_1(0) = 20$| and |$R_2(0) = 30$|⁠. When the average viscosity of the middle layer fluid is |$\mu _i = 2$|⁠, the optimal viscosity profile is constant at |$\mu \equiv 2$| (note that this is the only profile considered since we are taking |$\mu _i \leq \mu (r) \leq \mu _o$| for |$R_1 \leq r \leq R_2$|⁠). Therefore there is no viscosity jump at the inner interface and no instability in the layer itself. All of the instability occurs due to the jump in viscosity at the outer interface. As the average viscosity increases from there, the jump at the inner interface of the optimal viscosity profile remains zero while the viscosity gradient increases in order to decrease the viscosity jump at the outer interface. Eventually, as the average viscosity nears |$\mu = 3$|⁠, a point is reached in which the gradient stops increasing as illustrated by the flat portion of the graph in Fig. 5b. During this time, the viscosity jump at the inner interface increases while the viscosity jump at the outer interface decreases. The point at which the slope begins to increase again corresponds to the optimal viscosity profile over all values of average viscosity. After this point, the addition of polymer to increase viscosity would be detrimental to the stability of the system. To this point, only linear viscosity profiles have been considered, but many other viscosity profiles can also be used. The optimization procedure used for linear viscosity profiles in Fig. 3a has been repeated for several other types of viscosity profiles in Fig. 6. In addition to a viscosity profile which is linear at time |$\tau = 0$|⁠, viscosity profiles which are initially exponential and logarithmic are also considered. Recall that the viscosity profile changes with time in the |$r$|-coordinates. Therefore, it may be useful to consider viscosity profiles in the |$\zeta $|-coordinate system because they will remain fixed in time. Therefore, we also consider viscosity profiles which are linear and exponential with respect to |$\zeta $|⁠. A profile which is linear with respect to |$\zeta $| is quadratic with respect to |$r$| and a profile which is exponential in |$\zeta $| is proportional to |$e^{r^2}$|⁠. Figure 6a shows the dispersion relations of the optimal viscosity profiles of each type. The profile which is exponential in |$r$| is the most unstable while the profile which is exponential in |$\zeta $| is the least unstable. The corresponding optimal viscosity profiles are plotted in Fig. 6b. Notice that the profile which is least unstable, the one which is exponential in |$\zeta $|⁠, has both the smallest value of |$\mu (R_1)$| and the largest value of |$\mu (R_2)$|⁠. Therefore, it has the smallest viscosity jumps at the interfaces. Fig. 6. Open in new tabDownload slide (a) Plots of the dispersion relations for the optimal viscosity profiles, which are (i) linear with respect to |$r$|⁠, (ii) exponential with respect to |$r$|⁠, (iii) logarithmic with respect to |$r$|⁠, (iv) linear with respect to |$\zeta $| and (v) exponential with respect to |$\zeta $|⁠. (b) Plots of the corresponding viscosity profiles. Fig. 6. Open in new tabDownload slide (a) Plots of the dispersion relations for the optimal viscosity profiles, which are (i) linear with respect to |$r$|⁠, (ii) exponential with respect to |$r$|⁠, (iii) logarithmic with respect to |$r$|⁠, (iv) linear with respect to |$\zeta $| and (v) exponential with respect to |$\zeta $|⁠. (b) Plots of the corresponding viscosity profiles. 7.3 Time dependence In the previous sections, we considered the growth rate only at time |$\tau =0$|⁠. However, it is also important to understand how the growth rate changes with time. As time increases and the interfaces move outward, there are several physical factors at play. The curvature of the interfaces decreases which works to destabilize the flow while the velocity of the interfaces decreases which works to stabilize the flow. Gin & Daripa (2015) studied constant viscosity flows and found that the above two competing effects lead to non-monotonic behavior of the maximum growth rate. These results are in the original |$r$|-coordinate system. However, the maximum growth rate of disturbances in the |$\zeta $|-coordinate system is a monotonic function of time for constant viscosity flows. We demonstrate this analytically. For two-layer constant viscosity flows, the growth rate in the |$\zeta $|-coordinates (see equation (4.3)) is $$\begin{align*}& \sigma = \frac{Qn}{2 \pi R^2} \frac{\mu_o-\mu_i}{\mu_o+\mu_i} - \frac{T}{\mu_o+\mu_i} \frac{n^3-n}{R^3}, \end{align*}$$ where |$R(\tau )$| is the radius of the interface. Taking a derivative with respect to |$n$| and setting it equal to zero gives |$n_{max} = \sqrt{QR/(6 \pi T) (\mu _o-\mu _i) + 1/3}$|⁠. The most dangerous wave number is an integer value that is either the ceiling or floor of this number. Plugging this into (4.3) gives an upper bound on the maximum growth rate over all wave numbers. $$\begin{align*}& \begin{split} \sigma_{max} = &\frac{Q}{2 \pi R^2} \sqrt{\frac{QR}{6 \pi T}(\mu_o - \mu_i) + \frac{1}{3}} \left(\frac{\mu_o-\mu_i}{\mu_o+\mu_i}\right) \\ - &\frac{T}{\mu_o+\mu_i} \sqrt{\frac{QR}{6 \pi T}(\mu_o - \mu_i) + \frac{1}{3}} \left(\frac{QR}{6 \pi T}(\mu_o - \mu_i) - \frac{2}{3}\right) \frac{1}{R^3}. \end{split} \end{align*}$$ Taking the derivative with respect to |$R$| gives $$\begin{align*}& \frac{\partial \sigma_{max}}{\partial R} = - \frac{T \left(\frac{Q}{\pi T} (\mu_o-\mu_i)R + 2 \right) \left(\frac{Q}{\pi T} (\mu_o-\mu_i)R+4 \right)}{12 R^4 (\mu_i+\mu_o) n_{max}}. \end{align*}$$ If |$\mu _o> \mu _i$| then this expression is negative for all |$R$|⁠. Therefore, |$\sigma _{max}$| is a strictly decreasing function of |$R$| and hence time |$\tau $| in the |$\zeta $|-coordinates. For three-layer variable viscosity flows, there is an additional factor which affects the stability. The interfaces get closer together which makes the variable viscosity profile steeper and works to destabilize the flow. Despite this fact, the numerical results that follow show that |$\sigma _{max}$| is a decreasing function of time. This can be illustrated by the upper bound given by (5.9). The first two terms are strictly decreasing functions of |$R_1$| and |$R_2$| (and therefore of |$\tau $|⁠) while the third term is independent of time. Therefore, the upper bound is a decreasing function of time. In Fig. 7, the dispersion relation is plotted at several different times for a typical variable viscosity flow. Initially, the inner interface is at |$R_1 = 20$|⁠. As time increases, it moves outward. Note that the maximum value of |$\sigma $| decreases with time. However, the difference is more pronounced near the maximum value, which is mostly affected by the stability of the interfaces, than for short waves which are mostly affected by the layer instability. Fig. 7. Open in new tabDownload slide Several plots of the dispersion relation |$\sigma $| versus |$n$| at different times, as represented by the location of the inner interface. Parameter values are |$Q = 10$|⁠, |$\mu _i = 2$|⁠, |$\mu _o = 10$|⁠, |$\mu (R_1) = 5$|⁠, |$\mu (R_2) = 7$|⁠, |$T_1 = T_2 = 1$|⁠, |$R_1(0) = 20$| and |$R_2(0) = 30$|⁠. Fig. 7. Open in new tabDownload slide Several plots of the dispersion relation |$\sigma $| versus |$n$| at different times, as represented by the location of the inner interface. Parameter values are |$Q = 10$|⁠, |$\mu _i = 2$|⁠, |$\mu _o = 10$|⁠, |$\mu (R_1) = 5$|⁠, |$\mu (R_2) = 7$|⁠, |$T_1 = T_2 = 1$|⁠, |$R_1(0) = 20$| and |$R_2(0) = 30$|⁠. In order to shed more light on this time-dependent behavior, we investigate how |$\sigma _{max}$| and the most dangerous wave number |$n_{max}$| evolve in time. Figure 8a is a plot of |$\sigma _{max}$| versus |$R_1$|⁠, and Fig. 8b is a plot of |$n_{max}$| versus |$R_1$|⁠. Notice first that |$\sigma _{max}$| is a decreasing function of time and |$n_{max}$| is an increasing function of time. The fact that |$n_{max}$| increases with time is a well-known fact for constant viscosity flows (Cardoso & Woods, 1995). Also observe that there is a critical value |$R_1^*$| such that for |$R_1> R_1^*$|⁠, |$\sigma _{max}$| is constant and |$n_{max}$| is infinite. This is the point at which the layer instability comes to dominate the flow. For |$R_1 < R_1^*$|⁠, the instability of the interfaces dominates and the behavior or |$\sigma _{max}$| and |$n_{max}$| is similar to what happens for constant viscosity flow. For |$R_1> R_1^*$|⁠, the layer is more unstable than the interfaces, and therefore the short wave instability dominates and |$\sigma _{max} = {\lim _{n \to \infty } \sigma (n)}$|⁠. This behavior mimics what we see from the upper bound (5.9) in which the two terms related to the interfaces are decreasing functions while the term related to the layer instability is constant. Fig. 8. Open in new tabDownload slide Plots of (a) the maximum growth rate |$\sigma _{max}$| versus the radius of the inner interface |$R_1$| and (b) the most dangerous wave number |$n_{max}$| versus the radius of the inner interface |$R_1$|⁠. The parameter values are |$Q = 10$|⁠, |$\mu _i = 2$|⁠, |$\mu _o = 10$|⁠, |$\mu (R_1) = 5$|⁠, |$\mu (R_2) = 7$|⁠, |$T_1 = T_2 = 1$|⁠, |$R_1(0) = 20$| and |$R_2(0) = 30$|⁠. Fig. 8. Open in new tabDownload slide Plots of (a) the maximum growth rate |$\sigma _{max}$| versus the radius of the inner interface |$R_1$| and (b) the most dangerous wave number |$n_{max}$| versus the radius of the inner interface |$R_1$|⁠. The parameter values are |$Q = 10$|⁠, |$\mu _i = 2$|⁠, |$\mu _o = 10$|⁠, |$\mu (R_1) = 5$|⁠, |$\mu (R_2) = 7$|⁠, |$T_1 = T_2 = 1$|⁠, |$R_1(0) = 20$| and |$R_2(0) = 30$|⁠. 7.4 Validation of QSSA In Section 2, we invoke the QSSA which assumes that the basic solution changes slowly in comparison to the growth of the disturbances. We now examine the validity of that assumption. Consider a two-layer constant viscosity flow in which the growth rate of a disturbance with wave number |$n$| is given by equation (4.3). The maximum growth rate over all wave numbers can be written as $$\begin{align*}& \sigma_{max} = \frac{2T}{(\mu_o+\mu_i)R^3} \left( \frac{QR}{6 \pi T} (\mu_o-\mu_i) + \frac{1}{3} \right)^{\frac{3}{2}}. \end{align*}$$ The expression for the position of the interface is |$R(\tau ) = \sqrt{Q\tau /\pi + R(0)^2}$|⁠. In particular, |$R \propto \sqrt{\tau }$|⁠. Therefore, |$\sigma _{max} \propto \tau ^{-3/4}$| for |$\tau \gg 1$|⁠. In comparison, the interfacial position of the basic solution changes at a rate $$\begin{align*}& \frac{1}{R} \frac{dR}{d\tau} = \frac{Q}{2 \pi R^2} \propto \tau^{-1}. \end{align*}$$ Therefore, for large |$\tau $| the disturbances will grow faster than the basic solution. For three-layer variable viscosity flow, the QSSA has a more solid foundation. Recall from the previous subsection that the interfacial instability dominates at early times, but the instability of the middle layer dominates at later times. The upper bound (5.8) demonstrates that the layer instability is bounded by a constant term that depends only on |$\mu ^{\prime}(\zeta )$|⁠. This is further validated by the region of Fig. 8a in which |$\sigma _{max}$| is constant. Therefore, for variable viscosity flows the interfaces will be moving very slowly at later times while the growth of disturbances remains constant. Figure 9 is a numerical comparison of the growth rate of disturbances |$\sigma _{max}$| and the rates of change of each individual interface of the base flow. The parameters used are the same as in Fig. 8 and therefore the solid curve is the same as the curve in Fig. 8a. Notice that the growth rate of the disturbances is always greater than the rate of change of the interfaces, but that this is especially true for later times. The late time behavior will be true even if the interfaces are stabilized by very large interfacial tension. Fig. 9. Open in new tabDownload slide Plots of the maximum growth rate |$\sigma _{max}$| and the rates of change of the interfaces of the basic flow. All parameter values are the same as Fig. 8: |$Q = 10$|⁠, |$\mu _i = 2$|⁠, |$\mu _o = 10$|⁠, |$\mu (R_1) = 5$|⁠, |$\mu (R_2) = 7$|⁠, |$T_1 = T_2 = 1$|⁠, |$R_1(0) = 20$| and |$R_2(0) = 30$|⁠. Fig. 9. Open in new tabDownload slide Plots of the maximum growth rate |$\sigma _{max}$| and the rates of change of the interfaces of the basic flow. All parameter values are the same as Fig. 8: |$Q = 10$|⁠, |$\mu _i = 2$|⁠, |$\mu _o = 10$|⁠, |$\mu (R_1) = 5$|⁠, |$\mu (R_2) = 7$|⁠, |$T_1 = T_2 = 1$|⁠, |$R_1(0) = 20$| and |$R_2(0) = 30$|⁠. 7.5 Variable Injection Rate Recently Beeson-Jones & Woods (2015) and Gin & Daripa (2018) explored the idea of using a variable injection rate |$Q(t)$| to stabilize multi-layer constant viscosity flows. In these works, the maximum injection rate which results in a stable flow is investigated. Unfortunately, there is no injection rate which stabilizes a variable viscosity flow because short waves are always unstable. However, as an analogy, we can find the maximum injection rate that keeps the growth rate under a certain threshold. Figure 10 shows the maximum injection rate such that the maximum growth rate is below 0.001 for a certain constant viscosity flow and a certain variable viscosity flow. The constant viscosity flow has a viscosity of |$\mu = 6$| in the middle layer while the variable viscosity flow has a linear viscosity profile with |$\mu (R_1) = 5$| and |$\mu (R_2) = 7$|⁠. The variable viscosity flow allows for the fluid to be injected more quickly while maintaining the same level of instability. Fig. 10. Open in new tabDownload slide Plots of the maximum injection rate that results in a value of |$\sigma _{max} \leq 0.001$| for a constant viscosity flow with |$\mu = 6$| and a variable viscosity flow with |$\mu (R_1) = 5$| and |$\mu (R_2) = 7$|⁠. Other parameter values are |$\mu _i = 2$|⁠, |$\mu _o = 10$|⁠, |$T_1 = T_2 = 1$|⁠, |$R_1(0) = 20$| and |$R_2(0) = 30$|⁠. Fig. 10. Open in new tabDownload slide Plots of the maximum injection rate that results in a value of |$\sigma _{max} \leq 0.001$| for a constant viscosity flow with |$\mu = 6$| and a variable viscosity flow with |$\mu (R_1) = 5$| and |$\mu (R_2) = 7$|⁠. Other parameter values are |$\mu _i = 2$|⁠, |$\mu _o = 10$|⁠, |$T_1 = T_2 = 1$|⁠, |$R_1(0) = 20$| and |$R_2(0) = 30$|⁠. 8. Conclusions The stability of three-layer radial porous media flows with variable viscosity is an important issue in many applications. This work is the first to address this topic. First, the linear stability problem is formulated resulting in an eigenvalue problem with time-dependent coefficients and eigenvalue-dependent boundary conditions. This derivation requires an appropriate change of variables that fixes the positions of the interfaces and the viscosity profile of the middle layer fluid. Several important analytical aspects of the eigenvalue problem are studied. First, upper bounds on the spectrum are derived using a variational approach. Then it is shown that for a certain bandwidth of wave numbers, there is a countably infinite set of positive eigenvalues, and the corresponding eigenfunctions are complete in a certain |$L^2$| space. The eigenvalues have been computed numerically in order to investigate the effect of various parameters on the stability of variable viscosity flows. The following are some of the key numerical results: (i) Variable viscosity flows can reduce the maximum growth rate by reducing the jumps in viscosity at the interfaces, but short waves become unstable; (ii) Short wave instability is dominated by the viscosity gradient in the layer while long and intermediate wavelengths are dominated by the instability of the interfaces; (iii) The optimal viscosity profile is the one which balances the interfacial instability with the instability of the layer; (iv) increasing interfacial tension decreases the viscosity gradient of the optimal viscosity profile; (v) A viscosity profile which is exponential with respect to |$\zeta $| is optimal among the types of profiles considered; (vi) |$\sigma _{max}$| is a decreasing function of time. This is due to the instability of the interfaces decreasing with time while the layer instability remains relatively unchanged; and (vii) Variable viscosity flows allow for faster injection without making the flow more unstable. There are several directions for future work in this area. Theorem 6.2 gives a way to numerically simulate the motion of the interfaces based on the linear theory. Any initial perturbation of the interfaces and viscous profile can be expanded in terms of the complete set of orthonormal eigenfunctions using equation (6.9). This could then be compared with full nonlinear simulations and experiments of multilayer radial Hele-Shaw flows with variable viscosity. Another direction is to expand the scope of Theorem 6.2 to all wave numbers and parameter regimes. Finally, a weakly nonlinear analysis, such as Miranda & Widom (1998) for two-layer flows with constant viscosity, or fully nonlinear stability analysis could be performed on this system. Acknowledgements This work has been supported in part by the U.S. National Science Foundation grant DMS-1522782. Additionally, the authors would like to acknowledge the reviewers whose comments have helped improve the paper immensely. References Anjos , P. , Dias , E. & Miranda , J. ( 2015 ) Kinetic undercooling in Hele-Shaw flows . Phys. Rev. E (3) , 92 , 043019. Google Scholar OpenURL Placeholder Text WorldCat Beeson-Jones , T. & Woods , A. ( 2015 ) On the selection of viscosity to suppress the Saffman–Taylor instability in a radially spreading annulus . J. Fluid Mech. , 782 , 127 – 143 . Google Scholar Crossref Search ADS WorldCat Beeson-Jones , T. & Woods , A. ( 2017 ) Control of viscous instability by variation of injection rate in a fluid with time-dependent rheology . J. Fluid Mech. , 829 , 214 – 235 . Google Scholar Crossref Search ADS WorldCat Cardoso , S. & Woods , A. ( 1995 ) The formation of drops through viscous instability . J. Fluid Mech. , 289 , 351 – 378 . Google Scholar Crossref Search ADS WorldCat Dallaston , M. & McCue , S. ( 2013 ) Bubble extinction in Hele-Shaw flow with surface tension and kinetic undercooling regularization . Nonlinearity , 26 , 1639 – 1665 . Google Scholar Crossref Search ADS WorldCat Daripa , P. ( 2008 ) Hydrodynamic stability of multi-layer Hele-Shaw flows . J. Stat. Mech. Theory Exp. , P12005 . Google Scholar OpenURL Placeholder Text WorldCat Daripa , P. & Ding , X. ( 2012 ) A numerical study of instability control for the design of an optimal policy of enhanced oil recovery by tertiary dispalcement processes . Transp. Porous Media , 93 , 675 – 703 . Google Scholar Crossref Search ADS WorldCat Daripa , P. , Glimm , J., Lindquist , B. & McBryan , O. ( 1988 ) Polymer floods: a case study of nonlinear wave analysis and of instability control in tertiary oil recovery . SIAM J. Appl. Math. , 48 , 353 – 373 . Google Scholar Crossref Search ADS WorldCat Daripa , P. & Pasa , G. ( 2006 ) A simple derivation of an upper bound in the presence of a viscosity gradient in three-layer Hele-Shaw flows . J. Stat. Mech. Theory Exp. , P01014 . Google Scholar OpenURL Placeholder Text WorldCat Dias , E. & Miranda , J. ( 2010 ) Control of radial fingering patterns: a weakly nonlinear approach . Phys. Rev. E , 81 , 016312. Google Scholar OpenURL Placeholder Text WorldCat Gin , C. & Daripa , P. ( 2015 ) Stability results for multi-layer radial Hele-Shaw and porous media flows . Phys. Fluids , 27 , 012101. Google Scholar OpenURL Placeholder Text WorldCat Gin , C. ( 2015 ) Topics in stability analysis of multi-layer Hele-Shaw and porous media flows . Doctoral Dissertation . Gin , C. & Daripa , P. ( 2018 ) Time-dependent injection strategies and interfacial stability in multi-layer Hele-Shaw and porous media flows. arXiv e-prints, arXiv:1811.10721. Gorell , S. & Homsy , G. ( 1983 ) A theory of the optimal policy of oil recovery by the secondary displacement process . SIAM J. Appl. Math. , 43 , 79 – 98 . Google Scholar Crossref Search ADS WorldCat Ince , E. ( 1956 ) Ordinary Differential Equations . New York : Dover Publications . Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Kim , H. , Funada , T., Joseph , D. & Homsy , G. ( 2009 ) Viscous potential flow analysis of radial fingering in a Hele-Shaw cell . Phys. Fluids , 21 , 074106. Google Scholar OpenURL Placeholder Text WorldCat Kim , J. , Xu , F. & Lee , S. ( 2017 ) Formation and destabilization of the particle band on the fluid–fluid interface . Phys. Rev. Lett. , 118 , 074501. Google Scholar OpenURL Placeholder Text WorldCat Li , S. , Lowengrub , J., Fontana , J. & Palffy-Muhoray , P. ( 2009 ) Control of viscous fingering patterns in a radial Hele-Shaw cell . Phys. Rev. Lett. , 102 , 174501. Google Scholar OpenURL Placeholder Text WorldCat Luo , R. , Chen , Y. & Lee , S. ( 2018 ) Particle-induced viscous fingering: Review and outlook . Phys. Rev. Fluids , 3 , 110502. Google Scholar OpenURL Placeholder Text WorldCat Miranda , J. & Widom , M. ( 1998 ) Radial fingering in a Hele-Shaw cell: a weakly nonlinear analysis . Phys. D , 120 , 315 – 328 . Google Scholar Crossref Search ADS WorldCat Morrow , L. , Moroney , T. & McCue , S. ( 2019 ) Numerical investigation of controlling interfacial instabilities in non-standard Hele-Shaw configurations . J. Fluid Mech. , 877 , 1063 – 1097 . Google Scholar Crossref Search ADS WorldCat Muskat , M. ( 1934 ) Two fluid systems in porous media. The encroachment of water into an oil sand . Phys. Dokl. , 5 , 250 – 264 . Google Scholar OpenURL Placeholder Text WorldCat Muskat , M. ( 1946 ) The Flow of Homogeneous Fluids Through Porous Media . New York : McGraw Hill . Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Paterson , L. ( 1981 ) Radial fingering in a Hele–Shaw cell . J. Fluid Mech. , 113 , 513 – 529 . Google Scholar Crossref Search ADS WorldCat Saffman , P. & Taylor , G. ( 1958 ) The penetration of a fluid into a porous medium or Hele-Shaw cell containing a more viscous liquid . Proc. R. Soc. Lond. Ser. A , 245 , 312 – 329 . Google Scholar OpenURL Placeholder Text WorldCat Tang , H. , Grivas , W., Homentcovschi , D., Geer , J. & Singler , T. ( 2000 ) Stability considerations associated with the meniscoid particle band at advancing interfaces in Hele-Shaw suspension flows . Phys.Rev. Lett. , 85 , 2112 . Google Scholar Crossref Search ADS PubMed WorldCat Vaquero-Stainer , C. , Heil , M., Juel , A. & Pihler-Puzović , D. ( 2019 ) Self-similar and disordered front propagation in a radial Hele-Shaw channel with time-varying cell depth . Phys. Rev. Fluids , 4 . doi: 064002. Google Scholar OpenURL Placeholder Text WorldCat Walter , J . ( 1973 ) Regular eigenvalue problems with eigenvalue parameter in the boundary condition . Math. Z. , 133 , 301 – 312 . Google Scholar Crossref Search ADS WorldCat Xu , F. & Lee , S. ( 2019 ) The enhancement of viscous fingering with bidisperse particle suspension . J. Fluid Mech. , 860 , 487 – 509 . Google Scholar Crossref Search ADS WorldCat Zheng , Z. , Kim , H. & Stone , H. ( 2015 ) Controlling viscous fingering using time-dependent strategies . Phys. Rev. Lett. , 115 , 174501. Google Scholar OpenURL Placeholder Text WorldCat © The Author(s) 2021. Published by Oxford University Press on behalf of the Institute of Mathematics and its Applications. All rights reserved. This article is published and distributed under the terms of the Oxford University Press, Standard Journals Publication Model (https://academic.oup.com/journals/pages/open_access/funder_policies/chorus/standard_publication_model)
TI - Stability results on radial porous media and Hele-Shaw flows with variable viscosity between two moving interfaces
JF - IMA Journal of Applied Mathematics
DO - 10.1093/imamat/hxab001
DA - 2021-03-18
UR - https://www.deepdyve.com/lp/oxford-university-press/stability-results-on-radial-porous-media-and-hele-shaw-flows-with-Cy59vFcldd
SP - 1
EP - 1
VL - Advance Article
IS - 
DP - DeepDyve
ER -