FLUID INJECTION IN A POROUS MEDIUM: THE RADIAL SAFFMAN–TAYLOR INSTABILITY

Abstract

We consider planar flow involving two viscous fluids in a porous medium. One fluid is injected through a line source at the origin and moves radially outwards, pushing the second, ambient fluid outwards. There is an interface between the two fluids and if the inner injected fluid is of lower viscosity, the interface is unstable to small disturbances and radially directed unstable Saffman–Taylor fingers are produced. A linearized theory is presented and is compared with nonlinear results obtained using a numerical spectral method. An additional theory is also discussed, in which the sharp interface is replaced with a narrow diffuse interfacial region. We show that the nonlinear results are in close agreement with the linearized theory for small-amplitude disturbances at early times, but that large-amplitude fingers develop at later times and can even detach completely from the initial injection region.

1. Introduction

Suppose two different viscous fluids lie horizontally in a porous medium. There is a sharp interface between the two fluids and it necessarily forms a flat horizontal plane. The lower fluid has viscosity $\mu _1$ and the upper has viscosity $\mu _2$ . Now the lower fluid is caused to move vertically upwards with some constant speed, so that the top fluid 2 is effectively being pushed upwards by bottom fluid 1. In a classical paper, Saffman and Taylor [25] carried out a linearized stability analysis of this problem and showed that the horizontal interface between the two fluids would be unstable if $\mu _1 < \mu _2$ . Thus, when a less viscous fluid pushes a more viscous one, small disturbances at the interface are unstable and grow exponentially fast in time, eventually forming finger-like structures at the interface. Saffman and Taylor also investigated flow near the nose of one of these fingers and confirmed their calculations by comparison with experiment. Since their analysis made use of linearized theory, it is therefore predicated on the assumption that perturbations to the initial interface shape remain small as time progresses. Clearly this cannot remain true for unstable interfaces with rapidly growing fingers and so there must be some time beyond which linearized theory ceases to be valid. At later times, therefore, nonlinear effects will become evident at the interface and its shape would need to be computed numerically. An early attempt to do this was presented by Davidson [5]. Much has subsequently been written about the Saffman–Taylor instability, and good reviews of the subject are given by Homsy [14] and Couder [4].

The Saffman–Taylor instability, at least from the mathematical point of view, is similar to several other unstable flow types in which an interface separates two fluids with different material properties. Notable among these are Rayleigh–Taylor flows, in which one fluid of high density lies above a second fluid of lower density. The interface between them is again unstable and small perturbations to it rapidly develop into fingers of the heavy fluid moving down, alternating with bubbles of the lighter fluid floating upwards. Rayleigh–Taylor flows can form over vastly differing length and time scales, ranging from the microscopic to galactic (see Kelley et al. [17]). As a result, those types of interfacial instability have been the focus of intense research. Rayleigh–Taylor flows are also possible in porous media, where the fluid density may depend on factors such as temperature or the concentration of dissolved salts. A model of this type has been presented by Trevelyan et al. [28] who considered density to depend on the concentrations of two different solutes. More recently, Forbes et al. [9] compared two nonlinear models for Rayleigh–Taylor flow in a porous medium and demonstrated that including diffusion at a diffuse interface, with density dependent on the concentration of a single solute, allowed very large amplitude unstable fingers to form. Unlike the situation for regular fluids in contact at an interface, the fingers formed in a porous medium did not possess overturning sections at their heads. Because the Saffman–Taylor flow also involves fluid percolating through a porous medium, it is therefore to be expected that the way in which fingers evolve at its interface should be similar to those seen by Forbes et al. [9].

Saffman–Taylor flows are of interest in a number of practical situations. Experimentally, they have been studied using the Hele–Shaw apparatus, which consists of a viscous fluid confined between closely spaced parallel plates. This device enables laboratory visualization of the long fingers that form at the interface when a less-viscous fluid is forced into the apparatus that already contains an ambient fluid of higher viscosity. Some images of finger evolution are displayed in the article by Paterson [23], who found large-amplitude patterns when fluid was either injected or withdrawn radially in a Hele–Shaw cell. Instability and splitting of Saffman–Taylor fingers in a Hele–Shaw cell are demonstrated in the experimental work of Tabeling et al. [27]. Recently, Islam and Gandhi [15] have suggested that the fractal-like branching networks found in biological systems and respiratory networks can be fabricated in the laboratory by controlling spatial anisotropy in a Saffman–Taylor experiment.

The original paper by Saffman and Taylor [25] considered an initially planar interface between two horizontal fluid layers. Here, however, we instead focus on a radial outflow, in which an inner fluid is injected through a line source placed at the origin of some coordinate system. The flow through the porous medium is supposed to be two-dimensional, similar to conditions in a Hele–Shaw cell [23], and the injected fluid pushes the surrounding fluid outwards. Based on the original work of Saffman and Taylor [25], it is to be expected that the circular interface between the two fluids in radial flow should also be unstable if $\mu _1 < \mu _2$ , that is, if the inner injected fluid is less viscous than the outer ambient fluid. The interface will become progressively more distorted with time, forming large-amplitude radially directed waves around its circumference. Eventually, nonlinear effects will come to dominate the characteristics of the outflow and the shapes of all the fingers at the interface. Miranda and Widom [21] carried out a weakly nonlinear analysis of radial Saffman–Taylor flow, focussing on the second-order terms in the expression for the interface shape. Their calculations show that nonlinearity is responsible for three important qualitative changes to the predictions of linearized theory alone. First, nonlinear effects allow coupling between different Fourier modes, with the result that the fingers formed at the interface are broader than linearization suggests, since superharmonic and subharmonic modes are excited. Second, as a consequence of broadening of the fingers, splitting is to be expected at the fingertips. A third consequence of nonlinearity (evident from their Figure 4) is that regions of high curvature form at selected points of the interface profile and this will be discussed more fully in Section 6.

Experimental observations confirm the first two of these effects of nonlinearity; that is, as time progresses, the unstable fingers that form on the expanding interface broaden as they develop and fingertip splitting then takes place. These effects are particularly evident in the images presented by Chen [3], who injected either dyed water or dyed oil into a Hele–Shaw cell, displacing ambient glycerine. After some time, the injected fluid formed intricate patterns, with an interface that was almost fractal in appearance, caused by repeated widening and splitting of the fingers at smaller and smaller length scales.

The third consequence of nonlinearity, which we have identified above, relates to the formation of regions of very high curvature at the interface, when the two fluids are immiscible and a genuine sharp interface exists between them. Sorbie et al. [26] refer to the difficulty experienced by researchers aiming to compute immiscible viscous fingering by direct numerical solution of the governing partial differential equations. They stress that they do not see this difficulty as arising from issues relating purely to numerical accuracy or numerical instability, but rather that a deeper problem with the underlying physics is responsible. They were concerned to create a practical method for modelling oil recovery in situations where oil trapped underground in a porous medium was driven towards an extraction well using injected water. This creates a Saffman–Taylor scenario with an elaborately fingered interface, since the water has lower viscosity than the oil. They presented a semi-heuristic scheme that allowed them to compute quantities of interest with two immiscible fluids present.

This difficulty with computing the behaviour of a sharp interface between immiscible fluids is well known in other, similar, fluid instabilities. Moore [22] first showed that a sharp interface between two different fluids can develop regions of very high curvature and that, in fact, the curvature becomes infinite within some finite critical time $t = t_C$ . His analysis was carried out for planar Kelvin–Helmholtz instability, which occurs when two inviscid fluids move at different speeds either side of a sharp interface. Moore used asymptotic analysis to estimate the critical time $t_C$ . His formula indicates that $t_C$ decreases markedly as the amplitude of the initial disturbance at the interface increases. Forbes [7] demonstrated that Moore’s [22] formula also gave a good estimate for the critical time $t_C$ , at which a curvature singularity forms at the interface in the planar Rayleigh–Taylor instability. Forbes [8] later computed planar radial Rayleigh–Taylor outflow with less-dense fluid injected through a line source into a more dense ambient fluid, in the analogous situation to the viscous Saffman–Taylor instability discussed here. He showed that curvature singularities also formed at the interface within finite time, when the two fluids did not mix and a sharp interface was present. The physical model of the flow ceases to be valid beyond this critical time, $t> t_C$ . Nevertheless, by changing the physical model to one in which the two fluids were weakly miscible and replacing the sharp interface with a narrow interfacial region across which the fluid density could change rapidly but smoothly, Forbes [8] demonstrated that the miscible-fluid model could now continue beyond the critical time $t> t_C$ . Radially directed fingers formed around the interface and, after sufficient time, they formed overturning mushroom-shaped structures at their heads.

It is therefore apparent that purely immiscible models of Saffman–Taylor instability will not be appropriate for large times and that miscible models will be required. Pritchard [24] carried out a linearized analysis for radial injection into a Hele–Shaw cell, in the case when the viscosity of the fluid system varied continuously with both temperature and the concentration of some solute; an explicit interface was not included in this analysis and so the system is physically equivalent to a miscible one in which there are finite-width interfacial zones. He found that there were two moving fronts in his model, one where temperature changed rapidly from one state to another and a second across which the concentration of the solute changed. This gave two effective interfaces, either of which could exhibit unstable behaviour. Iwasaki et al. [16] demonstrated experimentally that viscous finger formation in radial Saffman–Taylor flow could be suppressed by using a system of two fluids in which partial mixing across the interfacial zone could occur. Viscous fingering was also shown by Hinton and Jyoti [13] to be suppressed in a mathematical model in which the two viscous fluids were allowed to form vertical buoyant plumes within a porous medium in three-dimensional geometry.

In this present paper, we review in Section 2 the governing equations for a system of two fluids in a porous medium in which a sharp interface is present. We develop a linearized solution in Section 3 for radial Saffman–Taylor flow driven by a line source at the origin. The base flow consists of a simple radial outflow with a circular interface that moves outwards as time progresses and linearization consists of small-amplitude perturbations to this flow. The linearized stability condition is presented in Section 3.1. Linearization can be expected to give a good approximation to the Saffman–Taylor flow so long as perturbations to the simple base flow remain small, but, as discussed above, perturbations grow exponentially rapidly in unstable cases, so that linearization soon loses validity and nonlinear effects then dominate. Accordingly, we present in Section 4 a semi-numerical method for solving the fully nonlinear equations appropriate to the case of a sharp interface between two immiscible fluids. Large-amplitude perturbations are allowed for in this model, but it is found that the numerical method fails within finite time; as discussed above and consistent with Sorbie et al. [26], we suggest this is a physical limitation in the immiscible model, rather than a numerical issue, and we argue that curvature singularity at the sharp interface is responsible. For this reason, we consider in Section 5 a third model in which the two fluids are now miscible. This description of the problem treats the system as a single fluid in which the viscosity varies continuously throughout the region. Fluid viscosity is supposed to vary linearly with temperature and an additional convection–diffusion equation for temperature is introduced into the problem description. There is therefore no longer a sharp interface, but instead a narrow interfacial zone exists in the fluid, across which the viscosity (and temperature) vary smoothly but rapidly from one value to another. Results of these three models are presented and discussed in Section 6, with each successive model having effectively evolved from the limitations of its prior. A discussion concludes the paper in Section 7.

2. Governing equations for sharp-interface model

This section reviews the basic equations needed to describe the fully nonlinear problem for the case in which two completely immiscible fluids meet at an infinitesimally thin interface, described by some curve $r = R (\theta ,t)$ in plane polar coordinates. The system under consideration involves a porous material between two parallel plates separated by a narrow gap of thickness $2h$ . The system therefore behaves as a porous Hele–Shaw cell. A hole is present in the upper plate at the origin of a coordinate system and one viscous fluid is injected through it into the medium, which is already suffused with a different viscous fluid. The flow is assumed to be steady and incompressible and the two fluids meet at an initially circular interface of radius a. Since the thickness of the cell is small and the cell is oriented horizontally, the gravitational force on the fluid is considered to be insignificant. We will assume that the fluid is slow moving and dominated by viscous friction. As a result, the fluid seepage velocity $\mathbf {q}$ is related to the pore pressure p through Darcy’s law,

(2.1) $$ \begin{align} \mathbf{q} = -\frac{k}{\mu}\, \nabla p , \end{align} $$

where k is the porosity of the medium (in this set-up: $k = h^2/3 $ [m²]). The quantity $\mu $ represents the viscosity of the fluid and p denotes the pore pressure of the fluid within the porous medium. The usual plane polar coordinates $(r,\theta )$ will be used, and these are related to Cartesian coordinates $(x,y)$ by means of the familiar relations $x=r\cos \theta $ and $y=r\sin \theta $ .

From this point onwards, all variables will be nondimensionalized. Let us denote the inner (injected) fluid as fluid 1 and the outer (displaced) one as fluid 2, with respective viscosities $\mu _1$ and $\mu _2$ . At the origin $r = 0$ , fluid 1 is injected at volume rate Q and this process is modelled using a simple line source. At the initial time $t = 0$ , fluid 1 occupies a cylinder of radius a, which will be taken to be the length scale in the new nondimensional coordinates. All times are referenced to the quantity $(2 \pi h a^2)/Q$ and the appropriate scale for speed is $Q / (2 \pi h a)$ . Pressures are nondimensionalized using the quantity $(Q \mu _2)/(2 \pi h k)$ . In these dimensionless coordinates, the initial cylinder containing inner fluid 1 has radius $1$ and an important parameter is

(2.2) $$ \begin{align} \beta = \mu_1 / \mu_2 , \end{align} $$

which is the dimensionless ratio of the viscosities of the inner and outer fluids. A definition sketch of this nondimensional system is given in Figure 1.

Figure 1 Schematic diagram of the dimensionless Saffman–Taylor system with radial outflow. A line source is present at the origin and injects fluid with unit flux. The unperturbed initial interface is the unit circle (sketched with a dashed line, red online). The location of the (sharp) interface is $r = R(\theta ,t)$ and $\beta $ represents the viscosity ratio.

In dimensionless coordinates and variables, Darcy’s law (2.1) for inner fluid 1 takes the form

(2.3) $$ \begin{align} \mathbf{q}_1 = - \frac{1}{\beta} \, \nabla p_1\quad\textrm{in } 0 < r < R(\theta ,t), \end{align} $$

and in outer fluid 2, the governing law is

(2.4) $$ \begin{align} \mathbf{q}_2 = - \nabla \, p_2\quad \textrm{in } r> R(\theta ,t). \end{align} $$

The two fluids each satisfy the mass-continuity equation

(2.5) $$ \begin{align} \nabla \cdot \mathbf{q}_j = 0 \end{align} $$

for $ j = 1,2$ , in their appropriate domains. Since a line source is assumed to be present at the origin, it follows that the fluid velocity there takes the form

(2.6) $$ \begin{align} \mathbf{q}_1 \rightarrow \frac{1}{2r} \, \mathbf{e_r} \quad \textrm{as } r \rightarrow 0. \end{align} $$

In this expression, $\mathbf {e_r}$ represents the unit vector pointing radially outwards from the origin. In each fluid, it follows from conservation of mass (2.5) and Darcy’s law (2.3) or (2.4) that the pressure within each fluid satisfies Laplace’s equation

(2.7) $$ \begin{align} \nabla^2 p_j = 0 \quad\text{for } j = 1,2. \end{align} $$

The interface, located at $r = R(\theta ,t)$ , is subject to three boundary conditions. Two of these are kinematic, expressing the fact that neither fluid is free to cross the interface and enter into the domain of the other fluid. This gives rise to conditions

(2.8) $$ \begin{align} u_j = \frac{\partial R}{\partial t} + \frac{v_j}{R}\frac{\partial R}{\partial \theta} \quad \textrm{on } r = R(\theta,t) ,\quad j=1,2. \end{align} $$

There is also a dynamic condition at the interface and it indicates that the jump in pressure across that surface is due to the surface tension there. This is written as

(2.9) $$ \begin{align} p_2 - p_1 = - \sigma \, \kappa\quad \textrm{on } r = R(\theta, t). \end{align} $$

The dimensionless constant

(2.10) $$ \begin{align} \sigma = \frac{2 \pi h k \gamma}{\mu_2 a Q} \end{align} $$

is a surface-tension parameter and the function

(2.11) $$ \begin{align} \kappa = \frac{R^2 + 2R_\theta^2 - R R_{\theta\theta}}{[R^2 + R^2 _\theta]^{3/2}} \end{align} $$

is the curvature of the interface.

3. Linearized solution

To approximate the behaviour of this nonlinear system with the solution to a simpler linearized problem, we assume that the system behaves as a small perturbation to a known, basic outflow. The obvious background flow in this situation is one in which the velocity in both fluid layers is simply the outflow (2.6) generated by the line source and presumed to hold everywhere. To complement this velocity vector, we assume that the background interface shape remains circular, with expanding radius $R_0 (t)$ . Thus,

(3.1) $$ \begin{align} u_j^{(0)}(r) = \frac{1}{2r};\quad v_j^{(0)}(r) = 0,\quad j=1,2, \end{align} $$

and the kinematic condition (2.8) reduces to

(3.2) $$ \begin{align} \frac{\partial R_0}{\partial t} = \frac{1}{2 R_0}. \end{align} $$

This simple equation (3.2) is integrated, under the assumption that the initial bubble radius was $R_0(0) = 1$ , and shows that the injection source causes the background bubble radius to increase with time, according to

(3.3) $$ \begin{align} R_0(t) = \sqrt{1 + t}. \end{align} $$

The pore pressures associated with this background flow are obtained by substituting the velocity components $u_j^{(0)}$ and $v_j^{(0)}$ from (3.1) into the two dimensionless forms (2.3) and (2.4) of Darcy’s law. This gives

(3.4) $$ \begin{align} \begin{aligned} p_1^{(0)}(r,t) & = -\frac{\beta}{2} \ln(r) + f_1(t) \quad \textrm{in } 0 < r < R_0(t) , \\ p_2^{(0)}(r,t) & = -\frac{1}{2} \ln(r) + f_2(t)\quad \textrm{in } r> R_0(t). \end{aligned} \end{align} $$

Here, the two quantities $f_1(t)$ and $f_2(t)$ are arbitrary functions of integration.

We suppose that some small perturbation is now made to this underlying background flow (3.1), (3.3), and is characterized by a small parameter $\epsilon $ . This constant $\epsilon $ typically arises from a description of initial disturbances to the background flow. The solution variables are expanded in series form to give

(3.5) $$ \begin{align} \begin{aligned} p_j(r,\theta,t) & = p_j^{(0)}(r,t) + \epsilon p_j^{(1)}(r,\theta,t) + \mathcal{O} ( \epsilon^2) \\ u_j(r,\theta,t) & = u_j^{(0)}(r,t) + \epsilon u_j^{(1)}(r,\theta,t) + \mathcal{O} ( \epsilon^2)\quad\textrm{for } j=1,2. \end{aligned} \end{align} $$

Here, the zeroth-order pore pressures and radial velocity components are given in (3.4) and (3.1), respectively. The radius of the perturbed interface $r = R(\theta ,t)$ is likewise written

(3.6) $$ \begin{align} R(\theta, t) = R_0(t) + \epsilon R_1(\theta,t) + \mathcal{O} ( \epsilon^2 ) , \end{align} $$

in which $R_0(t)$ is the known base radius (3.3) of the expanding axisymmetric background flow.

The perturbation expansions (3.5), (3.6) are substituted into the governing equations (2.7) and boundary conditions (2.8), (2.9). Furthermore, since the boundary conditions hold at the exact (unknown) interface $r = R(\theta , t)$ , Taylor series expansions are used to project these conditions onto the known location $r = R_0 (t)$ given in (3.3). Correct to the first order in the small parameter $\epsilon $ , the two kinematic boundary conditions (2.8) linearize to

(3.7) $$ \begin{align} u_j^{(1)} = \frac{\partial R_1}{\partial t} + \frac{R_1}{4 R_0^2} \quad \textrm{on } r = R_0(t) , \quad j=1,2. \end{align} $$

At zeroth order in $\epsilon $ , the dynamic condition (2.9) yields simply

(3.8) $$ \begin{align} f_2(t) - f_1(t) = - \frac{\sigma}{R_0(t)} , \end{align} $$

involving the two functions of integration $f_1$ and $f_2$ in (3.4). At first order, the linearized dynamic condition becomes

(3.9) $$ \begin{align} [ p_2^{(1)} - p_1^{(1)}] + R_1 \bigg[ \frac{\beta - 1}{2 R_0} \bigg] = \frac{\sigma}{R_0^2} \bigg[ R_1 + \frac{\partial^2 R_1}{\partial \theta^2} \bigg] \quad \textrm{on } r = R_0(t). \end{align} $$

The linearized pressure function $p_1^{(1)}$ satisfies Laplace’s equation (2.7) in the linearized interior region $0 < r < R_0 (t)$ and pressure function $p_2^{(1)}$ obeys the same equation in $r> R_0 (t)$ .

In view of the expansions (3.5), the desired forms of the linearized pressure functions that satisfy Laplace’s equation in their respective domains are

(3.10) $$ \begin{align} \begin{aligned} p_1^{(1)}(r, \theta ,t) & = \sum_{n=1}^{\infty} r^n [ A_n^{*}(t) \cos(n\theta) + B_n^{*}(t) \sin(n\theta)] \quad \textrm{in } 0 < r < R_0(t) , \\ p_2^{(1)}(r, \theta ,t) & = \sum_{n=1}^{\infty} r^{- n} [ C_n^{*}(t) \cos(n\theta) + D_n^{*}(t) \sin(n\theta)] \quad \textrm{in } r> R_0(t). \end{aligned} \end{align} $$

The form of the perturbed interface $R_1(\theta , t)$ is chosen to be consistent with (3.10) and so is

(3.11) $$ \begin{align} R_1(\theta, t) = \sum_{m=1}^{\infty} M_n^{*}(t) \cos(n\theta) + N_n^{*}(t) \sin(n\theta). \end{align} $$

In these expressions, the functions $A_n^{*}(t)$ and so on are so far arbitrary. The velocity components in each of the two regions are found by differentiating the pressure functions in the series representations (3.10), using Darcy’s law (2.3) in fluid 1 and (2.4) in outer fluid 2.

The linearized boundary conditions (3.7)–(3.9) determine the Fourier coefficients in the representations (3.10) and (3.11). The first kinematic condition in (3.7), with $j=1$ , gives the two sets of differential equations

(3.12) $$ \begin{align} \begin{aligned} \frac{{d} M_n^{*}(t)}{{d} t} & = - \frac{1}{\beta} n R_0^{n-1} A_n^{*}(t) - \frac{1}{4 R_0^2} M_n^{*}(t) , \\ \frac{{d} N_n^{*}(t)}{{d} t} & =- \frac{1}{\beta} n R_0^{n-1} B_n^{*}(t) - \frac{1}{4 R_0^2} N_n^{*}(t). \end{aligned} \end{align} $$

A similar system of differential equations can be derived in the same way from the second kinematic boundary condition in (3.7) with $j=2$ . However, it is easier to consider the difference between the two kinematic conditions, which at once gives $u_2^{(1)} = u_1^{(1)}$ on $r = R_0(t)$ . This results in the further two sets of algebraic equations

(3.13) $$ \begin{align} \begin{aligned} C_n^{*}(t) & = - \frac{1}{\beta} R_0^{2n} A_n^{*}(t) , \\ D_n^{*}(t) & = - \frac{1}{\beta} R_0^{2n} B_n^{*}(t). \end{aligned} \end{align} $$

Finally, the linearized dynamic condition (3.9) gives two more systems of algebraic equations for the Fourier coefficients. These are

(3.14) $$ \begin{align} \begin{aligned} C_n^{*}(t) R_0^{-n} - A_n^{*}(t) R_0^n - M_n(t) \frac{1 - \beta}{2} R_0^{-1} & = \frac{\sigma}{R_0^2} ( 1-n^2) M_n^{*}(t) , \\ D_n^{*}(t) R_0^{-n} - B_n^{*}(t) R_0^n - N_n(t) \frac{1 - \beta}{2} R_0^{-1} & = \frac{\sigma}{R_0^2}( 1-n^2) N_n^{*}(t). \end{aligned} \end{align} $$

To solve this linearized problem, the systems of equations (3.12)–(3.14) must be solved for the time-dependent Fourier coefficients $A_n^{*}$ , $B_n^{*}$ , $C_n^{*}$ and $D_n^{*}$ in the series (3.5), along with coefficients $M_n^{*}$ and $N_n^{*}$ in the representation (3.6) of the interface shape. Taking the first equation in each of (3.12)–(3.14) leads to a first-order differential equation for the coefficients $M_n^{*}$ . We obtain

$$ \begin{align*} \frac{d M_n^{*}}{dt} + \frac{M_n^{*}}{R_0^2} \bigg[ \frac{1}{4} + \frac{n(\beta-1)}{2(\beta+1)} \bigg] + \frac{M_n^{*}}{R_0^3} \frac{\sigma n ( n^2 -1)}{\beta+1} = 0. \end{align*} $$

An identical equation is also obtained for the coefficients $N_n^{*}$ from the second equation in each of the three sets of boundary conditions (3.12)–(3.14). This differential equation is simplified very significantly by regarding the base radius $R_0(t)$ in (3.3) as the independent variable, rather than the time t itself. Application of the chain rule of calculus then leads to

(3.15) $$ \begin{align} \frac{{d} M_n^{*}}{{d} R_0} = M_n^{*} \bigg( \frac{\lambda_n}{R_0} - \frac{\delta_n}{R_0^2} \bigg). \end{align} $$

Now, the solution to (3.15) can be obtained at once, in the form

(3.16) $$ \begin{align} M_n^{*} ( R_0) = \mathcal{C}_n R_0^{\lambda_n} {e}^{\delta_n/R_0}. \end{align} $$

In this expression, the quantities $\mathcal {C}_n$ are constants of integration and we have defined additional constants

(3.17) $$ \begin{align} \begin{aligned} \lambda_n & = - \frac{1}{2} + \frac{n(1-\beta)}{1 + \beta} , \\ \delta_n & = \frac{2 \sigma n(n^2-1)}{1 + \beta} , \end{aligned} \end{align} $$

for convenience. A similar solution to (3.16) is likewise obtained for the coefficients $N_n^{*},$ so that the linearized solution for the interface shape in (3.11) becomes

(3.18) $$ \begin{align} R_1(\theta, t) = \sum_{m=1}^{\infty} R_0^{\lambda_n} {e}^{\delta_n/R_0} [ K_n \cos(n\theta) + H_n \sin(n\theta) ]. \end{align} $$

The function $R_0(t)$ is the expanding base radius defined in (3.3), and $K_n$ , $H_n$ are constant Fourier coefficients that would be determined from initial conditions.

3.1. Stability condition

The result (3.18) gives rise to an important stability condition for the linearized solution. It follows from (3.3) that $R_0(t)> 1$ and, consequently, the interface perturbation $R_1(\theta ,t)$ in (3.18) grows with time at the nth Fourier mode if $\lambda _n> 0$ . This shows that the interface becomes immediately unstable if

(3.19) $$ \begin{align} \beta < \frac{n - 1/2}{n + 1/2}. \end{align} $$

For large mode number n, we find that equation (3.19) reduces to suggesting unstable behaviour when $\beta < 1$ , or in dimensional variables, $\mu _1 < \mu _2$ . From the physical viewpoint, instability thus occurs when the low-viscosity fluid forces the more viscous fluid to recede. Conversely, in the case where $\mu _2 < \mu _1$ (that is, $\beta> 1$ ) and the low-viscosity fluid is displaced by the more viscous one, the interface is always stable.

To see the effect of mode number and viscosity ratio upon stability, we write the amplitude of the nth mode in (3.18) in the form

$$ \begin{align*} R_0^{\lambda_n} {e}^{\delta_n/R_0} \equiv \exp \{ G_n(t) \} , \end{align*} $$

in which we have defined a growth function

(3.20) $$ \begin{align} G_n(t) & = \lambda_n \ln ( R_0) + \frac{\delta_n}{R_0} \nonumber \\ & = \frac{2 \sigma n(n^2-1)}{1 + \beta} \frac{1}{R_0(t)} + \bigg[ \frac{n(1-\beta)}{1 + \beta} - \frac{1}{2} \bigg] \ln ( R_0(t)). \end{align} $$

Here, we have used the definitions in (3.17). As time becomes large, the expression (3.20) behaves as

(3.21) $$ \begin{align} \frac{G_n(t)}{\ln ( R_0(t))} \sim \bigg[ \frac{n(1-\beta)}{1 + \beta} - \frac{1}{2} \bigg], \end{align} $$

and we have plotted contours of this quantity in Figure 2 for varying $\beta $ and n. The numerical values of the function (3.21) are shown on the diagram for the first ten contours. When this term is positive, the linearized solution (3.18) grows without bound and so is unstable, and this is indicated on the figure. Conversely, negative values result in a decaying, stable mode. The critical value $G_n = 0$ divides the parameter plane in Figure 2 into stable and unstable regions, as indicated on the diagram, and gives rise to the stability criterion (3.19). Figure 2 is slightly novel in the context of the Saffman–Taylor instability, although it is qualitatively similar to a linearized stability criterion derived by Grodzki and Szymczak [12] for injection of a corrosive acid into a porous substrate.

Figure 2 Stability diagram for modes $n=1, \ldots ,20$ over viscosity ratios $0 \leq \beta \leq 1$ . Coloured curves represent contours of the growth function given in (3.14). (Colour available online.)

4. Numerical solution for sharp interface

To understand the effects of nonlinearity, in the nonlinear model that consists of two distinct fluids separated by a sharp interface, a numerical semi-analytical approach is used to solve the full radial nonlinear problem. Exact solutions of Laplace’s equations (2.7) in each fluid domain are used and their time-dependent Fourier coefficients are found numerically as the solutions to coupled systems of differential equations. This builds on the “basic” scheme of Forbes et al. [10], extended to planar cylindrical geometry by Forbes [8].

In the two different fluid regions, the solutions to Laplace’s equation follow the overall mathematical structure (3.10) used for the linearized solution in Section 3. However, terms such as $r^n$ and $r^{-n}$ evaluated on an increasing radius can become arbitrarily large or small when the mode number n becomes large. This results in a numerical method that is too ill-conditioned to be of much practical use. To avoid this problem, the radius must be scaled with respect to some function of comparable size to the radius of the expanding interface and so we choose to take the scale function to be (3.3), which for clarity we rewrite here as

(4.1) $$ \begin{align} R_S (t) = \sqrt{1 + t}. \end{align} $$

Accordingly, the two pressures in this nonlinear problem are now taken to be

(4.2) $$ \begin{align} \begin{aligned} p_1 (r, \theta ,t) & = - \frac{\beta}{2} \ln(r) + A_0(t) + \sum_{n=1}^N A_n(t) \bigg( \frac{r}{R_S(t)} \bigg)^n \cos(n\theta) \quad \textrm{in } 0 < r < R(\theta ,t) , \\ p_2 (r, \theta ,t) & = - \frac{1}{2} \ln(r) + C_0(t) + \sum_{n=1}^N C_n(t) \bigg( \frac{r}{R_S(t)} \bigg)^{- n} \cos(n\theta) \quad \textrm{in } r> R(\theta ,t) \end{aligned} \end{align} $$

with the scaling function given by (4.1). Other forms for $R_S$ could doubtlessly be chosen, although these have not been investigated here.

The fluid boundary, located at $r = R(\theta ,t)$ , likewise has a form based on the linearized function in (3.18). We assume that

(4.3) $$ \begin{align} R(\theta,t) = R_0(t) + \sum_{n=1}^N R_n(t) \cos(n\theta). \end{align} $$

We observe that, geometrically, we are free to rotate the coordinates about the origin, so that the interface and the two nonlinear pressures can be chosen to be symmetric about the x-axis ( $\theta =0$ in polar coordinates). This is why only the (even) cosine functions are present in (4.2) and (4.3). Although the two pressures in (4.2) are in the forms of solutions to Laplace’s equation, which is a linear partial differential equation, the problem is nevertheless highly nonlinear, since the shape of the interface in (4.3) is unknown; therefore, the two fluid solution domains in (4.2) are not known in advance. The aim here is to solve for the Fourier coefficients $A_n(t)$ , $C_n(t)$ and $R_n(t)$ in these representations, and thus to determine the three unknown functions $p_1$ , $p_2$ and R in (4.2) and (4.3).

In the inner region 1, the fluid satisfies Darcy’s law (2.3), and so it has velocity $\mathbf {q}_1 = u_1 \mathbf {e_r} + v_1 \mathbf {e_{\theta }}$ with radial component $u_1$ and azimuthal component $v_1$ given by

(4.4) $$ \begin{align} \begin{aligned} &&u_1(r,\theta,t) &= -\frac{1}{\beta} \frac{\partial p_1}{\partial r} = \frac{1}{2r} - \frac{1}{\beta R_S(t)} \sum_{n=1}^N n A_n(t) \bigg( \frac{r}{R_S(t)} \bigg)^{n-1} \cos(n\theta) , \\ &&v_1(r,\theta,t) &= -\frac{1}{\beta r} \frac{\partial p_1}{\partial \theta} = \frac{1}{\beta R_S(t)} \sum_{n=1}^N n A_n(t) \bigg( \frac{r}{R_S(t)} \bigg)^{n-1} \sin(n\theta). \end{aligned} \end{align} $$

Similarly, the velocity $\mathbf {q}_2$ in outer fluid 2 is obtained by differentiation of the pressure $p_2$ according to Darcy’s law (2.4), resulting in components

(4.5) $$ \begin{align} \begin{aligned} &&u_2(r,\theta,t) &= - \frac{\partial p_2}{\partial r} = \frac{1}{2r} + \frac{1}{R_S(t)} \sum_{n=1}^N n C_n(t) \bigg( \frac{r}{R_S(t)} \bigg)^{-n-1} \cos(n\theta) , \\ &&v_2(r,\theta,t) &= -\frac{1}{r} \frac{\partial p_2}{\partial \theta} = \frac{1}{R_S(t)} \sum_{n=1}^N n C_n(t) \bigg( \frac{r}{R_S(t)} \bigg)^{-n-1} \sin(n\theta). \end{aligned} \end{align} $$

These forms are used to satisfy the nonlinear boundary conditions on the interface $r = R(\theta ,t)$ .

It follows from the kinematic boundary conditions that the zeroth-mode Fourier coefficient $R_0(t)$ in the representation (4.3) for the interface shape can be written in terms of the remaining coefficients in that expression. The expressions (4.4) for the two velocity components in fluid 1 are evaluated on the (unknown) interface $r = R(\theta ,t)$ and substituted into the first kinematic condition in (2.8). The expression is multiplied by R and solved for $R \partial R / \partial t$ to give

$$ \begin{align*} R \frac{\partial R}{\partial t} = \frac{1}{2} - \frac{1}{\beta} \sum_{n=1}^N \frac{n A_n(t)}{[ R_S(t)]^n} \bigg[ R^{n} \cos(n\theta) + R^{n-1} \frac{\partial R}{\partial \theta} \sin(n\theta) \bigg]. \end{align*} $$

This equation can be written in the simplified form

$$ \begin{align*} \frac{\partial}{\partial t} \bigg( \frac{1}{2} R^2 \bigg) = \frac{1}{2} - \frac{1}{\beta} \sum_{n=1}^N A_n(t) \frac{\partial}{\partial \theta} \bigg( \bigg[ \frac{R}{R_S(t)} \bigg]^n \sin(n\theta) \bigg). \end{align*} $$

The zeroth Fourier mode can be extracted from this at once, simply by integration. We obtain

$$ \begin{align*} \frac{{d}}{{d} t} \bigg( \int_{-\pi}^\pi R^2 \,d\theta \bigg) = 2 \pi. \end{align*} $$

Integration with respect to time now gives

(4.6) $$ \begin{align} \int_{-\pi}^\pi R^2(\theta,t) \, d\theta = 2 \pi t + \int_{-\pi}^\pi R^2(\theta,0)\, d\theta. \end{align} $$

The representation (4.3) is now substituted into this expression and the integrals are evaluated exactly by means of the orthogonality of the trigonometric functions. As an example, the term on the left-hand side of (4.6) becomes

$$ \begin{align*} \int_{-\pi}^\pi R^2(\theta,t) \,d\theta = 2 \pi \bigg[ R_0^2(t) + \frac{1}{2} \sum_{\ell =1}^N R_{\ell}^2(t) \bigg]. \end{align*} $$

As a result, the identity (4.6) yields

(4.7) $$ \begin{align} R_0(t) = \sqrt{ R_0^2(0) + t + \frac{1}{2} \sum_{\ell =1}^N [ R_{\ell}^2(0) - R_{\ell}^2(t)] }. \end{align} $$

This equation is the full nonlinear equivalent of the simple result (3.3).

It is now necessary to calculate the remaining Fourier coefficients $R_1 , \ldots , R_N$ , and ordinary differential equations (ODEs) are derived for these quantities from the first kinematic condition (2.8) with $j=1$ . Substituting the representations (4.4) for the two velocity components in fluid 1 and evaluating these at the interface $r = R(\theta ,t)$ reduces the first kinematic boundary condition to

(4.8) $$ \begin{align} & \frac{1}{2 R} - \frac{1}{\beta R_S(t)} \sum_{n=1}^N n A_n(t) \bigg[ \frac{R}{R_S(t)} \bigg]^{n-1} \cos(n\theta) \nonumber \\ &\quad = \frac{{d} R_0}{{d} t} + \sum_{n=1}^N \frac{{d} R_n}{{d} t} \cos(n\theta) + \frac{1}{\beta R_S(t)} \sum_{n=1}^N n A_n(t) \frac{R^{n-2}}{[ R_S(t)]^{n-1}} \frac{\partial R}{\partial \theta} \sin(n\theta). \end{align} $$

A differential equation for $R_0 (t)$ can be derived from this relation by integrating over $-\pi < \theta < \pi $ . This adds no new information, however, since the formula (4.7) already gives $R_0$ in terms of the higher coefficients. Consequently, this is not pursued further here. Differential equations for the higher coefficients $R_{\ell }(t)$ , $\ell = 1, 2, \ldots , N$ , however, are obtained by multiplying the first kinematic boundary condition (4.8) by the basis functions $\cos (\ell \theta )$ and integrating over $\theta $ . Integration by parts is used on the last term in (4.8) and, after some algebra, the desired ordinary differential equations are obtained in the form

(4.9) $$ \begin{align} \frac{{d} R_{\ell}}{{d} t} & = \frac{1}{2\pi} J_{\ell}^{(1)}(t) + \frac{A_1(t)}{\pi \beta R_S(t)} [ K^{(1)}_{\ell,1}(t) - \ell S^{(1)}_{\ell,1}(t) -\pi \delta_{\ell ,1} ] \nonumber \\ & \quad + \frac{1}{\pi \beta R_S(t)} \sum_{n=2}^N A_n(t) [K^{(1)}_{\ell ,n}(t) - \ell S^{(1)}_{\ell ,n}(t)]. \end{align} $$

In this expression, we have employed the usual Kronecker delta symbol

$$ \begin{align*} \delta_{\ell,1} = \begin{cases} 0 \quad\text{if } \ell \neq 1, \\ 1 \quad\text{if } \ell = 1. \end{cases} \end{align*} $$

In addition, we have defined intermediate functions

(4.10) $$ \begin{align} \begin{aligned} K^{(1)}_{\ell ,1}(t) & = \int_{-\pi}^\pi \ln ( R(\theta,t) ) \cos(\theta) \cos(\ell\theta) \, d\theta , \\ K^{(1)}_{\ell ,n}(t) & =\bigg(\frac{n}{n-1} \bigg) \int_{-\pi}^\pi \bigg[ \frac{R(\theta,t)}{R_S(t)} \bigg]^{n-1} \cos(n\theta)\cos(\ell\theta) \, d\theta , \quad n = 2, \ldots , N , \end{aligned} \end{align} $$

for $\ell = 1, \ldots , N$ . We observe that the structure of the first kinematic condition (4.8) requires that the case $n=1$ must be treated separately. Similarly, the further intermediate functions

(4.11) $$ \begin{align} \begin{aligned} S^{(1)}_{\ell ,1}(t) & = \int_{-\pi}^\pi \ln ( R(\theta,t) ) \sin(\theta) \sin(\ell\theta) \, d\theta , \\ S^{(1)}_{\ell ,n}(t) & = \bigg(\frac{n}{n-1} \bigg) \int_{-\pi}^\pi \bigg[ \frac{R(\theta,t)}{R_S(t)} \bigg]^{n-1} \sin(n\theta) \sin(\ell\theta) \, d\theta ,\quad n = 2, \ldots , N \end{aligned} \end{align} $$

have been defined here, along with the functions

(4.12) $$ \begin{align} J_{\ell}^{(1)}(t) = \int_{-\pi}^\pi \frac{1}{R(\theta,t)} \cos(\ell \theta) \, d\theta \end{align} $$

for $\ell = 1, 2, \ldots , N$ .

The second kinematic condition in (2.8), with $j=2$ , may be treated in the same way as for the first condition above. The representations (4.5) for the two velocity components in outer fluid 2 are incorporated into the second kinematic condition and evaluated at the interface $r = R(\theta ,t)$ . The resulting equation is again multiplied by basis functions $\cos (\ell \theta )$ , $\ell = 1, 2, \ldots , N$ , and integrated over $- \pi < \theta < \pi $ . The algebra is simpler than for the first kinematic condition, since the case $n=1$ in the sums does not need to be treated separately. We obtain

(4.13) $$ \begin{align} \frac{{d}R_{\ell}}{{d}t} = \frac{1}{2\pi} J_{\ell}^{(1)}(t) + \frac{1}{\pi R_S(t)} \sum_{n=1}^N C_n(t)[K^{(2)}_{\ell ,n}(t) + \ell S^{(2)}_{\ell ,n}(t)] , \quad \ell = 1, \ldots ,N. \end{align} $$

Here, we have defined extra intermediate functions

(4.14) $$ \begin{align} \begin{aligned} K^{(2)}_{\ell ,n}(t) & = \bigg(\frac{n}{n+1}\bigg) \int_{-\pi}^{\pi} \bigg[ \frac{R(\theta ,t)}{R_S(t)} \bigg]^{-n-1} \cos(n\theta) \cos(\ell\theta) \, d\theta \\ S^{(2)}_{\ell ,n}(t) & = \bigg(\frac{n}{n+1}\bigg) \int_{-\pi}^{\pi} \bigg[ \frac{R(\theta ,t)}{R_S(t)} \bigg]^{-n-1} \sin(n\theta) \sin(\ell\theta) \, d\theta \end{aligned} \end{align} $$

for $\ell = 1, 2, \ldots ,N$ .

The second kinematic condition in the form (4.13) gives a system of differential equations for $R_{\ell }(t)$ . However, this information is already provided by (4.9) and so, to extract as much independent information as possible from these two kinematic conditions, we effectively replace the first condition (4.9) with the difference of that equation and the second kinematic condition (4.13). This now generates the algebraic system

(4.15) $$ \begin{align} 0 &= \sum_{n=1}^N [K_{\ell ,n}^{(2)}(t) + \ell S_{\ell ,n}^{(2)}(t)] C_n(t) - \frac{1}{\beta} [ K^{(1)}_{\ell ,1}(t) -\ell S^{(1)}_{\ell ,1}(t) -\pi \delta_{\ell ,1} ] A_1(t)\nonumber \\ & \quad - \frac{1}{\beta} \sum_{n=2}^N [K^{(1)}_{\ell ,n}(t) - \ell S^{(1)}_{\ell ,n}(t)] A_n(t). \end{align} $$

This is now a matrix equation that relates the coefficients $C_1 , \ldots , C_N$ to the coefficients $A_1 , \ldots , A_N$ at each value of time t.

Finally, the dynamic condition (2.9) is subjected also to Fourier analysis, by multiplying by basis functions $\cos (\ell \theta )$ , $\ell = 1, 2, \ldots , N$ and integrating over $\theta $ , as previously. We ignore the zeroth mode $\ell = 0$ which effectively provides extra information about the difference between the two coefficients $A_0$ and $C_0$ in the representations (4.2) for pore pressures; this information is not needed here because pressure is arbitrary within an additive constant. The $\ell $ th Fourier mode for the dynamic condition (2.9) can be written as

(4.16) $$ \begin{align} \sum_{n=1}^N \mathcal{M}_{\ell ,n}(t) C_n(t) - \sum_{n=1}^N \mathcal{N}_{\ell ,n}(t) A_n(t) = \frac{1}{2}(1-\beta) \mathcal{L}_\ell(t) -\sigma \mathcal{F}_\ell (t) \end{align} $$

for $\ell = 1, 2, \ldots , N$ . Again, in view of the intricacy of these conditions, it has proven very convenient to introduce further intermediate matrix quantities

(4.17) $$ \begin{align} \begin{aligned} \mathcal{M}_{\ell ,n}(t) & = \int_{-\pi}^\pi \bigg[ \frac{R(\theta ,t)}{R_S(t)} \bigg]^{-n} \cos(n\theta) \cos(\ell\theta) \, d\theta \\[.2cm] \mathcal{N}_{\ell ,n}(t) & = \int_{-\pi}^\pi \bigg[ \frac{R(\theta ,t)}{R_S(t)} \bigg]^{n} \cos(n\theta) \cos(\ell\theta) \, d\theta \end{aligned} \end{align} $$

for $ \ell = 1, 2, \ldots , N$ and $n= 1, 2, \ldots , N$ , along with the vector terms

(4.18) $$ \begin{align} \begin{aligned} \mathcal{L}_{\ell}(t) & = \int_{-\pi}^\pi \ln [ R(\theta ,t)] \cos(\ell\theta) \, d\theta \\ \mathcal{F}_{\ell}(t) & = \int_{-\pi}^\pi \kappa (\theta ,t) \cos(\ell \theta) \, d\theta , \quad \ell = 1, 2, \ldots , N. \end{aligned} \end{align} $$

The constant $\sigma $ is the surface-tension parameter in (2.10) and $\kappa $ is the curvature in (2.11).

4.1. Computational algorithm

The spectral representations (4.2) already solve the governing (Laplace) equations (2.7) exactly for the pore pressure and it is only required to calculate the coefficients in these Fourier series along with the coefficients in (4.3). The three sets of Fourier coefficients in these representations are found from the three sets of equations that derive from the boundary conditions on the interface.

In our solution algorithm, we regard the Fourier coefficients in (4.3) for the interface shape as the fundamental unknown quantities to be found. We therefore create the $(N \times 1)$ vector of unknowns

(4.19) $$ \begin{align} \mathbf{X}(t) = [ R_1(t) , R_2(t) , \ldots , R_N(t)]^T , \end{align} $$

and seek to write the kinematic and dynamic conditions at the interface as an N-dimensional dynamical system of the form

(4.20) $$ \begin{align} \frac{d\mathbf{X}}{d t}= \mathrm{F} ( t , \mathbf{X}). \end{align} $$

This system (4.20) is marched forward in time using a matlab routine for solving ODEs. The vector (4.19) is thus obtained at later times t, starting from some initial configuration $\mathbf {X}(0)$ , and the solution variables are then reconstructed from their spectral representations. For this purpose, a dedicated matlab routine is written that takes the vector $\mathbf {X}$ in (4.19) at time t as its input and returns the vector of derivatives $\mathbf {F}$ in (4.20) as output.

This dedicated routine therefore takes the vector $\mathbf {X}$ of coefficients $R_1 , \ldots , R_N$ at time t and, from these, it immediately calculates the average interface radius $R_0(t)$ according to (4.7). Next, the interface shape $R(\theta; t)$ and the first two spatial derivatives $\partial R / \partial \theta $ and $\partial ^2 R / \partial \theta ^2$ are re-constructed from (4.3), and the curvature is obtained from (2.11). Following this, the matrices and vectors of intermediate quantities (4.10)–(4.12) and (4.14), (4.17), (4.18) are calculated and stored. All these functions involve integrations over the spatial domain $-\pi < \theta < \pi $ and we evaluate these to extremely high accuracy using the Gaussian quadrature routine lgwt.m written by von Winckel [29].

It remains for this interior routine to calculate the Fourier coefficients $A_j(t)$ , $C_j(t)$ , $j=1, \ldots ,N$ , in the representations (4.2) for the two pressures. These are obtained from the two conditions (4.15) and (4.13), which constitute a matrix-vector equation of the form

(4.21) $$ \begin{align} \mathbf{M}_L \, \mathbf{v} = \mathbf{R}_H , \end{align} $$

in which $\mathbf {v}$ denotes the augmented vector

(4.22) $$ \begin{align} \mathbf{v} = [ A_1(t) , \ldots , \ A_N(t); C_1(t), \ldots , C_N(t)]^T \end{align} $$

of size $(2N \times 1)$ . The coefficient matrix $\mathbf {M}_L$ in the matrix system (4.21) is of size $(2N \times 2N)$ and has the partitioned form

$$ \begin{align*} \mathbf{M}_L \equiv \left[ \begin{array}{@{}ll@{}} \mathbf{m^{KA}} & \mathbf{m^{KC}} \\ \mathbf{m^{DA}} & \mathbf{m^{DC}} \end{array}\right] , \end{align*} $$

in which the two submatrices coming from the kinematic conditions (4.15) have components

$$ \begin{align*} [ \mathbf{m^{KA}} ]_{\ell ,1} & = - \frac{1}{\beta} [K^{(1)}_{\ell ,1}(t) -\ell S^{(1)}_{\ell ,1}(t) - \pi \delta_{\ell ,1}]\quad \textrm{if } n=1 ,\nonumber \\ [ \mathbf{m^{KA}} ]_{\ell ,n} & = - \frac{1}{\beta} [K^{(1)}_{\ell ,n}(t) - \ell S^{(1)}_{\ell ,n}(t)] \quad \textrm{if } n=2, \ldots, N, \nonumber \\ [ \mathbf{m^{KC}} ]_{\ell ,n} & =K_{\ell ,n}^{(2)}(t) + \ell S_{\ell ,n}^{(2)}(t), \nonumber \end{align*} $$

and the two submatrices that arise from the dynamic condition (4.16) have components

$$ \begin{align*} [ \mathbf{m^{DA}}]_{\ell ,n} & = - \mathcal{N}_{\ell ,n}(t) , \nonumber \\ [ \mathbf{m^{DC}}]_{\ell ,n} & = \mathcal{M}_{\ell ,n}(t). \nonumber \end{align*} $$

The quantities $\mathcal {M}_{\ell ,n}$ and $\mathcal {N}_{\ell ,n}$ are defined in (4.17). The right-hand side vector $\mathbf {R}_H$ in the matrix equation (4.21) is of size $(2N \times 1)$ . Its first N elements are all zero so that it has the structure

$$ \begin{align*} [ \mathbf{R}_H ]_{\ell} = [ 0, \ldots , 0; \tfrac{1}{2}(1-\beta)\mathcal{L}_{\ell} - \sigma \mathcal{F}_{\ell} ]^T \end{align*} $$

with $\ell = 1, \ldots , N$ . The matrix equation (4.21) is solved for the vector $\mathbf {v}$ in (4.22) to give the coefficients $A_n$ and $C_n$ , and the two pressures are re-constructed from (4.2). Finally, the vector

$$ \begin{align*} \mathbf{F} = [ R^{\prime}_1(t) , R^{\prime}_2(t) , \ldots , R^{\prime}_N(t)]^T \end{align*} $$

in (4.20) is calculated from (4.13) and the system is integrated forward in time.

4.2. Smoothing algorithm

When a function is represented by a Fourier series, as is the case for the interface function $R(\theta ,t)$ in (4.3), for example, the convergence of the series becomes a concern if the original function R changes rapidly with $\theta $ . In particular, the Fourier series representation of a discontinuous function R does not converge to the original function, but instead to a function which changes rapidly at the location of the discontinuity and possesses spurious oscillations either side of it. This is known as the Gibbs phenomenon (see Kreyszig [20]) and it can be responsible for producing small-amplitude ripples near regions of interest in our numerical solutions. Where appropriate, we therefore smooth the discontinuous original function using Lanczos smoothing, which can be carried out simply by multiplying each Fourier coefficient $R_{n} (t)$ calculated in (4.3) by a damping factor

(4.23) $$ \begin{align} \frac{\sin (\lambda_T n)}{\lambda_T n} , \quad n = 1, 2, \ldots , N. \end{align} $$

Here, the small positive constants $\lambda _T$ are suitably chosen Lanczos parameters. Further details are given by Duchon [6]. It can be shown that this is equivalent to replacing the discontinuous original function $R(\theta ,t)$ in (4.3), at a given location $\theta ^{*}$ , with a continuous function that is the average of (4.3) over the moving interval $\theta ^{*} - \lambda _T < \theta < \theta ^{*} + \lambda _T$ . We typically choose the Lanczos parameter to have a value of approximately $\lambda _T \approx 0.05$ .

5. Numerical solution for diffuse interface

The nonlinear model described in Section 4 assumes that there exists an infinitesimally thin interface between the fluids. As a consequence of this approximation, that nonlinear model can develop Moore curvature singularities at its interface, resulting in the failure of the entire model within finite time [22]. For real viscous fluids, however, there is likely to be a narrow interfacial layer across which the viscosity changes rapidly but smoothly, rather than a precise interface of zero thickness at which the viscosity jumps discontinuously between the two fluids. For this reason, we consider a further nonlinear model in this section, in which the system is treated as a single fluid with viscosity that changes continuously in space. A narrow interfacial region will be present, across which the viscosity changes rapidly from one region to another, representing an interfacial zone where mixing of the two fluids could occur.

Possibly the simplest model of the type wanted in this section is one in which the viscosity depends upon temperature. In dimensionless variables, the temperature $T(r,\theta ,t)$ satisfies the heat equation

(5.1) $$ \begin{align} \frac{\partial T}{\partial t} + \mathbf{q}\cdot \nabla T = \sigma_H \nabla^2 T. \end{align} $$

This equation introduces a new dimensionless temperature-diffusion parameter ${\sigma _H = \kappa _H (2 \pi h)/Q}$ , in which the dimensional diffusion constant $\kappa _H$ has units of length squared per time. To some extent, this diffusion constant $\sigma _H$ plays a similar role to the surface-tension parameter $\sigma $ defined in (2.10), since it also is responsible for diffusive effects in the vicinity of the interface. We suppose that fluid is still injected through a line source at the origin, as previously, and that it has viscosity $\mu _1$ and temperature $T_1$ in dimensional variables; far away, the fluid has viscosity $\mu _2$ and temperature $T_2$ . We now return to dimensionless variables by referencing all temperatures with respect to the outer temperature $T_2$ so that the nondimensional temperature T in the heat equation (5.1) has the value $1$ far away and the dimensionless value

(5.2) $$ \begin{align} \tau = T_1/T_2 \end{align} $$

at the origin $r=0$ . Now the dimensionless form of Darcy’s law in the porous medium, replacing (2.3) and (2.4), becomes

(5.3) $$ \begin{align} \mathbf{q} = -\frac{1}{\mu(T)} \nabla p. \end{align} $$

In the simplest model, the viscosity is taken to vary linearly with temperature and we assume the simple “state” equation

(5.4) $$ \begin{align} \mu(T) = \beta - \frac{\beta -1}{1 - \tau} (T - \tau). \end{align} $$

Here, $\beta = \mu _1/\mu _2$ is the same viscosity ratio (2.2) used previously and $\tau $ is the ratio (5.2) of inner to outer fluid temperatures. The single fluid is taken to be incompressible, so that

(5.5) $$ \begin{align} \nabla \cdot \mathbf{q} = 0 , \end{align} $$

replacing (2.5) and, since a line source is present at the origin, the condition (2.6) still holds there.

As for the previous nonlinear model in Section 4, we again make use of a spectral representation to solve this nonlinear problem. Now, however, we are forced to restrict the computational domain by placing some outer cylinder $r = R_{\infty }$ appropriately far from the origin, and expressing the temperature as a double Fourier series in r and $\theta $ , with time-dependent doubly subscripted Fourier coefficients that are to be determined. The boundary conditions for temperature then become $T = \tau $ at $r=0$ and $T = 1$ at $r = R_\infty $ . After some experimentation, it was determined that the temperature could be well represented by the spectral form

(5.6) $$ \begin{align} T(r,\theta ,t) & = \frac{1}{2}(\tau+1) + \frac{1}{2}(\tau -1) \cos \bigg( \frac{\pi r}{R_{\infty}} \bigg) + \sum_{n=1}^N A_{0,n}(t) \sin \bigg( \frac{n \pi r}{R_{\infty}} \bigg)\nonumber \\ & \quad + \sum_{m=1}^M \sum_{n=1}^N A_{m,n}(t) J_m(\gamma_{m,n}r) \cos(m\theta). \end{align} $$

In this expression, we have defined constants

(5.7) $$ \begin{align} \gamma_{m,n} = j_{m,n} / R_{\infty} , \end{align} $$

in which the constant $j_{m,n}$ denotes the nth positive zero of the first-kind Bessel function $J_m (z)$ . With the choice of constants in (5.7), the assumed form (5.6) for the temperature function clearly satisfies the required boundary conditions at $r = 0$ and $r = R_{\infty }$ .

A representation is also needed for the pressure in the single fluid. This function must behave like the base pressure in (3.4) as $r \rightarrow 0$ , to take account of the line source at the origin. Accordingly, we have chosen

(5.8) $$ \begin{align} p(r,\theta;t) = - \frac{\beta}{2} \ln(r) + \sum_{n=1}^N B_{0,n}(t) J_0(\gamma_{0,n} r) + \sum_{m=1}^M \sum_{n=1}^N B_{m,n}(t) J_m(\gamma_{m,n} r) \cos(m\theta). \end{align} $$

This again makes use of the constants defined in (5.7) and the coefficient functions $B_{m,n}$ are to be determined. A system of ODEs for the coefficients $A_{m,n}$ is obtained by substituting the form (5.6) into the heat equation (5.1) to give

(5.9) $$ \begin{align} &\sum_{n=1}^N A_{0,n}^{\prime}(t) \sin \bigg(\frac{n \pi r}{R_\infty} \bigg) + \sum_{m=1}^M \sum_{n=1}^N A_{m,n}^{\prime}(t) J_m(\gamma_{m,n} r) \cos(m\theta)\nonumber \\ &\quad = - \frac{\sigma_H}{2}(\tau -1) \frac{\pi}{R_\infty} \frac{1}{r} \bigg[ \sin \bigg( \frac{\pi r}{R_\infty} \bigg)+ \frac{\pi r }{R_\infty} \cos \bigg( \frac{\pi r}{R_\infty} \bigg) \bigg]\nonumber \\ & \qquad + \sigma_H \sum_{n=1}^N \frac{\pi }{R_\infty} A_{0,n}(t) \frac{n}{r} \bigg[ \cos \bigg( \frac{n \pi r}{R_\infty} \bigg) - \frac{n \pi r}{R_\infty} \sin \bigg(\frac{n \pi r}{R_\infty} \bigg) \bigg]\nonumber \\ &\qquad - \sigma_H \sum_{m=1}^M \sum_{n=1}^N A_{m,n}(t) \gamma_{m,n}^2 J_m(\gamma_{m,n} r) \cos(m\theta)- \bigg( u \frac{\partial T}{\partial r}+ \frac{v}{r} \frac{\partial T}{\partial \theta} \bigg). \end{align} $$

To obtain this expression, we have used the Bessel differential equation given by Abramowitz and Stegun [1, p. 358].

Differential equations for the zeroth (azimuthal) mode coefficients are now obtained, by integrating (5.9) over $- \pi < \theta < \pi $ . The resulting equation is then multiplied by basis functions $\sin ( (\ell \pi r)/R_\infty )$ and integrated over $0 < r < R_\infty $ . We obtain, after rearrangement,

(5.10) $$ \begin{align} \frac{{d} A_{0,\ell}(t)}{{d} t} &= - \sigma_H (\tau -1) \bigg( \frac{\pi}{R_{\infty}^2} \bigg) \mathbf{V}_{\ell}^s - \sigma_H (\tau-1) \bigg( \frac{\pi^2}{R_\infty^3}\bigg) \mathcal{S}_{\ell} + \frac{2 \sigma_H}{R_\infty} \sum_{n=1}^N \mathbf{M}_{\ell ,n}^c A_{0,n}(t) \nonumber \\[4pt] & \quad - \sigma_H \bigg( \frac{\ell \pi}{R_{\infty}} \bigg)^2 A_{0,\ell}(t) - \frac{1}{\pi R_{\infty}} \int_0^{R_{\infty}} \int_{-\pi}^\pi \bigg(r u \frac{\partial T}{\partial r} + v \frac{\partial T}{\partial \theta} \bigg) \frac{1}{r} \sin \bigg( \frac{\ell \pi r}{R_{\infty}} \bigg) \, d\theta \, dr. \end{align} $$

Here, we have defined the intermediate vector

(5.11) $$ \begin{align} \mathcal{S}_{\ell} & \equiv \int_0^{R_\infty} \cos \bigg( \frac{\pi r}{R_\infty} \bigg) \sin \bigg(\frac{\ell \pi r}{R_\infty} \bigg) \, dr \nonumber \\[4pt] &= \begin{cases} 0 & \text{if } \ell = 1,\\ [ R_{\infty} \ell ] [1 + \cos(\ell \pi)] / [ \pi (\ell^2 -1 ) ] & \text{if }\ell \neq 1 , \end{cases} \end{align} $$

and a further $(N \times 1)$ vector $\mathbf {V}_{\ell }^s$ and $(N \times N$ ) matrix $\mathbf {M}_{\ell ,n}^c$ defined as

$$ \begin{align*} \mathbf{V}_{\ell}^s & = \int_0^{R_\infty} \frac{1}{r} \sin \bigg( \frac{\pi r}{R_\infty} \bigg) \sin \bigg( \frac{\ell\pi r}{R_\infty} \bigg)\,dr\quad \textrm{for } \ell = 1,2,\ldots ,N , \\[4pt]\mathbf{M}_{\ell ,n}^c & = \int_0^{R_\infty}\bigg( \frac{n \pi}{R_\infty} \bigg) \frac{1}{r}\sin \bigg( \frac{\ell \pi r}{R_\infty} \bigg)\cos \bigg( \frac{n \pi r}{R_\infty} \bigg)\, dr \quad \textrm{for } n, \ell = 1,2,\ldots ,N. \end{align*} $$

Differential equations for the kth azimuthal mode, for $k = 1,2,\ldots ,M$ , are similarly calculated by multiplying (5.9) by $\cos (k\theta )$ and integrating over $\theta $ . The equation thus obtained is then multiplied by $r J_k(\gamma _{k,\ell } r)$ and integrated over r. We make use of the orthogonality condition

(5.12) $$ \begin{align} \int_0^{R_\infty} r J_k(\gamma_{k,n} r) J_k(\gamma_{k,\ell} r) \, dr = \begin{cases} 0 & \textrm{if } n \neq \ell ,\\ (1/2) [ R_\infty J_{k+1}(\gamma_{k,\ell} R_\infty)]^2& \text{if } n = \ell \end{cases} \end{align} $$

for Bessel functions, which may be derived from [1, p. 485, (11.4.5)], and obtain

(5.13) $$ \begin{align} \frac{{d} A_{k,\ell}(t)}{{d} t} = - \sigma_H \gamma_{k,\ell}^2 A_{k,\ell}(t) - \frac{2}{\mathcal{C}_{k,\ell}} \int_0^{R_\infty} \int_{-\pi}^\pi \bigg( r u \frac{\partial T}{\partial r} + v \frac{\partial T}{\partial \theta} \bigg) J_k(\gamma_{k,\ell} r) \cos(k \theta) \, d\theta \, d r. \end{align} $$

In this expression, it is convenient to write

(5.14) $$ \begin{align} \mathcal{C}_{k,\ell} = \pi [ R_\infty J_{k+1}(\gamma_{k,\ell} R_\infty)]^2 , \end{align} $$

which comes from the orthogonality condition (5.12).

The system of ordinary differential equations (5.10), (5.13) can in principle be integrated forward in time, once the Fourier coefficients $B_{m,n}$ in the expression (5.8) for pressure are known. To obtain them, use must be made of Darcy’s law (5.3), which we here write in the form $\nabla p = - \mu (T) \mathbf {q}$ . The divergence of this expression is taken, leading to

$$ \begin{align*} \nabla^2 p = - \bigg[\frac{{d} \mu}{{d} T} \nabla T \cdot \mathbf{q} + \mu(T) \textrm{div}(\mathbf{q}) \bigg]. \end{align*} $$

The second term on the right-hand side vanishes, by the incompressibility condition (5.5), and the viscosity–temperature relation (5.4) gives

$$ \begin{align*} \nabla^2 p = - \bigg( \frac{1- \beta}{1 - \tau} \bigg) \nabla T \cdot \mathbf{q}. \end{align*} $$

Darcy’s law (5.3) is used again to eliminate $\mathbf {q}$ , resulting in the equation

$$ \begin{align*} \nabla^2 p = \bigg( \frac{1- \beta}{1 - \tau} \bigg) \frac{1}{\mu (T)} \bigg[ \frac{\partial p}{\partial r} \frac{\partial T}{\partial r} + \frac{1}{r^2} \frac{\partial p}{\partial \theta} \frac{\partial T}{\partial \theta} \bigg]. \end{align*} $$

The representation (5.8) for pressure is substituted into this equation and yields

(5.15) $$ \begin{align} \begin{aligned} & - \sum_{n=1}^N \gamma_{0,n}^2 B_{0,n}(t) J_0(\gamma_{0,n}r) - \sum_{n=1}^N \sum_{m=1}^M \gamma_{m,n}^2 B_{m,n}(t) J_m(\gamma_{m,n}r) \cos(m\theta) \\ &\quad = \bigg( \frac{1 -\beta}{1 -\tau} \bigg) \frac{1}{\mu(T)} \bigg[ \frac{\partial p}{\partial r} \frac{\partial T}{\partial r} + \frac{1}{r^2} \frac{\partial p}{\partial \theta}\frac{\partial T}{\partial \theta} \bigg]. \end{aligned} \end{align} $$

This new expression (5.15) replaces the simple Laplace equation for pressure (2.7) that was used in Section 4 for the nonlinear model with a sharp interface.

Equation (5.15) is again subjected to Fourier analysis, in a similar manner to the heat equation discussed above, making use of the orthogonality of the simple trigonometric functions and the orthogonality relation (5.12) for Bessel functions. The zeroth Fourier azimuthal modes are obtained first by integrating over $\theta $ then by multiplying by $r J_0 (\gamma _{0,\ell } r)$ , $\ell = 1, \ldots , N$ , and integrating over r. This results in

(5.16) $$ \begin{align} \!\!\!\!B_{0,\ell}(t) = - \bigg( \frac{\beta - 1}{\tau - 1} \bigg) \frac{1}{\gamma_{0,\ell}^2 \mathcal{C}_{0,\ell}} \int_0^{R_\infty} \!\int_{-\pi}^\pi \frac{1}{\mu(T)} \bigg[ r\frac{\partial p}{\partial r} \frac{\partial T}{\partial r} + \frac{1}{r} \frac{\partial p}{\partial \theta}\frac{\partial T}{\partial \theta} \bigg] J_0(\gamma_{0,\ell} r) \, d\theta \, dr.\quad \end{align} $$

The constants $\mathcal {C}_{0,\ell }$ in this formula are obtained from (5.14) with $k=0$ . In a similar way, the equations for the higher azimuthal Fourier modes $\cos (k\theta )$ with $k = 1, \ldots , M$ are

(5.17) $$ \begin{align} B_{k,\ell}(t) & = - \bigg( \frac{\beta -1}{\tau -1 } \bigg) \frac{2}{\gamma_{k,\ell}^2 \mathcal{C}_{k,\ell}} \int_0^{R_\infty} \int_{-\pi}^{\pi} \frac{1}{\mu(T)} \bigg[ r\frac{\partial p}{\partial r} \frac{\partial T}{\partial r} \nonumber \\ &\quad + \frac{1}{r} \frac{\partial p}{\partial \theta}\frac{\partial T}{\partial \theta} \bigg] J_k(\gamma_{k,\ell}r) \cos(k\theta) \, d\theta \, d r. \end{align} $$

5.1. Iterative solution

Before the differential equations (5.10), (5.13) can be integrated forward in time, the Fourier coefficients $B_{m,n}$ must be determined from (5.16) and (5.17). These equations essentially constitute a straightforward system of linear algebraic equations that give $B_{m,n}$ in terms of the coefficients $A_{m,n}$ . In principle, therefore, they could be solved for the $B_{m,n}$ coefficients using a standard package that implements Gaussian elimination. The practical difficulty, however, is that the vector of unknown coefficients $B_{m,n}$ , $m = 0, 1, 2, \ldots , M$ , $n = 1, 2, \ldots , N$ would contain $(M + 1) N$ elements; its coefficient matrix would be of size $(M+1)N \times (M+1)N$ and, for larger values of M and N, overwhelms both the memory capacity and speed of many current computers. A similar difficulty, when solving problems of this general type, was encountered by Walters et al. [30] and Forbes et al. [11]. To overcome this difficulty, those authors solve their matrix equation using a type of fixed-point iteration scheme and we follow that approach here too.

An initial guess is made for the coefficients $B_{0,\ell }^{(0)}(t)$ and $B_{k,\ell }^{(0)}(t)$ , $k = 1, 2, \ldots , M$ , $\ell = 1, 2, \ldots , N$ . Based on this guess, the pressure $p^{(0)}$ and its first partial derivatives are re-constructed from the spectral representation (5.8). The system of linear equations (5.16), (5.17) is then solved iteratively in the form

(5.18) $$ \begin{align} \begin{aligned} B_{0,\ell}^{(q+1)}(t) & = - \bigg( \frac{\beta - 1}{\tau - 1} \bigg) \frac{1}{\gamma_{0,\ell}^2 \mathcal{C}_{0,\ell}} \int_0^{R_\infty} \int_{-\pi}^\pi \frac{1}{\mu(T)} \bigg[ r\frac{\partial p^{(q)}}{\partial r} \frac{\partial T}{\partial r} \\ & \quad + \frac{1}{r} \frac{\partial p^{(q)}}{\partial \theta}\frac{\partial T}{\partial \theta} \bigg] J_0(\gamma_{0,\ell} r) \, d\theta \, dr , \\ B_{k,\ell}^{(q+1)}(t) & = - \bigg( \frac{\beta -1}{\tau -1 } \bigg) \frac{2}{\gamma_{k,\ell}^2 \mathcal{C}_{k,\ell}} \int_0^{R_\infty} \int_{-\pi}^\pi \frac{1}{\mu(T)} \bigg[ r \frac{\partial p^{(q)}}{\partial r} \frac{\partial T}{\partial r} \\ & \quad + \frac{1}{r} \frac{\partial p^{(q)}}{\partial \theta}\frac{\partial T}{\partial \theta} \bigg] J_k(\gamma_{k,\ell}r) \cos(k\theta) \, d\theta \, d r. \end{aligned} \end{align} $$

This iteration scheme $q = 1, 2, 3, \ldots $ is continued until convergence is achieved, in the form

$$ \begin{align*} \| B^{(q+1)} - B^{(q)} \| < \delta \end{align*} $$

for some acceptable convergence criterion defined by the small constant $\delta $ .

5.2. Initial condition for temperature

The initial temperature in the region is chosen so as to mimic a disturbance to the initial circle $r = 1$ at some azimuthal Fourier mode $m^{*}$ . This then enables comparison with the predictions of the linearized solution in Section 3 and the sharp-interface nonlinear solution in Section 4. To do this, we create the discontinuous initial temperature

(5.19) $$ \begin{align} T(r,\theta ,0) = \begin{cases} \tau &\textrm{if } r < R_{in} (\theta ), \\ 1 &\textrm{if } R_{in} (\theta ) < r < R_{\infty} , \end{cases} \end{align} $$

in which

$$ \begin{align*} R_{in} = 1 + \epsilon_T \cos ( m^{*} \theta ). \end{align*} $$

Here, the constant $\epsilon _T$ gives the amplitude of the disturbance in physical space.

The initial values of the Fourier coefficients $A_{m,n} (0)$ are now calculated using Fourier analysis of the representation (5.6), as described above. After use of the orthogonality criterion (5.12),

(5.20) $$ \begin{align} A_{0,\ell} (0) & = - \frac{\tau +1}{\ell\pi} [1 - \cos(\ell\pi)] - \frac{\tau - 1}{R_{\infty}} \mathcal{S}_{\ell} \nonumber \\ & \quad + \frac{1}{\pi R_\infty} \int_0^{R_\infty} \int_{-\pi}^\pi T(r,\theta ,0) \sin \bigg( \frac{\ell\pi r}{R_\infty} \bigg) \, d\theta \, d r \end{align} $$

for the zeroth azimuthal modes and

(5.21) $$ \begin{align} A_{k,\ell} (0) = \frac{2}{\mathcal{C}_{k,\ell}} \int_0^{R_\infty} \int_{-\pi}^\pi r T(r,\theta ,0) J_k (\gamma_{k,\ell} r) \cos(k\theta) \, d\theta \, d r \end{align} $$

for the higher azimuthal modes $k = 1, 2, \ldots , M$ . These formulae make use of the chosen initial temperature in (5.19). The constants $S_{\ell }$ in (5.20) are as defined previously in (5.11) and the quantities $\mathcal {C}_{k,\ell }$ in (5.21) are given in (5.14). As discussed in Section 4.2, however, the discontinuous initial profile (5.19) is not particularly desirable, since its Fourier representation is subject to spurious oscillations caused by the Gibbs phenomenon [20]. We avoid this difficulty here using Lanczos smoothing [6] in both spatial variables r and $\theta $ . The coefficients $A_{k,\ell } (0)$ calculated in (5.21) are multiplied by damping factors

(5.22) $$ \begin{align} \bigg[ \frac{\sin (\lambda_T k)}{\lambda_T k} \bigg] \bigg[ \frac{\sin (\lambda_R \ell )}{\lambda_R \ell} \bigg] , \end{align} $$

in which the small positive constants $\lambda _T$ and $\lambda _R$ are suitably chosen Lanczos parameters. This generalizes the single-variable case (4.23) and is equivalent to replacing the discontinuous initial condition (5.19), at a given point $(r^{*} , \theta ^{*} )$ , with a continuous function that is the moving average of (5.19) over the rectangular patch $r^{*} - \lambda _R < r < r^{*} + \lambda _R$ , $\theta ^{*} - \lambda _T < \theta < \theta ^{*} + \lambda _T$ in the $(r,\theta )$ plane. Thus, the two Lanczos parameters $\lambda _T$ and $\lambda _R$ determine the effective width of the initial interfacial zone. This Lanczos smoothing strategy is also applied approximately a further ten times over the lifetime of our numerical solution, so as to regularize the results.

5.3. Computational algorithm

The computational process used to simulate this nonlinear varying-viscosity problem is similar in philosophy to that described in Section 4.1 for the previous nonlinear model in which there was a sharp discontinuous interface between the two fluids. In this model, however, there are many more Fourier coefficients that must be found and so computer run-times can be lengthy. To reduce as much as possible the demands on computer time, we cache as many functions as possible in large arrays, so that they may be calculated only once, stored and not re-computed. As examples, trigonometric functions such as $\cos ( m\theta _q )$ evaluated for modes $m = 1, 2, \ldots ,M$ and at mesh points $\theta _q$ with $q = 1, 2, \ldots , Q$ are stored in $(M \times Q)$ matrices. Similarly, Bessel functions $J_m ( \gamma _{m,n} r_p )$ evaluated at grid points $p = 1, 2, \ldots , P$ are stored in $(M \times N \times P)$ arrays.

As with the previous nonlinear problem in Section 4.1, the overall algorithm here treats the differential equations (5.10), (5.13) as a dynamical system for the coefficients $A_{m,n}$ in which the additional coefficients $B_{m,n}$ must be determined at each new time step. This system again has an overall form equivalent to that in (4.20), except that, here, the unknown vector of coefficients is very large, of size $( (M+1)N \times 1 )$ in matrix notation, and has the structure

$$ \begin{align*} \mathbf{X}(t) = [ A_{0,1}, \ldots , A_{0,N} \; \ \ldots\ldots \; \ A_{m,1}, \ldots ,\ A_{m,N} \; \ \ldots\ldots \; \ A_{M,1}, \ldots ,\ A_{M,N}]^T. \end{align*} $$

The algorithm begins by calculating the initial coefficients for temperature, as described in Section 5.2, and storing these in a vector $\mathbf {X}(0)$ . Then, at each time step, a dedicated inner code takes the vector $\mathbf {X}$ of coefficients as its input and uses it to solve for the coefficients $B_{m,n}$ at that time step by the iteration process described in Section 5.1. The routine then calculates a vector $\mathbf {X}^{\prime } = \mathbf {F}$ that consists of all the derivatives of these Fourier coefficients, by evaluating the expressions (5.10), (5.13). This vector $\mathbf {F}$ is also of dimension $( (M+1)N \times 1 )$ , having $(M+1)N$ components; it is passed to a matlab routine that solves the differential equations and advances the solution to the next time step.

6. Presentation of results

This section presents some of the results of the numerical computation, primarily to study the effects of nonlinearity upon the evolving interface profiles. In addition, the effect of a finite-width interfacial zone, rather than a sharp interface, is also of interest. The linearized stability criterion sketched in Figure 2 gives an important guide as to the anticipated behaviour of solutions; furthermore, the contours of the growth function (3.21) in Figure 2 show precisely the relationship between mode number n and viscosity ratio $\beta $ that allows two different solutions to have the same growth rate $G_n$ . For simplicity, in this section, we therefore illustrate solutions at the one value $\beta = 0.5$ of the viscosity ratio and generate different growth rates by varying the initial mode number $m^{*}$ .

The first task is to establish the credibility of the nonlinear solution algorithm in Section 4 by comparing its predictions with those of the linearized solution from Section 3, since the two should be in good agreement if started from a perturbation of appropriately small amplitude. This comparison is illustrated in Figure 3 for a solution started from a pure fifth-mode perturbation ( $m^{*} = 5$ ) and one with a ninth-mode disturbance ( $m^{*} = 9$ ). In both cases, the initial amplitude had the small value $\epsilon = 0.05$ , and interface profiles are sketched at the two (dimensionless) times $t=1$ and $t=3$ . The linearized solution is obtained from (3.18) and drawn with a dashed line (red online), and the solid (blue) line represents the nonlinear result for the sharp-interface model.

Figure 3 A comparison of the linearized interface (dashed lines, red online) with the predictions of the nonlinear sharp-interface model of Section 4 (solid line, blue online) for $\beta = 0.5$ , $\sigma = 1.E-4$ , $\epsilon = 0.05$ . The nonlinear solution used $N=101$ Fourier coefficients, with Lanczos parameter $0.05$ . Results are shown for initial mode (a) $m^{*}=5$ and (b) $m^{*}=9$ .

Figure 4 Interfacial curvature $\kappa $ calculated for the nonlinear sharp-interface model of Section 4, for the same parameters as in Figure 3. The solution used $N=101$ Fourier coefficients, with Lanczos parameter $\lambda _T = 0.05$ . Results are shown for initial mode (a) $m^{*}=5$ and (b) $m^{*}=9$ .

At early times, there is good agreement between the linearized and nonlinear interface profiles, as would be expected. This is apparent in both the solutions at the earlier time $t = 1$ in both Figure 3(a) for $m^{*} = 5$ and in Figure 3(b) for $m^{*} = 9$ . For the $m^{*} = 5$ solution at $t= 1$ , the linearized and nonlinear interface profiles are nearly identical, and this is to be expected since, from Figure 2, the instability grows more slowly when $n=5$ than when $n=9$ . Figure 3(a) therefore helps establish the reliability of the numerical results in the nonlinear case. These were generated here using $N=101$ Fourier coefficients, and we have also verified that they are consistent with results at $N=81$ and $N=51$ , to a very high degree of accuracy. For the $m^{*} =9$ solution shown in Figure 3(b), there is still reasonable agreement between the linearized and nonlinear solutions at the earlier time $t=1$ , but differences become evident by the later time $t = 3$ . This is because the linearized solution (3.18), when started with the monochromatic disturbance $m^{*}=9$ , can only produce the single mode $n=9$ with an amplitude that must grow exponentially with time. By contrast, however, the nonlinear solution involves all the Fourier modes and, as a result, nonlinear saturation takes place, as energy is transferred to higher modes. Thus, the nonlinear disturbance by time $t=3$ has developed a nonsinusoidal shape with a much reduced overall wave amplitude.

Nonlinear curvatures of the interface are shown in Figure 4 for the same conditions $\beta = 0.5$ , $\sigma = 10^{-4}$ , $\epsilon = 0.05$ as in Figure 3, and for the three times $t=1$ , $2$ and $3$ . The curvatures were calculated from (2.11). The initial condition stipulates that, at $t=0$ , the interface is a circle of unit radius with an added perturbation; when that perturbation is small, then the linearization assumption (3.6) shows that the curvature behaves like

(6.1) $$ \begin{align} \kappa (\theta, t) = \frac{1}{\sqrt{1+t}} - \frac{\epsilon}{1+t} \bigg[ \frac{\partial^2 R_1}{\partial \theta^2} + R_1 \bigg] + \mathcal{O}( \epsilon^2). \end{align} $$

Thus, the mean curvature decreases with time like $\kappa \sim 1/\sqrt {t}$ , because the interface is being forced outwards by the injection of fluid at the origin; however, the linearized theory nevertheless predicts that the perturbations at the interface, along with their curvature, grow exponentially with time.

For the fifth-mode solution $m^{*} = 5$ in Figure 4(a), the curvature at the earliest time $t=1$ does indeed consist of an essentially sinusoidal disturbance to the mean curvature $1/\sqrt {2}$ , as predicted by (6.1). However, by time $t=3$ , the oscillations in curvature are noticeably nonsinusoidal, with sharply pointed troughs and broad crests, caused by the growing influence of nonlinear energy transfer to higher Fourier modes. A similar description applies also to the $m^{*} = 9$ curvatures displayed in Figure 4(b), except that the higher-mode perturbation in this case grows more rapidly than for the $m^{*} =5$ solution in Figure 4(a), and so the effects of nonlinearity are far more developed. At the last time $t=3$ shown, the curvature consists of nearly constant portions around the interface, punctuated by sudden, sharp downward spikes.

As the sharp-interface model of Section 4 is run further forward in time, nonlinear effects become more pronounced, as discussed in relation to Figure 4. The curvature of the interface, in particular, develops sharp spikes at selected points and, eventually, the entire model fails, due to the formation of curvature singularities. This was first predicted by Moore [22], and it occurs at some finite critical time that depends on the perturbation mode number and the initial disturbance amplitude. This is illustrated in Figure 5, where a ninth-mode solution, with initial mode number $m^{*} =9$ and initial amplitude $\epsilon = 0.05$ , has been allowed to evolve in time until failure occurred, in this instance, at a dimensionless time slightly larger than $t = 5.6$ . In Figure 5, we therefore show the nonlinear interface computed with $N=101$ coefficients at time $t=5.6$ (Figure 5(b)) and, for comparison, the diagram in Figure 5(a) shows the interface at $t=2$ . These nonlinear solutions are drawn with heavy solid lines (blue online) and, in addition, the linearized interface shape from Section 3 is indicated with dashed lines (red online).

Figure 5 A comparison of the linearized interface (dashed lines, red online) with the predictions of nonlinear sharp-interface model of Section 4 (solid line, blue online) for a ninth-mode perturbation ${m^{*} =9}$ and parameters $\beta = 0.5$ , $\sigma = 1\times 10^{-4}$ , $\epsilon = 0.05$ , at the two times (a) $t=2$ and (b) $t=5.6$ . The nonlinear solution used $N=101$ Fourier coefficients, here with Lanczos parameter $\lambda _T =0.02$ .

As expected from Figure 3, the agreement in Figure 5 between the linearized and nonlinear results at time $t=2$ is reasonably good. However, at time $t = 5.6$ , which is close to the failure time, the linearized solution substantially over-estimates the amplitude of the perturbations at the interface, because it assumes that all the energy remains confined within the monochromatic mode $n=9$ for all time. By contrast, the nonlinear interface develops high-curvature dimples and broad, flat tips that are reminiscent of the “tip-splitting” event illustrated in a schematic diagram in the paper by Kim et al. [18, Figure 4(b)].

The curvature is shown in Figure 6 for the same parameters as illustrated in Figure 5 and at the last time $t=5.6$ for which solutions could be computed. Although the mean curvature decreases slowly with time, as indicated in the linearized result (6.1), the profile in Figure 6 nevertheless also contains large positive and negative curvature spikes that punctuate almost flat portions in the profile where the curvature is roughly constant. This is entirely consistent with the nonlinear interface shape shown in Figure 5(b) at time $t=5.6$ . Our computer code failed to integrate significantly beyond this time, which is therefore anticipated to be close to the critical time at which points of infinite curvature form on the interface and the sharp-interface mode of Section 4 ceases to be valid.

Figure 6 The curvature for the nonlinear sharp-interface model of Section 4, for a ninth-mode perturbation $m^{*} =9$ and parameters $\beta = 0.5$ , $\sigma = 1\times 10^{-4}$ , $\epsilon = 0.05$ . The nonlinear solution used $N=101$ Fourier coefficients, here with Lanczos parameter $\lambda _T =0.02$ .

The single-fluid model in Section 5, with a diffuse interfacial zone of finite width across which some degree of fluid mixing can occur, was developed to overcome the formation of Moore curvature singularities [22] like those developing in Figure 6. Accordingly, we show in Figure 7 a comparison of the predictions of both nonlinear models, from Sections 4 and 5, at two different times. The background images in both cases are contours of the temperature T computed from the variable-viscosity model in Section 5 for the case $\beta =0.5$ , $\tau =1.5$ and diffusion constant $\sigma _H = 10^{-4}$ in the heat equation (5.1). This is a ninth-mode solution ( $m^{*} =9$ ) and was started from a perturbation with amplitude $\epsilon _T =0.05$ . The numerical results were run with $(M,N)=(101,105)$ Fourier coefficients and were subject to Lanczos smoothing at regular time intervals $\Delta t = 0.2$ using smoothing parameters $\lambda _T = \lambda _R = 0.05$ in the damping strategy (5.22). As outlined in Section 5.1, each new time step calculated by the solution algorithm requires the solution of a large matrix equation to obtain the coefficients $B_{m,n}$ . This is done iteratively and we have found that the use of eleven iterations at each time step is sufficient to reduce the iteration error in the solution of the equations (5.18) to below $\delta = 10^{-11}$ .

Figure 7 A comparison of the results of the nonlinear diffuse-interface approach of Section 5 with the predictions of the sharp-interface model of Section 4 (solid line, red online), for a ninth-mode perturbation $m^{*} =9$ with parameters $\beta = 0.5$ , $\sigma = 1\times 10^{-4}$ , $\epsilon _T = 0.05$ . The diffuse-interface solution used $\tau = 1.5$ .

Overlaid on Figure 7, in a heavy solid line (red online), is a nonlinear solution for a sharp interface for this same case $m^{*} =9$ , $\beta =0.5$ , with initial amplitude $\epsilon = 0.05$ . The surface tension in the dynamic boundary condition (2.9) has the value $\sigma = 10^{-4}$ , but it is important to observe that this is not the same as the diffusion parameter $\sigma _H = 10^{-4}$ also used in Figure 7 for the single-fluid model, even although both dimensionless parameters have the same numerical value. In this respect, then, the comparison between the two nonlinear theories in Figure 7 is not exact; nevertheless, it is about the best that can be done and the difference is not expected to detract greatly from the overall comparison.

For the diagram in Figure 7(a) at the earlier time $t=2$ , there is reasonably good overall agreement between the two nonlinear theories, despite clear differences in detailed features. In particular, the sharp interface from Section 4 lies very precisely over the crests of the outward-moving fingers of fluid 1 obtained with the diffuse model of Section 5. For the diagram in Figure 7(b), calculated from the two nonlinear theories at time $t=5.6$ , there is again reasonably close agreement between the two results at the crests of the fingers. Nevertheless, this is about the last time at which the sharp-interface model of Section 4 could yield a solution, due to the formation of curvature singularities as illustrated in Figures 5 and 6. The variable-viscosity model of Section 5, at this time, has formed fingers that have largely detached from the central region, within which fluid 1 continues to be injected at the origin. The heads of these fingers have a blunt shape that is very similar to that predicted by the nonlinear sharp-interface model, but they also possess pointed regions on each side of the finger-tips, forming small, overhanging regions. The variable-viscosity solution does not encounter a problem at $t = 5.6$ , but can continue to be integrated forward in time beyond this value, until issues of numerical accuracy eventually cause the integration package to stall at significantly later times.

Perhaps the most striking difference between the nonlinear sharp-interface model and the diffuse-interface model of Section 5 concerns the way in which the diffuse-interface results show outward-moving fingers that have separated almost completely from the central core. This often happens quite early and Figure 8 presents a sequence of solutions at eight early times, so as to illustrate the way in which detachment and separation of the fingers occurs. At the first time $t=0.1$ shown, the central core containing fluid 1 simply possesses nine little oscillations arising from the ninth-mode perturbation $m^{*} =9$ given to it initially. At some time close to $t \sim 1$ , the fingers can be seen to start the process of detaching from the central region, and by the last time $t=1.5$ shown, they have broken away almost completely.

Figure 8 Temperature profiles showing the detachment and shedding of fingers at eight early times, for a ninth-mode ( $m^{*} =9$ ) perturbation. Here, $\beta = 0.5$ and $\tau = 1.05$ . The temperature diffusion coefficient was $\sigma _H = 1\times 10^{-5}$ , the initial temperature amplitude was $\epsilon _T = 0.1$ and the algorithm was run with ${(M,N)=(81,85)}$ coefficients.

This is illustrated further in Figure 9, where the mode-nine solution from Figure 8 has been re-computed with $(M,N)=(101,105)$ Fourier coefficients and subject to a small amount of Lanczos smoothing with $\lambda _T = \lambda _R = 0.02$ . The vertical axis in this diagram represents the temperature, computed with the spectral representation (5.6) with temperature ratio $\tau = 1.05$ and starting from initial temperature amplitude $\epsilon _T = 0.1$ . At this later time, the fingers have effectively detached completely from the inner core region (fluid 1) and this region has consequently returned to a cylindrical shape. This is an interesting outcome since, although the linearized solution in Figure 2 indicates that the solution is unstable (at early times and for small finger amplitude), this is apparently not the case at very late times; instead, the solution begins as an unstable flow, but then the unstable fingers break off and move far away, leaving behind a circular cylinder of fluid 1 injected at the origin. This is now a stable geometry, with $n=0$ , and unless further disturbances are given to the system, it will remain a stable flow with a circular inner core that has radius $R_0 (t)$ growing slowly according to (3.3).

Figure 9 A temperature profile showing detached fingers at time $t=3.6$ for the ninth-mode ( $m^{*} =9$ ) perturbation in Figure 8. Here, $\beta = 0.5$ , $\sigma _H = 1\times 10^{-5}$ , $\tau = 1.05$ , $\epsilon _T = 0.1$ and the number of modes was increased to $(M,N)=(101,105)$ .

7. Conclusion

In this article, we have presented numerical solutions for two highly nonlinear models of Saffman–Taylor flow occurring in radial geometry; that is, an ambient viscous fluid was originally present in a porous medium and then a second fluid of lower viscosity was injected through a line source at the origin. This creates a circular interface between the two fluids that moves radially outwards. An initial perturbation is given to the otherwise circular interface and this perturbation grows unstably because the inner fluid is less viscous than the outer one.

The first model considered was one in which there was assumed to be an infinitesimally thin interface between the two viscous fluids, so that the viscosity jumps discontinuously as we cross the interface. A classical type of linearized solution to this sharp-interface model has been presented in Section 3, leading to the stability condition (3.19) from which it follows that high-mode perturbations are unstable if the viscosity ratio $\beta < 1$ so that the inner fluid is less viscous than the outer. Carrying out such a solution for radial geometry is somewhat novel in the context of Saffman–Taylor instability and Figure 2 summarizes the findings of the linearized solution in a compact form. Similar results, in other contexts, have also been presented by Grodzki and Szymczak [12] and Beeson-Jones and Woods [2].

The fully nonlinear thin-interface model was then solved numerically in Section 4, using an adaptation of the “basic” spectral method introduced by Forbes et al. [10]. Here, we have used an exact Fourier series-type solution (4.2) of Laplace’s equation for the two fluid pore pressures within the porous medium, taking account of the line source of fluid 1 at the origin. This solution involves Fourier coefficients that are unknown functions of time, and these must be found by imposing the exact nonlinear conditions on a moving interface $r = R(\theta ,t)$ of unknown shape and location. The structure of the equations that arise in this porous medium problem requires a matrix equation to be solved at each new time step and, at first, we found that the system was so ill-conditioned that it could not be done to satisfactory accuracy or with sufficient total number N of Fourier modes. This difficulty was overcome here by scaling the radial coordinate r with a time-dependent function $R_S (t)$ in the representations (4.2) for pressure. It was found that the simple circular radius function in (4.1) was sufficient to make the equations well conditioned, although eventually, a critical time is attained at which the interface function $r = R(\theta ,t)$ develops Moore curvature singularities [22] at certain points $\theta $ . This is a feature of interfacial fluid models that involve infinitesimally thin interfaces, as shown by Moore [22], but it is possible that a different choice of scaling function $R_S(t)$ than that used in Section 4 might allow numerical solutions to be computed more accurately close to the critical failure time. This may be worth future study.

To avoid the curvature singularities encountered with a sharp interface, a second model was investigated in Section 5, in which the viscosity varies smoothly in space and there is a narrow interfacial region across which the change in viscosity from one region to the other is rapid. This was accomplished here by allowing the viscosity to depend solely on temperature and adding a convection-diffusion energy equation (5.1) to the mathematical model of the system. We postulated the simplest possible “equation of state” for the fluid, consisting of a linear relation (5.4) between viscosity and temperature. It was shown here that this new variable-viscosity model agreed well with both the linearized and nonlinear results from the first sharp-interface model, at least at early times; furthermore, the variable-viscosity model continues to produce meaningful solutions at times significantly after the critical time when the sharp-interface model fails because of the curvature singularities it creates at the interface.

The two nonlinear theories produce outward-moving fingers that agree quite closely near their tips. However, the major difference is that the diffuse-interface model in Section 5 creates long fingers that soon detach from the central region in which fluid 1 is injected. These fingers continue to move outwards as time increases, leaving behind a circular central region that grows slowly with time according to the simple background shape (3.3). Of course, this only happens because the injection of fluid 1 at the origin occurs in the steady fashion (2.6) and, if this were replaced by a source strength that varied with time, it would be expected that further fingers would continue to form and detach as time progressed. In some studies, such as the recently published results of Kim et al. [19], long fingers with somewhat bulbous heads also form, similar to those found here, although most do not appear to detach from the central injection region in quite the same way as here. Those authors used a different set of model equations, based on a phase-field approach with the Cahn–Hilliard equation; in addition, their viscosity depended exponentially on the concentration of a dissolved solute rather than the temperature. Our model in Section 5 avoids the use of phase-field methods, since our intention was to mimic conditions in the sharp-interface model of Section 4 as much as possible, and these differences doubtlessly influence the details of the solutions that are obtained. It would be of interest to consider more complicated relations between viscosity and temperature than the simple linear formula (5.4), and in addition to study the effects of dissolved solutes, such as in Kim et al. [19], Trevelyan et al. [28] and Forbes et al. [9]. This work is on-going.