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 [Reference Grodzki and Szymczak12] 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. [Reference Forbes, Chen and Trenham10], extended to planar cylindrical geometry by Forbes [Reference Forbes8].

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 [Reference von Winckel29].

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 [Reference Kreyszig20]) 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 [Reference Duchon6]. 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 [Reference Moore22]. 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 [Reference Abramowitz and Stegun1, 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 [Reference Abramowitz and Stegun1, 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. [Reference Walters, Turner and Forbes30] and Forbes et al. [Reference Forbes, Walters and Turner11]. 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 [Reference Kreyszig20]. We avoid this difficulty here using Lanczos smoothing [Reference Duchon6] 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 [Reference Moore22], 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. [Reference Kim, Funada, Joseph and Homsy18, 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 [Reference Moore22] 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)$ .