FULLY 3D FLUID OUTFLOW FROM A SPHERICAL SOURCE

Abstract

We consider fully three-dimensional time-dependent outflow from a source into a surrounding fluid of different density. The source is distributed over a sphere of finite radius. The nonlinear problem is formulated using a spectral approach in which two streamfunctions and the density are represented as a Fourier-type series with time-dependent coefficients that must be calculated. Linearized theories are also discussed and an approximate stability condition for early stages in the outflow is derived. Nonlinear solutions are presented and different outflow shapes adopted by the fluid interface are investigated.

1 Introduction

When a compact source of incompressible fluid is placed within a second incompressible fluid, then at some appropriate distance away, the outflow often behaves as if the source were simply a mathematical point singularity. Thus, if the source produces fluid at some rate $Q_S$ (volume per time), then at some distance r from the source, there will be a radial outflow of (approximate) speed $Q_S / ( 4\pi r^2 )$ . However, if the fluid issued from the source differs in density from the ambient fluid, then some type of interfacial region will be formed between the two and will move outwards with some shape which might not necessarily be a simple sphere.

Modelling the idealized problem of a point source embedded in a different ambient fluid, in the presence of an interface between the fluids, is in a sense a fundamental problem in free-surface hydrodynamics. Although it is simple to conceptualize and can be modelled with a straightforward system of well-known equations, it is nevertheless highly nonlinear, due to the presence of the interface, so its solutions are difficult to calculate and its behaviour may be unexpected. This outflow problem merits study in its own right and is also related to a range of situations of scientific interest, particularly in astrophysics and geophysics where there may also be a gravitational field directed inwards towards the source. In that case, the situation is one in which a light fluid is driven, by the source, out into a heavier ambient fluid, against the inward acceleration of gravity. Thus, it constitutes a type of spherical Rayleigh–Taylor instability.

In the original problem of Rayleigh [23] and Taylor [25], two horizontal sheets of inviscid fluid were pictured, separated by an infinitesimally thin flat interface. A density jump was present at the interface, and the entire fluid system was subject to the acceleration of gravity. A small perturbation to the interface is stable if the lower fluid is the heavier of the two. However, if the upper fluid is the heavier, then the system experiences the Rayleigh–Taylor (RT) instability and the small disturbance grows rapidly. In the linearized theory of Rayleigh and Taylor, in which disturbances to the interface are supposed to be small, the theory nevertheless predicts that the unstable growth occurs exponentially rapidly; clearly, the linearized theory fails to be valid within some finite time and nonlinear effects then dominate. An enormous volume of literature now exists on this unstable flow, and the review article by Andrews and Dalziel [3] identifies three stages of the development of RT unstable flow for density ratios close to one. In the first stage, for early times, the growth is indeed exponential, as predicted by the linearized theory discussed above. However, a second stage in RT mixing exists, in which the amplitude of the disturbance saturates, and instead, Fourier modes involving longer wavelengths are entrained. Andrews and Dalziel [3] also identify a third stage in RT flows, in which the growing bubble tip becomes self-similar. A very detailed review of Rayleigh–Taylor and related instabilities is given by Zhou [28, 29]. Furthermore, the RT instability can occur over distances that differ across hundreds of orders of magnitude, as indicated by Kelley et al. [16].

The RT instabilities can occur in geometries other than the original planar case studied by Rayleigh [23] and Taylor [25]. In addition, fully three-dimensional RT instabilities can also occur on a surface, and these are reviewed by Zhou [29]. A question of interest concerns the extent to which initial conditions affect the flow properties at later times and the degree to which memory of initial perturbations might be lost if, indeed, RT outflow is ultimately dominated by the formation of self-similar structures. Zhou [29] presents a detailed discussion of this point, and reviews the differing opinions and the extensive experimental and numerical findings underpinning them. Forbes [9] considered a cylindrical RT instability, in which a line source lying along the z-axis produced an outflow of less-dense fluid into a denser ambient fluid, and investigated the growth of modal periodic disturbances around the otherwise cylindrical interface. He found a strong dependence on initial conditions and observed that a mode-n perturbation on the cylinder would ultimately lead to the formation of n jets of light fluid around the cylinder, each of which formed an overturning mushroom-shaped structure at its tip. Forbes [10] then considered the analogous problem in spherical geometry, in which a point source of less-dense fluid was embedded in a heavier ambient fluid, and there he obtained a very different outcome. He undertook a linearized stability analysis of small perturbations to the spherical interface and obtained a stability equation that had been derived earlier by Mikaelian [19]. He demonstrated that this stability equation, for many choices of parameters, predicted that the first spherical mode had the highest growth rate and was therefore the most unstable. Forbes [10] argued that the mode-coupling effects of nonlinearity would eventually result in the first mode being excited, regardless of initial conditions, so that the long-term disturbance would be dominated by the first spherical mode, giving a one-sided jet emerging from the sphere. This was indeed observed in many of Forbes’ numerical solutions of a viscous model based on the Boussinesq equations. Later, Forbes and Brideson [12] and Forbes [11] investigated the extent to which straining and rotation of the outer ambient fluid could alter these one-sided morphologies. Nevertheless, one-sided outflows in astrophysical circumstances have been predicted numerically (Lovelace et al. [18]) and observed experimentally by Gómez et al. [14] and Bietenholz et al. [6] (who also found that the outflow region behaved as a point source in their measurements).

This key difference between planar and spherical RT outflows is believed to be of importance in geophysical and astrophysical applications. In an early model of motion within the Earth’s core, Jeffreys and Bland [15] considered linearized equations describing flow within a sphere of viscous self-gravitating fluid, heated at its centre. They found that the lowest spherical mode $n=1$ is indeed the most unstable, so that asymmetrical features could develop within the Earth, and they saw this as a possible explanation for the presence of continents at the Earth’s surface. Chandrasekhar [7] developed this theme by studying flow of a (highly) viscous fluid within a sphere, in which the fluid is incompressible and subject to radially directed acceleration. He found that for a fluid of sufficiently high viscosity, maximum instability again occurs at the first spherical mode $n=1$ . It is interesting that Chandrasekhar’s [7] conclusion, concerning the $n=1$ mode, is qualitatively similar to that obtained by Forbes [10] whose linearized analysis ignored viscosity altogether. More recently, Mondal and Korenaga [20] have presented an analysis of the full dispersion relation for the linearized RT instability of a two-fluid Newtonian viscous fluid system in spherical geometry, and Terrones and Heberling [26] contrast the behaviour of such a spherical arrangement of fluids with that of the corresponding planar flow. They also demonstrate that the first spherical mode $n=1$ is often the most unstable, even when the two fluids have extremely different viscosities and the inner “core” fluid is allowed to be of arbitrarily large viscosity. A recent paper by Oren and Terrones [21] has indicated that the presence of an outer boundary can affect the stability of the system, as predicted by linearized theory, however, and for a finite-depth outer fluid, the lowest mode $n=1$ may no longer necessarily be the most unstable. These results have obvious relevance to geophysical modelling of the Earth’s mantle and core.

Approximating the outflow region as a mathematical point source poses problems for any numerical solution scheme, since this creates a singularity within the flow. In the axisymmetric outflows studied by Forbes [10, 11] based on the use of a single streamfunction, this singularity was able to be incorporated directly into a numerical spectral representation of the streamfunction and the fluid velocity components, and the effects of the singularity at the source could be dealt with in a fairly straightforward manner. In fully three-dimensional flow, however, a minimum of two streamfunctions is required, as will be discussed in Section 3, and the effects of the singularity at the point source turn out to be much more difficult to ameliorate. Consequently, we consider a finite-sized source region, over which the fluid speeds remain finite, but nevertheless produces a finite flux, as for the point source. This formulation is discussed in Section 2 and an approximate stability criterion for a source of finite radius is derived in Section 5. This paper therefore aims to contribute to the discussion concerning the influence of initial conditions, by taking into account two extra factors, namely, the effects of finite-sized source regions and the influence of fully three-dimensional geometry upon the eventual outflow morphologies. This is a computationally intensive aim, and we have made use of a high-performance computer that can accommodate large-scale parallelization of our solution algorithm on a graphics processor unit. Some results are displayed in Section 7. A discussion in Section 8 concludes this paper.

2 Boussinesq outflow model

We consider a finite-sized fluid source distributed over some sphere of radius $R_S$ . Its total flux is $Q_S$ (volume per time) and it also possesses a total mass $M_S$ . Outside that source radius, the fluid is regarded as incompressible and is viscous, with dynamic viscosity $\mu $ . We therefore postulate a background flow, produced by the distributed source, of the form

(2.1) $$ \begin{align} \begin{cases} \displaystyle \mathbf{q}_S = \frac {Q_S \mathbf{r}}{4\pi r^3} & \mbox{for} \ r> R_S\\[10pt] \displaystyle \mathbf{q}_S = \frac {Q_S r \mathbf{r}}{4\pi R_S^4} & \mbox{for} \ r < R_S. \end{cases} \end{align} $$

Here, the source velocity vector is denoted $\mathbf {q}_S$ . The vector $\mathbf {r} = x\mathbf {i} + y\mathbf {j} + z\mathbf {k}$ is the usual position vector of a point measured from the origin of a Cartesian coordinate system (positioned at the centre of the source sphere) and $r =\! \sqrt{x^2 + y^2 + z^2}$ is its length. Since the divergence of the vector $\mathbf {q}_S$ in (2.1) is positive within the source sphere $r < R_S$ , then mass is being created in that region, as desired. This background velocity (2.1) is directed purely radially, and is continuous at the edge surface $r = R_S$ . The mass density of the fluid is denoted $\rho ( x,y,z,t )$ and is assumed to be constant at the same reference value $\rho _1$ inside the source sphere $r < R_S$ . By integrating the divergence $\textrm {div} \mathbf {q}_S$ of the source velocity in (2.1) over the source sphere, it is straightforward to show that the mass outflow per time generated by this source is $\rho _1 Q_S$ , as is to be expected. Finally, we suppose that at the initial time $t = 0$ , the inner fluid with density $\rho _1$ occupies a spherical region $r < a$ which completely contains the distributed source $r < R_S$ so that, necessarily, $R_S < a$ . There is an interface at $r = a$ and some outer ambient fluid with density $\rho _2$ beyond the interface in the domain $r> a$ .

In the Boussinesq approximation, the two-fluid medium with its interface between the fluids is modelled as a single fluid with a density that is a smoothly varying continuous function. The interface is therefore replaced by a narrow zone in which the density changes smoothly and rapidly from $\rho _1$ in the inner fluid to $\rho _2$ in the outer layer. Furthermore, it is assumed that the density variation is only slight, so that the total density function can be written as

$$ \begin{align*} \rho ( x,y,z,t ) = \rho_1 + \overline{\rho} ( x,y,z,t ) , \end{align*} $$

in which the perturbation function $\overline {\rho }$ is small relative to the density $\rho _1$ of the inner fluid, which is now used as a reference quantity.

Dimensionless variables are now introduced. All lengths are scaled relative to the radius a of the initial sphere of inner fluid. Rather than use the source flux $Q_S$ in (2.1) as a scaling quantity, we instead introduce some arbitrary flux $Q_0$ (with the dimensions of volume per time), since this then allows the source to be removed simply by setting $Q_S = 0$ if desired. Speeds are nondimensionalized relative to the quantity $Q_0 / a^2$ and time is measured as a fraction of $a^3 / Q_0$ . The density $\rho $ and pressure p are scaled using $\rho _1$ and $\rho _1 Q_0^2 / a^4$ , respectively. It will be seen that the solution to this problem will be governed by the six dimensionless constants:

(2.2) $$ \begin{align} \begin{array}{lll} f_{\ell} = \dfrac {Q_S}{Q_0} ,& D = \dfrac {\rho_2}{\rho_1},& r_S = \dfrac {R_S}{a} , \\[12pt] R_e = \dfrac {\rho_1 Q_0}{\mu a} ,& F^2 = \dfrac {Q_0^2}{G \mathcal{M} a^3} ,& \sigma = \dfrac{K a}{Q_0}. \end{array} \end{align} $$

Clearly, the first parameter $f_{\ell }$ is the source flux made dimensionless with respect to the (arbitrary) flux $Q_0$ ; in practice, this quantity would simply be one unless there is no source, when it would become zero. The density ratio is D and the quantity $r_S$ gives the nondimensional radius of the source sphere, and it is assumed that $r_S < 1$ so that the source is wholly contained within the inner fluid. The Reynolds number for this viscous flow is $R_e$ and it varies inversely with the fluid viscosity $\mu $ . For the results to be presented in this paper, the dynamic viscosity $\mu $ is taken to be constant throughout the fluid domain so that both the source and the ambient fluids have the same viscosity. We have also defined a Froude number F based on the gravitational constant G and the mass $\mathcal {M}$ of the distributed source, and finally, it is convenient to introduce a parameter $\sigma $ that measures the density diffusion rate K in a dimensionless manner.

In these dimensionless coordinates, it is convenient to represent the velocity vector $\mathbf {q}$ as the sum

(2.3) $$ \begin{align} \mathbf{q} = \mathbf{q}_S + \mathbf{Q}, \end{align} $$

in which the source velocity is obtained from (2.1) in the dimensionless form

(2.4) $$ \begin{align} \begin{cases} \displaystyle \mathbf{q}_S = \frac {f_{\ell} \mathbf{r}}{4\pi r^3} \quad \mbox{for} \ r> r_S\\[8pt] \displaystyle \mathbf{q}_S = \frac {f_{\ell} r \mathbf{r}}{4\pi r_S^4} \quad \mbox{for} \ r < r_S, \end{cases} \end{align} $$

and the correction velocity $\mathbf {Q}$ is incompressible, so that

(2.5) $$ \begin{align} \textrm{div}\, \mathbf{Q} = 0. \end{align} $$

It is convenient to represent the source velocity and the correction velocity vectors in terms of their Cartesian components and accordingly,

$$ \begin{align*} \mathbf{q_S} & = u_S \mathbf{i} + v_S \mathbf{j} + w_S \mathbf{k}, \\ \mathbf{Q} & = U \mathbf{i} + V \mathbf{j} + W \mathbf{k}. \end{align*} $$

In Boussinesq theory, the density perturbation function $\overline {\rho }$ satisfies the transport equation

(2.6) $$ \begin{align} \frac {\partial \overline{\rho}}{\partial t} + \mathbf{q} \cdot \nabla \overline{\rho} = \sigma \nabla^2 \overline{\rho} , \end{align} $$

and the Navier–Stokes–Boussinesq equation that approximates the conservation of linear momentum in this theory takes the form

(2.7) $$ \begin{align} \frac {\partial \mathbf{q}}{\partial t} + ( \mathbf{q} \cdot \nabla ) \mathbf{q} + \nabla p = ( 1 + \overline{\rho} ) \mathbf{G}_S + \frac {1}{Re} \nabla^2 \mathbf{q}. \end{align} $$

In this expression, $\mathbf {G}_S$ represents the gravitational force per mass due to the distributed source (of dimensional mass $\mathcal {M}$ , as in (2.2)). To simplify matters, we have approximated this gravitational acceleration with the form

(2.8) $$ \begin{align} \begin{cases} \displaystyle \mathbf{G}_S = - \frac {1}{F^2} \frac {\mathbf{r}}{r^3} \quad \mbox{for}\ r> R_S\\[8pt] \displaystyle \mathbf{G}_S = - \frac {1}{F^2} \frac {r \mathbf{r}}{r_S^4} \quad \mbox{for} \ r < R_S, \end{cases} \end{align} $$

mimicking to some extent the form (2.4) assumed for the outflow velocity vector produced by the source.

We now follow Forbes and Brideson [13] and Walters and Forbes [27], and satisfy the requirement (2.5) identically, using two streamfunctions $\Psi ( x,y,z,t )$ and $\chi ( x,y,z,t )$ and expressing the correction velocity vector as $\mathbf {Q} = \textrm {curl} \, \mathbf {A}$ with vector potential $\mathbf {A} = \Psi \mathbf {i} + \chi \mathbf {j}$ . Thus, the three correction velocity components can be written as

(2.9) $$ \begin{align} U = - \frac {\partial \chi}{\partial z} , \quad V = \frac {\partial \Psi}{\partial z} , \quad W = \frac {\partial \chi}{\partial x} - \frac {\partial \Psi}{\partial y}. \end{align} $$

The pressure p is now removed from further consideration by taking the vector curl of the momentum equation (2.7). We define the vorticity vector $\boldsymbol {\zeta } = \textrm {curl}\, \mathbf {q}$ , and making use of (2.4) shows that

(2.10) $$ \begin{align} \boldsymbol{\zeta} = \textrm{curl} \, \mathbf{q} = \textrm{curl} \, \mathbf{Q} \equiv \zeta^X \mathbf{i} + \zeta^Y \mathbf{j} + \zeta^Z \mathbf{k}. \end{align} $$

Now the curl operator applied to the Navier–Stokes–Boussinesq equation (2.7) gives rise to the Boussinesq vorticity equation

(2.11) $$ \begin{align} \frac {\partial \boldsymbol{\zeta}}{\partial t} + ( \textrm{div} \, \mathbf{q} ) \boldsymbol{\zeta} + (\mathbf{q} \cdot \nabla ) \boldsymbol{\zeta} - ( \boldsymbol{\zeta} \cdot \nabla ) \mathbf{q} = \nabla \overline{\rho} \times \mathbf{G}_S + \frac {1}{Re} \nabla^2 \boldsymbol{\zeta}. \end{align} $$

Although this vorticity equation (2.11) has three components, it turns out that only two of these are independent; if the first component is differentiated with respect to x and the second with respect to y and these are added, the result is the z-derivative of the third component. This requires a significant amount of algebra to accomplish, but it is nevertheless true because of the vector identity

$$ \begin{align*} \textrm{div}\, \boldsymbol{\zeta} = \textrm{div}\, \textrm{curl}\, \mathbf{Q} = 0. \end{align*} $$

Thus the vorticity equation (2.11) reduces to a system of two partial differential equations essentially for the two streamfunctions $\Psi $ and $\chi $ . The transport equation (2.6) then represents a third equation to be solved for the density perturbation function $\overline {\rho }$ .

3 Numerical solution technique

The numerical solution of a nonlinear system of fluid-mechanical equations in three space and one time variable is a demanding, computer-intensive task, and can only succeed if the numerical algorithm is both efficient and parallelizable. Originally, we intended to use spherical polar coordinates as in Forbes [10], but the coordinate singularity at the origin, in these coordinates, quickly renders the numerical algorithm too inefficient in three-dimensional space. Accordingly, we have instead elected to use three-dimensional Cartesian coordinates and base our numerical approach around a spectral representation of the fluid variables in this simpler geometry. A Cartesian-based technique has also been used by Lin et al. [17] to model the motion of the core and mantle of planets, although their method is based around the use of a finite-difference package. In addition, they consider two-dimensional geometry only and neglect the inertial terms in the equations of motion, thus effectively solving a linearized approximate model.

Therefore, to solve our system of three nonlinear partial differential equations, we first impose a cubical Cartesian “computational box” over the domain $-L < x < L$ , $-L < y < L$ , $-L < z < L$ centred at the origin. Since the inner fluid initially lies in a sphere of dimensionless radius $1$ and the distributed source is contained within that sphere, we require $ r_S < 1 < L$ . The two streamfunctions in (2.9) are now represented spectrally as

(3.1) $$ \begin{align} \begin{aligned} \chi ( x,y,z,t ) &= \displaystyle\sum_{m=1}^M \sum_{n=1}^N \sum_{p=1}^P A_{mnp} (t) \cos ( \alpha_m (x + L) ) \cos ( \alpha_n (y + L) ) \\ & \quad \times \sin ( \alpha_p (z + L) ), \\[.2cm] \Psi ( x,y,z,t ) & = \displaystyle \sum_{m=1}^M \sum_{n=1}^N \sum_{p=1}^P B_{mnp} (t) \sin ( \alpha_m (x + L) ) \sin ( \alpha_n (y + L) ) \\ & \quad \times \sin ( \alpha_p (z + L) ), \end{aligned} \end{align} $$

with

(3.2) $$ \begin{align} \alpha_m = \frac {m\pi}{2 L} ,\quad \alpha_n = \frac {n\pi}{2 L} , \quad \alpha_p = \frac {p\pi}{2 L}. \end{align} $$

To aid computational efficiency, the numbers of Fourier modes in each variable are taken to be equal, $M = N = P$ , so that only the one vector $\alpha $ needs to be calculated, to accommodate all the terms in (3.2) required to evaluate the series (3.1). An appropriate spectral representation for the density perturbation function $\overline {\rho }$ can be difficult to find since it must be completely general. After some experimentation, it was decided to choose

(3.3) $$ \begin{align} \begin{aligned} \overline{\rho} ( x,y,z,t ) &= \sum_{m=0}^M \sum_{n=0}^N \sum_{p=0}^P C_{mnp} (t) \cos ( \alpha_m (x + L) ) \cos ( \alpha_n (y + L) )\\ & \quad \times \cos ( \alpha_p (z + L) ). \end{aligned} \end{align} $$

From (2.9) and (3.1), the three correction velocity components now become

(3.4) $$ \begin{align} \begin{aligned} U ( x,y,z,t ) = & - \sum_{m=1}^M \sum_{n=1}^N \sum_{p=1}^P \alpha_p A_{mnp} (t) \cos ( \alpha_m (x + L) ) \cos ( \alpha_n (y + L) ) \\ & \times \cos ( \alpha_p (z + L) ), \\[.2cm] V ( x,y,z,t ) = & \sum_{m=1}^M \sum_{n=1}^N \sum_{p=1}^P \alpha_p B_{mnp} (t) \sin ( \alpha_m (x + L) ) \sin ( \alpha_n (y + L) ) \\ & \times \cos ( \alpha_p (z + L) ) \end{aligned} \end{align} $$

and

(3.5) $$ \begin{align} W ( x,y,z,t ) & = - \sum_{m=1}^M \sum_{n=1}^N \sum_{p=1}^P [\alpha_m A_{mnp} (t) + \alpha_n B_{mnp} (t) ] \sin ( \alpha_m (x + L) ) \nonumber \\ & \quad \quad \times \cos ( \alpha_n (y + L) ) \sin ( \alpha_p (z + L) ). \end{align} $$

The three components of the vorticity vector in (2.10) can also be obtained by straightforward differentiation. This gives

(3.6) $$ \begin{align} \zeta^X ( x,y,z,t ) & = \sum_{m=1}^M \sum_{n=1}^N \sum_{p=1}^P Z^X_{mnp} (t) \sin ( \alpha_m (x + L) ) \sin ( \alpha_n (y + L) ) \sin ( \alpha_p (z + L) ), \nonumber \\[.2cm] \zeta^Y ( x,y,z,t ) & = \sum_{m=1}^M \sum_{n=1}^N \sum_{p=1}^P Z^Y_{mnp} (t) \cos ( \alpha_m (x + L) ) \cos ( \alpha_n (y + L) ) \sin ( \alpha_p (z + L) ), \nonumber\\[.2cm] \zeta^Z ( x,y,z,t ) & = \sum_{m=1}^M \sum_{n=1}^N \sum_{p=1}^P Z^Z_{mnp} (t) \cos ( \alpha_m (x + L) ) \sin ( \alpha_n (y + L) ) \cos ( \alpha_p (z + L) ), \end{align} $$

in which we have defined the auxiliary Fourier coefficients

(3.7) $$ \begin{align} \begin{aligned} Z^X_{mnp} (t) & = \alpha_m \alpha_n A_{mnp} (t) + ( \alpha_n^2 + \alpha_p^2 ) B_{mnp} (t), \\[3pt] Z^Y_{mnp} (t) & = \alpha_m \alpha_n B_{mnp} (t) + ( \alpha_m^2 + \alpha_p^2 ) A_{mnp} (t), \\[3pt] Z^Z_{mnp} (t) & = \alpha_p ( \alpha_m B_{mnp} (t) - \alpha_n A_{mnp} (t) ). \end{aligned} \end{align} $$

The components of the fluid velocity vector are obtained from (2.3) by adding the contributions $( u_S , v_S , w_S )$ from the distributed source to the correction components calculated in (3.4) and (3.5).

The first two components of the vorticity equation (2.11) are now subjected to Fourier analysis. From the computational viewpoint, it is convenient to write these in the forms

(3.8) $$ \begin{align} \begin{aligned} \frac {\partial \zeta^X}{\partial t} & = \mathcal{N}^X + \frac {1}{R_e} \bigg( \frac {\partial^2 \zeta^X}{\partial x^2} + \frac {\partial^2 \zeta^X}{\partial y^2} + \frac {\partial^2 \zeta^X}{\partial z^2} \bigg), \\[.2cm] \frac {\partial \zeta^Y}{\partial t} & = \mathcal{N}^Y + \frac {1}{R_e} \bigg( \frac {\partial^2 \zeta^Y}{\partial x^2} + \frac {\partial^2 \zeta^Y}{\partial y^2} + \frac {\partial^2 \zeta^Y}{\partial z^2} \bigg) , \end{aligned} \end{align} $$

in which the nonlinear convective terms have been collected in the two auxiliary functions

(3.9) $$ \begin{align} \begin{aligned} \mathcal{N}^X & \equiv - ( \textrm{div}\, \mathbf{q} ) \zeta^X + \bigg( \zeta^X \frac {\partial u}{\partial x} - u \frac {\partial \zeta^X}{\partial x} \bigg) \\ & \quad + \bigg( \zeta^Y \frac {\partial u}{\partial y} - v \frac {\partial \zeta^X}{\partial y} \bigg) + \bigg( \zeta^Z \frac {\partial u}{\partial z} - w \frac {\partial \zeta^X}{\partial z} \bigg) + \Gamma^X, \\[.2cm] \mathcal{N}^Y & \equiv - ( \textrm{div}\, \mathbf{q} ) \zeta^Y + \bigg( \zeta^X \frac {\partial v}{\partial x} - u \frac {\partial \zeta^Y}{\partial x} \bigg) \\ & \quad + \bigg( \zeta^Y \frac {\partial v}{\partial y} - v \frac {\partial \zeta^Y}{\partial y} \bigg) + \bigg( \zeta^Z \frac {\partial v}{\partial z} - w \frac {\partial \zeta^Y}{\partial z} \bigg) + \Gamma^Y. \end{aligned} \end{align} $$

The third component of the vorticity equation (2.11) is linearly dependent on the first two components in (3.8) and so adds no new information. In these expressions, the quantity $\textrm {div}\, \mathbf {q}$ is calculated from (2.4) and (2.5). It becomes

$$ \begin{align*} \begin{cases} \displaystyle \textrm{div}\, \mathbf{q} = 0 \quad \mbox{for} \ r> r_S\\[8pt] \displaystyle \textrm{div}\, \mathbf{q} = \frac {f_{\ell}}{\pi r_S^4} \sqrt{x^2 + y^2 + z^2} \quad \mbox{for} \ r < r_S. \end{cases} \end{align*} $$

We have also defined the vector

$$ \begin{align*} \Gamma^X \mathbf{i} + \Gamma^Y \mathbf{j} + \Gamma^Z \mathbf{k} = \nabla \overline{\rho} \times \mathbf{G}_S \end{align*} $$

with the gravitational acceleration $\mathbf {G}_S$ defined in (2.8). The first component of this vector is

$$ \begin{align*} \begin{cases} \displaystyle \Gamma^X = - \frac {1}{F^2} \frac {z ( \partial \overline{\rho} / \partial y ) - y ( \partial \overline{\rho} / \partial z )}{[ x^2 + y^2 + z^2 ]^{3/2}} \quad \mbox{for} \ r> r_S,\\[8pt] \displaystyle \Gamma^X = - \frac {1}{F^2 r_S^4} [z ( \partial \overline{\rho} / \partial y ) - y ( \partial \overline{\rho} / \partial z ) ] \sqrt{x^2 + y^2 + z^2} \quad \mbox{for} \ r < r_S, \end{cases} \end{align*} $$

and the other two components are obtained similarly. The first equation in the system (3.8) is multiplied by the basis function

$$ \begin{align*} \sin ( \alpha_k (x + L) ) \sin ( \alpha_{\ell} (y + L) ) \sin ( \alpha_q (z + L) ), \end{align*} $$

and integrated over the “computational cube” $-L < x,y,z < L$ . This results in the system of differential equations

(3.10) $$ \begin{align} \frac {d Z^X_{k \ell q} (t)}{dt} &= - \frac {1}{R_e} \Delta_{k \ell q}^2 Z^X_{k \ell q} (t) \nonumber \\ & \quad + \frac {1}{L^3} \int_{-L}^L \int_{-L}^L \int_{-L}^L \mathcal{N}^X \sin ( \alpha_k (x + L) ) \sin ( \alpha_{\ell} (y + L) ) \sin ( \alpha_q (z + L) ) \, dx \, dy \, dz , \nonumber \\ k &= 1 , \ldots , M , \quad \ell = 1 , \ldots , N , \quad q = 1 , \ldots , P. \end{align} $$

Similarly, the second of the equations in the system (3.8) is multiplied by the basis function

$$ \begin{align*} \cos ( \alpha_k (x + L) ) \cos ( \alpha_{\ell} (y + L) ) \sin ( \alpha_q (z + L)), \end{align*} $$

and integrated to give

(3.11) $$ \begin{align} \frac {d Z^Y_{k \ell q} (t)}{dt} &= - \frac {1}{R_e} \Delta_{k \ell q}^2 Z^Y_{k \ell q} (t) \nonumber \\ & \quad + \frac {1}{L^3} \int_{-L}^L \int_{-L}^L \int_{-L}^L \mathcal{N}^Y \cos ( \alpha_k (x + L) ) \cos ( \alpha_{\ell} (y + L) ) \sin ( \alpha_q (z + L) ) \, dx \, dy \, dz , \nonumber \\ k &= 1 , \ldots , M , \quad \ell = 1 , \ldots , N , \quad q = 1 , \ldots , P. \end{align} $$

In these equations, the Fourier coefficients $Z^X_{k \ell q} (t)$ and $Z^Y_{k \ell q} (t)$ for the vorticity vector are those defined in (3.7), and it is also convenient to define constants

(3.12) $$ \begin{align} \Delta_{k \ell q}^2 = \alpha_k^2 + \alpha_{\ell}^2 + \alpha_q^2. \end{align} $$

Systems of differential equations for the original Fourier coefficients $A_{k \ell q} (t)$ and $B_{k \ell q} (t)$ are recovered from the differential equations in (3.10), (3.11) by differentiating and rearranging the first two equations in (3.7). This yields

(3.13) $$ \begin{align} \begin{aligned} \frac {d A_{k \ell q} (t)}{dt} & = \frac {1}{\alpha_q^2 \Delta_{k \ell q}^2} \bigg[ ( \alpha_{\ell}^2 + \alpha_q^2 ) \frac {d Z^Y_{k \ell q} (t)}{dt} - \alpha_k \alpha_{\ell} \frac {d Z^X_{k \ell q} (t)}{dt} \bigg], \\ \frac {d B_{k \ell q} (t)}{dt} & = \frac {1}{\alpha_q^2 \Delta_{k \ell q}^2} \bigg[ ( \alpha_k^2 + \alpha_q^2 ) \frac {d Z^X_{k \ell q} (t)}{dt} - \alpha_k \alpha_{\ell} \frac {d Z^Y_{k \ell q} (t)}{dt} \bigg]. \end{aligned} \end{align} $$

It remains now to undertake a similar Fourier analysis of the transport equation (2.6). It is first multiplied by the basis function

$$ \begin{align*} \cos ( \alpha_k (x + L) ) \cos ( \alpha_{\ell} (y + L) ) \cos ( \alpha_q (z + L) ) \end{align*} $$

with $k = 0,1, \ldots , M$ , $\ell = 0,1, \ldots , N$ and $q = 0,1, \ldots , P$ , and integrated over the “computational cube,” as previously. After some algebra we obtain

(3.14) $$ \begin{align} \frac {d C_{k \ell q} (t)}{dt} &= - \sigma \Delta_{k \ell q}^2 C_{k \ell q} (t) \nonumber \\ & \quad - \frac {1}{L^3 [ 1 + \delta_{k,0} ] [ 1 + \delta_{\ell ,0} ] [ 1 + \delta_{q,0} ]} \int_{-L}^L \int_{-L}^L \int_{-L}^L \mathcal{N}^R \cos ( \alpha_k (x + L) ) \nonumber \\[.1cm] & \quad \times \cos ( \alpha_{\ell} (y + L) ) \cos( \alpha_q (z + L) ) \, dx \, dy \, dz , \nonumber \\ k &= 0, 1 , \ldots , M , \quad \ell = 0, 1 , \ldots , N , \quad q = 0, 1 , \ldots , P. \end{align} $$

In this expression, the convective nonlinear terms in (3.3) have been combined into the single function

(3.15) $$ \begin{align} \mathcal{N}^R = u \frac {\partial \overline{\rho}}{\partial x} + v \frac {\partial \overline{\rho}}{\partial y} + w \frac {\partial \overline{\rho}}{\partial z} , \end{align} $$

following the notation used in (3.9). The constants $\Delta _{k \ell q}^2$ are those defined previously in (3.12). In addition, the quantity $\delta _{0,k}$ is the usual Kronecker delta symbol and takes the value $1$ if $k = 0$ , but zero otherwise.

The expressions (3.13) and (3.14) represent a coupled system of $3 M N P$ ordinary differential equations for the three sets of Fourier coefficients $A_{mnp} (t)$ , $B_{mnp} (t)$ , $C_{mnp} (t)$ . These can now be integrated forward in time from some suitable initial state. Once these coefficients have been determined, all the solution variables can be reconstructed from their spectral representations (3.1)–(3.7). However, since the source has been distributed over the sphere $r < r_S$ , it is required to make the perturbation density $\overline {\rho }$ zero in that region. Therefore, at every new time step in the numerical integration process, the function $\overline {\rho } ( x,y,z; t )$ is first created from coefficients $C_{mnp} (t)$ and its Fourier representation (3.3), and then replaced with the new perturbation density function

(3.16) $$ \begin{align} \begin{cases} \displaystyle {\overline{\rho}}^{*} ( x,y,z; t ) = \overline{\rho} ( x,y,z; t ) & \mbox{for} \ r> r_S\\[8pt] \displaystyle {\overline{\rho}}^{*} ( x,y,z; t ) = 0 &\mbox{for} \ r < r_S. \end{cases} \end{align} $$

New Fourier coefficients

(3.17) $$ \begin{align} C_{mnp}^{*} (t) & = \frac {1}{L^3 [ 1 + \delta_{m,0} ] [ 1 + \delta_{n ,0} ] [ 1 + \delta_{p,0} ]} \int_{-L}^L \int_{-L}^L \int_{-L}^L {\overline{\rho}}^{*} ( x,y,z; t ) \nonumber \\[6pt] & \quad \times \cos ( \alpha_m (x + L) ) \cos ( \alpha_n (y + L) ) \cos ( \alpha_p (z + L) ) \, dx \, dy \, dz \end{align} $$

are obtained from Fourier analysis of the new perturbation density function (3.16). The previous density coefficients $C_{mnp} (t)$ are now replaced with these new ones $C_{mnp}^{*} (t)$ and the time-integration of the system of differential equations continues. In practice, we sometimes apply a technique such as Lanczos smoothing to these Fourier coefficients; this is equivalent to replacing the discontinuous function (3.16) with a function that varies rapidly but continuously across the region of possible discontinuity.

Initial conditions for this computation usually consist of assuming that the starting velocity is simply the background produced by the distributed source, and this is achieved simply by setting $A_{mnp} (0) = 0$ , $B_{mnp} (0) = 0$ . The initial “shape” of the interfacial region is usually imposed by creating some density perturbation function

(3.18) $$ \begin{align} \begin{cases} \displaystyle \overline{\rho} ( x,y,z,0 ) = D - 1 & \mbox{for} \ r> \mathcal{R}_0\\[8pt] \displaystyle \overline{\rho} ( x,y,z,0 ) = 0 & \mbox{for} \ r < \mathcal{R}_0, \end{cases} \end{align} $$

and Fourier-analysing it, exactly as in (3.17), to obtain initial coefficients $C_{mnp} (0)$ . The function $\mathcal {R}_0$ can be chosen to give the desired starting profile and we have often used

(3.19) $$ \begin{align} \mathcal{R}_0 ( \phi, \theta ) = \epsilon P_n^m ( \cos\phi ) \cos (m\theta ) \end{align} $$

in spherical polar coordinates with elevation angle $\phi = \arctan (\!\sqrt {x^2 + y^2} / z )$ and azimuthal angle $\theta = \arctan ( y / x )$ . The function $P_n^m$ is the associated Legendre function of the first kind, with integer indices m and n (Abramowitz and Stegun [1, page 332]).

For a fully three-dimensional computation such as this, it is generally a nontrivial task to visualize a function such as the density perturbation $\overline {\rho } ( x,y,z; t )$ at some time t, since the function still depends on three spatial variables. Our approach is to draw a surface at the mid-range value $\overline {\rho } = (D - 1) / 2$ and to regard this as the approximate location of the (diffuse) interface. This is done using linear interpolation. On a grid of mesh points $( x_a , y_b , z_c )$ used to evaluate the expressions and to carry out the numerical quadratures in expressions such as (3.10), (3.11), (3.14), we fix the values of $y_b$ and $z_c$ and consider the one-dimensional function

$$ \begin{align*} \overline{\rho}_{\textrm{draw}} (x) = \overline{\rho} ( x , y_b , z_c; t ) - (D - 1) / 2. \end{align*} $$

The aim is to locate the point $x^{*}$ at which $\overline {\rho }_{\textrm {draw}} ( x^{*} ) = 0$ . This is done by considering the product

$$ \begin{align*} \overline{\rho}_{\textrm{draw}} ( x_{a+1} ) \overline{\rho}_{\textrm{draw}} ( x_a ) \end{align*} $$

at neighbouring mesh points. If this is negative, then the perturbation density $\overline {\rho }$ takes its mid-range value between the two points at the approximate location

(3.20) $$ \begin{align} x^{*} = \frac {x_a \overline{\rho}_{\textrm{draw}} ( x_{a+1} ) - x_{a+1} \overline{\rho}_{\textrm{draw}} ( x_a )} {\overline{\rho}_{\textrm{draw}} ( x_{a+1} ) - \overline{\rho}_{\textrm{draw}} ( x_a )}. \end{align} $$

A marker is then drawn at the point $( x^{*} , y_b , z_c )$ , and this test is continued until every neighbouring pair of points in the lattice has been considered. We repeat this process in the other two variables y and z, to ensure that the full surface is represented.

4 Stability analysis for a point source

In this section, we re-cast the inviscid equations in spherical polar coordinates and then carry out a linearized stability analysis of an interface which is now allowed to become an infinitesimally thin surface, across which the density perturbation jumps discontinuously from its inner value $\overline {\rho } = 0$ to the value $\overline {\rho } = D - 1$ outside. For purely axi-symmetric outflow, this has been undertaken in part by Forbes [10, 11] who found the rather remarkable result that, at least for infinite Froude number $1/F = 0$ , the Fourier–Legendre mode of lowest nonspherical order $n = 1$ is always unstable, with the largest growth rate; consequently, any initial shape of the interface can be expected to evolve into a one-sided jet through either one of the two poles of a sphere centred at the origin. Forbes’ numerical results evidently bore out that prediction in the fully nonlinear case too. The purpose of this section is to determine what effect fully three-dimensional geometry has on that prediction.

In the inviscid approximation $1 / R_e = 0$ , the term involving $\nabla ^2 \mathbf {q}$ disappears from the Navier–Stokes equation (2.7). The source region $0 < r < r_S$ is allowed to shrink to a point source on the origin that nevertheless still produces the same outflow flux $f_{\ell }$ as previously. The flow becomes irrotational everywhere in the fluid outside the singular point source, $r> 0$ , and so

(4.1) $$ \begin{align} \mathbf{q}_1 = \nabla \Phi_1 , \quad \mathbf{q}_2 = \nabla \Phi_2 , \end{align} $$

in which $\mathbf {q}_1$ is the velocity vector inside the interface $r < R ( \phi , \theta , t )$ , and the scalar function $\Phi _1 ( r , \phi , \theta , t )$ is its velocity potential. Close to the point source at the origin, the flow becomes singular with

(4.2) $$ \begin{align} \mathbf{q}_1 \rightarrow \frac {f_{\ell} \mathbf{r}}{4 \pi r^3} \quad\textrm{as } r \rightarrow 0. \end{align} $$

The velocity outside the interface $r> R ( \phi , \theta , t )$ is $\mathbf {q}_2$ and it has velocity potential $\Phi _2$ as indicated in (4.1).

On the interface $r = R ( \phi , \theta , t )$ , each fluid obeys its own kinematic boundary condition

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

Here, the velocity components are denoted $( u , w , v )$ in the radial r-direction, the declination angle $\phi $ and the azimuthal $\theta $ -direction, respectively. The dynamic condition across the interface is that the pressures either side must be equal there. Since pressures in each fluid can be evaluated from the irrotational form of the Bernoulli equation, it follows that the dynamical condition is

(4.4) $$ \begin{align} D \frac {\partial \Phi_2}{\partial t} - \frac {\partial \Phi_1}{\partial t} + \frac {1}{2} D \| \mathbf{q}_2 \|^2 - \frac {1}{2} \| \mathbf{q}_1 \|^2 - \frac {( D - 1 )}{F^2} \frac {1}{R} = \frac {(D - 1) f_{\ell}^2}{32 \pi^2 R_0^4 (t)} - \frac {(D - 1)}{F^2 R_0 (t)} \quad\textrm{on } r = R , \end{align} $$

where the spherical radius function

(4.5) $$ \begin{align} R_0 (t) = \bigg[ 1 + \frac {3 f_{\ell} t}{4 \pi} \bigg]^{1/3} \end{align} $$

is the radius of a purely spherical bubble of inner fluid that started with unit radius and grew with time under the effect of the outflow flux $f_{\ell }$ .

Since the two fluids inside and outside the interface are both assumed to be incompressible, then Laplace’s equation $\nabla ^2 \Phi _j = 0$ holds for each velocity potential $j = 1 , 2$ in (4.1), over its domain of validity. It follows therefore that the inviscid model has an exact solution for spherical outflow from a point source at the origin, with purely radial velocity components

$$ \begin{align*} u_1 = u_2 = \frac {f_{\ell}}{4 \pi r^2} , \end{align*} $$

and interface location $r = R_0 (t)$ obtained from (4.5). We now seek to linearize about this base spherical outflow, postulating perturbations to it of the forms

(4.6) $$ \begin{align} \begin{aligned} \Phi_1 ( r , \phi , \theta , t ) & = -f_{\ell} / ( 4 \pi r ) + \varepsilon \Phi^L_1 ( r , \phi , \theta , t ) + \mathcal{O} ( \varepsilon^2 ), \\[.1cm] \Phi_2 ( r , \phi , \theta , t ) & = -f_{\ell} / ( 4 \pi r ) + \varepsilon \Phi^L_2 ( r , \phi , \theta , t ) + \mathcal{O} ( \varepsilon^2 ), \\ R ( \phi, \theta , t ) & = R_0 (t) + \varepsilon R^L ( \phi, \theta , t ) + \mathcal{O} ( \varepsilon^2 ). \end{aligned} \end{align} $$

The dimensionless constant $\varepsilon $ is a small parameter related to the magnitude of the initial perturbation made to the base spherical outflow.

As a result of this linearization process (4.6), the linearized potentials $\Phi ^L_1$ and $\Phi ^L_2$ still satisfy Laplace’s equation, but now in the linearized domains $r < R_0 (t)$ and $r> R_0 (t)$ , respectively. In spherical polar coordinates, we therefore seek solutions

(4.7) $$ \begin{align} \begin{aligned} \Phi^L_1 ( r , \phi , \theta , t ) & = P^L_1 (t) r^n P_n^m ( \cos\phi ) \cos (m\theta ) , & r < R_0 (t), \\[.1cm] \Phi^L_2 ( r , \phi , \theta , t ) & = P^L_2 (t) r^{- (n + 1)} P_n^m ( \cos\phi ) \cos (m\theta ) , & r> R_0 (t), \end{aligned} \end{align} $$

in which the two functions $P^L_1 (t)$ and $P^L_2 (t)$ must be determined. The perturbation $R^L$ to the interface location in (4.6) is sought in a similar form

(4.8) $$ \begin{align} R^L ( \phi, \theta , t ) = A^L (t) P_n^m ( \cos\phi ) \cos (m\theta ) , \end{align} $$

where the amplitude function $A^L (t)$ must likewise be found from the boundary conditions on the interface.

Under the linearization process (4.6), the two kinematic boundary conditions (4.3) become

(4.9) $$ \begin{align} U^L_1 = U^L_2 = \frac {f_{\ell}}{2\pi} \frac {R^L}{R_0^3 (t)} + \frac {\partial R^L}{\partial t} \quad\textrm{on } r = R_0 (t). \end{align} $$

In this expression, the two functions $U^L_1 = \partial \Phi ^L_1 / \partial r$ and $U^L_2 = \partial \Phi ^L_2 / \partial r$ are the linearized perturbations to the radial component of the fluid velocity vector. The dynamic condition (4.4) likewise linearizes to give

(4.10) $$ \begin{align} & D \frac {\partial \Phi^L_2}{\partial t} - \frac {\partial \Phi^L_1}{\partial t} + \frac {f_{\ell}}{4\pi R_0^2 (t)} [ D U^L_2 - U^L_1 ] - (D - 1) \frac {f_{\ell}^2}{8 \pi^2 R_0^5 (t)} R^L \nonumber \\ & \quad+ \frac {( D - 1 )}{F^2} \frac {R^L}{R_0^2 (t)} = 0 \quad\textrm{on } r = R_0 (t). \end{align} $$

The two linearized kinematic conditions (4.9) at once permit the unknown functions $P^L_1 (t)$ and $P^L_2 (t)$ in the expressions (4.7) to be eliminated in favour of the amplitude function $A^L (t)$ in (4.8). We obtain

$$ \begin{align*} n R_0^{n-1} (t) P^L_1 (t) & = \frac {f_{\ell}}{2\pi} \frac {A^L (t)}{R_0^3 (t)} + \frac {d A^L (t)}{dt}, \nonumber \\[.1cm] - (n + 1) R_0^{-n-2} (t) P^L_2 (t) & = \frac {f_{\ell}}{2\pi} \frac {A^L (t)}{R_0^3 (t)} + \frac {d A^L (t)}{dt}. \nonumber \end{align*} $$

These expressions are now incorporated into the linearized dynamic condition (4.10). After some algebra,

(4.11) $$ \begin{align} &( D n + n + 1 ) \bigg[ R_0^3 (t) \frac {d^2 A^L (t)}{dt^2} + \frac {3 f_{\ell}}{4\pi} \frac {d A^L (t)}{dt} \bigg] \nonumber \\ & \quad + \bigg\{ 2\bigg( \frac {f_{\ell}}{4\pi}\bigg)^2 R_0^{-3} (t) [ Dn(n-1) - (n+1) (n+2) ] - \frac {( D - 1 )}{F^2} n (n+1) \bigg\} A^L (t) = 0. \end{align} $$

This is a difficult ordinary differential equation of second order for the amplitude function $A^L (t)$ that describes the growth or decay of the linearized perturbation (4.8) to the basic spherical outflow.

Although it is linear, the differential equation (4.11) has coefficients that are functions of time, and this makes it difficult to treat. This is addressed by making the spherical radius function $R_0 (t)$ the independent variable, in place of time t. Using the chain rule of calculus and the definition (4.5) then converts (4.11) to the more manageable form

(4.12) $$ \begin{align} R_0^2 \frac {d^2 A^L}{d R_0^2} + R_0 \frac {d A^L}{d R_0} + ( a_n - b_n R_0^3 ) A^L = 0 , \end{align} $$

in which the two additional constants have been defined to be

(4.13) $$ \begin{align} \begin{aligned} a_n & = \frac {2 [ Dn (n-1) - (n+1) (n+2) ]}{( Dn + n + 1 )}, \\[.1cm] b_n & = \frac {( D - 1 )}{F^2} \bigg( \frac {4\pi}{f_{\ell}} \bigg)^2 \frac {n (n+1)}{( Dn + n + 1 )}. \end{aligned} \end{align} $$

Several cases of interest now present themselves based on this equation.

4.1 No gravitation

If gravitational effects are entirely absent, then the flow is no longer a genuine RT outflow; instead, convergence effects due to the spherical geometry are dominant. This is known as the Bell–Plesset effect and is a feature of situations in which the resting interface between the fluids possesses finite curvature. A critical analysis of Bell–Plesset flows is given by Epstein [8]. In the present inviscid analysis, in which the two fluids are separated by an interface of infinitesimal thickness, vorticity deposition at the interface occurs as time progresses (see Peng et al. [22]), so that the interface becomes a vortex sheet in the fluid.

When there is no gravity, then $1 / F^2 = 0$ . In that case, $b_n = 0$ in (4.13) and the differential equation (4.12) becomes simply

$$ \begin{align*} R_0^2 \frac {d^2 A^L}{d R_0^2} + R_0 \frac {d A^L}{d R_0} + a_n A^L = 0. \end{align*} $$

This is an Euler–Cauchy equation and its general solution takes the form

(4.14) $$ \begin{align} A^L ( R_0 ) = C_1 R_0^{\sqrt{- a_n}} + C_2 R_0^{- \sqrt{- a_n}}, \end{align} $$

in which $C_1$ and $C_2$ are arbitrary constants and $a_n$ is as defined in (4.13).

Whenever

$$ \begin{align*} D> \frac {(n+1) (n+2)}{n (n - 1)}, \end{align*} $$

then $a_n> 0$ and (4.14) can be written in the equivalent real form

$$ \begin{align*} A^L ( R_0 ) = D_1 \cos ( \sqrt{a_n} \log R_0 ) + D_2 \sin ( \sqrt{a_n} \log R_0 ) \end{align*} $$

with $D_1$ and $D_2$ arbitrary real constants. In this case, the solution for the amplitude function $A^L$ is clearly bounded so that the spherical outflow is stable.

However, when

(4.15) $$ \begin{align} D < \frac {(n+1) (n+2)}{n (n - 1)} \end{align} $$

with positive outflow $f_{\ell }> 0$ , then $a_n < 0$ and $R_0$ is a monotonically increasing function of time t. In this case, the solution (4.14) shows that the $(m,n)$ solution mode in (4.8) will be unstable for all azimuthal modes m and perturbations to the spherical outflow will grow as t increases. The condition (4.15) for instability was also obtained by Forbes [10]. He observed, in particular, that the instability criterion (4.15) would always be satisfied by the first mode $n = 1$ , so that these solutions are always unstable regardless of the density ratio D. Forbes [10] concluded, therefore, that the first mode $n = 1$ would always come to dominate the outflow morphology, so that eventually, one-sided outflows would be observed.

4.2 The $\boldsymbol {n = 1}$ modes

In the special case of the first modes $n = 1$ (for which the two possibilities $m = 0$ and $m = 1$ exist), the two constants in (4.13) become

(4.16) $$ \begin{align} \begin{aligned} a_1 & = - \alpha_1^2 \quad\textrm{with } \alpha_1 = \sqrt{ \frac{12}{D + 2}}, \\[.1cm] b_1 & = \frac {( D - 1 )}{F^2} \bigg( \frac {4\pi}{f_{\ell}} \bigg)^2 \frac {2}{( D + 2 )}. \end{aligned} \end{align} $$

We make the change of independent variable $\xi = R_0^{3/2}$ in the differential equation (4.12) and obtain

$$ \begin{align*} \xi^2 \frac {d^2 A^L}{d \xi^2} + \xi \frac {d A^L}{d \xi} - \bigg( \frac {4}{9} \alpha_1^2 + \frac {4}{9} b_1 \xi^2 \bigg) A^L = 0. \end{align*} $$

This is a modified Bessel equation and it has the general solution

(4.17) $$ \begin{align} A^L ( R_0 ) = C_1 I_{\nu} \big( \tfrac {2}{3} \sqrt{ b_1 R_0^3} \big) + C_2 K_{\nu} \big( \tfrac {2}{3} \sqrt{ b_1 R_0^3} \big), \end{align} $$

in which $I_{\nu }$ and $K_{\nu }$ are modified Bessel functions of the first and second kinds, respectively, of order

$$ \begin{align*} \nu = \frac {2}{3} \alpha_1 = \frac {4}{\sqrt{3 ( D + 2 )}}. \end{align*} $$

The two constants $C_1$ and $C_2$ are again arbitrary, and $\alpha _1$ and $b_1$ are as defined in (4.16).

For $D> 1$ , the constant $b_1$ in (4.16) is positive, and so the amplitude $A^L$ in (4.17) increases exponentially with $R_0$ . Using the definition (4.5) and the asymptotic form for $I_{\nu }$ in Abramowitz and Stegun [1, page 377] shows that

$$ \begin{align*} A^L (t) \sim t^{- 1/4} \exp ( \gamma_1 \sqrt{t} ) \quad\textrm{with } \gamma_1 = \frac {2}{3} \sqrt{\frac {3 f_{\ell}}{4\pi} b_1} \quad\textrm{as } t \rightarrow \infty. \end{align*} $$

Thus the first Legendre mode $n = 1$ grows exponentially for any $D> 1$ and Forbes [10] found numerically that this $n = 1$ mode ultimately dominated the outflow shape in the axi-symmetric case $m=0$ .

4.3 Higher modes $\boldsymbol {n> 1}$

In the general case $n> 1$ , similar analysis to that in Section 4.2 shows that the general solution for the amplitude of the nth linearized mode can be written as

(4.18) $$ \begin{align} A^L ( R_0 ) = C_1 I_{\nu} \big( \tfrac {2}{3} \sqrt{ b_n R_0^3} \big) + C_2 K_{\nu} \big( \tfrac {2}{3} \sqrt{ b_n R_0^3} \big) , \end{align} $$

which is of the same form as (4.17), except that now

$$ \begin{align*} \nu = \tfrac {2}{3} \sqrt{- a_n} \end{align*} $$

with $a_n$ defined in (4.13), and this becomes imaginary for sufficiently large n. Furthermore, if the inner fluid is more dense, $0 < D < 1$ , then $\sqrt {b_n}$ becomes imaginary and the modified Bessel functions $I_{\nu }$ , $K_{\nu }$ in the solution (4.18) may be replaced with ordinary Bessel functions $J_{\nu }$ , $Y_{\nu }$ of a real argument. Thus higher modes are stable for $D < 1$ .

5 Stability analysis for a finite source

In the previous section, an inviscid stability analysis was carried out, under the assumption that the source could be regarded as a mathematical point singularity. However, as the focus of this present article is on distributed sources, it is appropriate here to investigate the effect of source radius $r_S$ on stability.

We proceed as in Section 4, again considering two inviscid fluids separated by a sharp interface $r = R ( \phi , \theta , t )$ . There are again two velocity potentials $\Phi _1$ and $\Phi _2$ each satisfying Laplace’s equation in their respective fluid domains, and the fluid velocities $\mathbf {q}_1$ and $\mathbf {q}_2$ are calculated from their gradients, as in (4.1). Each fluid obeys a kinematic boundary condition (4.3) on its interface and a dynamic condition of the form (4.4) holds also on the interface. However, since a source of finite radius $r_S$ is now of interest, the singular boundary condition (4.2) for a point source at the origin is replaced with the finite constraint

(5.1) $$ \begin{align} \mathbf{q}_1 = \frac {f_{\ell}}{4 \pi r_S^2} \mathbf{e}_r \quad\textrm{on } r = r_S , \end{align} $$

in which $\mathbf {e}_r$ is the spherical unit vector directed radially outward from the origin.

In view of the choice (5.1) of boundary condition at the finite spherical source, the same base spherical outflow used in Section 4 applies here also, with spherical radius function (4.5) as previously. Accordingly, the same linearization procedure (4.6) is undertaken, and leads to the same linearized kinematic and dynamic boundary conditions (4.9), (4.10) to be applied on the moving sphere $r = R_0 (t)$ . However, because of the source boundary condition (5.1) here, the new form for the inner velocity potential $\Phi ^L_1$ , replacing (4.7), is

$$ \begin{align*} \Phi^L_1 ( r , \phi , \theta , t ) = P^L_1 (t) \bigg[ r^n + \frac {n r_S^{2n+1}}{( n+1) r^{n+1}} \bigg] P_n^m ( \cos\phi ) \cos (m\theta ) , \quad r < R_0 (t). \end{align*} $$

The outer potential $\Phi ^L_2$ and the linearized interface perturbation $R^L$ retain the same forms as in (4.7) and (4.8).

After some considerable algebra, a new differential equation for the amplitude function $A^L (t)$ in the linearized perturbation $R^L$ in (4.8) may be derived to replace (4.11). This equation has the form

$$ \begin{align*} R_0 (t) &\bigg[ D n + \frac {\mathcal{L}_0 (t)}{n \mathcal{D}_0 (t)} \bigg] \frac {d^2 A^L (t)}{dt^2} + \frac {f_{\ell}}{4\pi R_0^2 (t)} \bigg[ 3 D + n + 1 - \frac {\mathcal{N}_0 (t) \mathcal{L}_0 (t)}{n \mathcal{D}_0^2 (t)} \bigg] \frac {d A^L (t)}{dt} \\[4pt] & + \bigg[ \frac {f_{\ell}^2}{8 \pi^2 R_0^5 (t)} \bigg\{ D (n - 1) - \frac {\mathcal{M}_0 (t) \mathcal{L}_0 (t)}{n \mathcal{D}_0^2 (t)} \bigg\} - \frac {( D - 1 )}{F^2} \frac {(n+1)}{R_0^2 (t)} \bigg] A^L (t) = 0 , \end{align*} $$

in which it is convenient to define the four intermediate functions

$$ \begin{align*} \begin{aligned} \mathcal{L}_0 (t) & = (n+1) R_0^{2n + 1} (t) + n r_S^{2n + 1} , \\[3pt] \mathcal{M}_0 (t) & = (n+2) R_0^{2n + 1} (t) + (n-1) r_S^{2n + 1} , \\[3pt] \mathcal{N}_0 (t) & = (n-3) R_0^{2n + 1} (t) + (n+4) r_S^{2n + 1} , \\[3pt] \mathcal{D}_0 (t) & = R_0^{2n + 1} (t) - r_S^{2n + 1}. \end{aligned} \end{align*} $$

This is a linear second-order differential equation that determines the stability of the outflow and it reduces to (4.11) in the limit $r_S \rightarrow 0$ . However, it is extremely difficult to analyse in closed form and even the changes of variable considered in Section 4 appear to offer no assistance.

To make progress, we have therefore considered an approximate “quasi-stationary” type of analysis, in which it is assumed that the spherical base radius $R_0 (t)$ in (4.5) is approximately constant. A possible justification for this assumption might come from the fact that $R_0$ behaves as $t^{1/3}$ for large t, and so might be considered to operate on a longer time scale than that over which any instability might develop. To simplify the analysis as much as possible, we have considered the early stages of the flow and made the (crude) approximation $R_0 \approx 1$ . This then yields the approximate stability equation

(5.2) $$ \begin{align} \mathcal{A}_n \frac {d^2 A^L (t)}{dt^2} + \mathcal{B}_n \frac {d A^L (t)}{dt} + \mathcal{C}_n A^L (t) = 0 , \end{align} $$

where we have defined the three additional constants

(5.3) $$ \begin{align} \begin{aligned} \mathcal{A}_n & = [ 1 - r_S^{2n + 1} ] [ (Dn + n + 1) - (D - 1) n r_S^{2n + 1} ], \\[3pt] \mathcal{B}_n & = \frac {f_{\ell}}{4\pi} [ 3(Dn + n + 1) -2 ( 3Dn + 2n^2 + 2n + 2 ) r_S^{2n + 1} + 3(D-1) n r_S^{4n + 2} ], \\[3pt] \mathcal{C}_n & = - \frac {(D - 1)}{F^2} n (n+1) [ 1 - r_S^{2n + 1} ]^2 + \frac {f_{\ell}^2}{8\pi^2}\{ D n (n-1) [ 1 - r_S^{2n + 1} ]^2 \\[3pt] & \quad - [ (n+2) + (n-1) r_S^{2n + 1} ] [ (n+1) + n r_S^{2n + 1} ] \}. \end{aligned} \end{align} $$

Since the approximate equation (5.2) involves only the constant coefficients defined in (5.3), it has simple exponential solutions behaving like $\exp (\lambda t)$ in which the exponent $\lambda $ satisfies a quadratic equation. Furthermore, because $r_S < 1$ , it follows that $\mathcal {A}_n$ is always positive and so the growth-rate constant $\lambda $ can be guaranteed to be positive if $- \mathcal {C}_n> 0$ .

For the first mode $n=1$ , we observe that

$$ \begin{align*} - \mathcal{C}_1 = \frac {2 (D - 1)}{F^2} [ 1 - r_S^3 ]^2 + \frac {3 f_{\ell}^2}{8\pi^2} [ 2 + r_S^3 ] \end{align*} $$

and since this is positive, this approximate analysis therefore predicts that the first mode is always unstable for any permissible source radius $r_S$ .

Two different stability situations are illustrated in Figure 1, according to this simple quasi-stationary analysis. In Figure 1(a), we display stability curves for parameter values that reflect the situation discussed by Forbes [10], in which there is essentially no gravity ( $1/F^2 = 10^{-10}$ ) and the density ratio $D = 6$ is large. The quantity $- \mathcal {C}_n$ given in (5.3) is plotted against source radius $r_S$ and the results in Figure 1(a) concur with the more precise analysis in Section 4.1. In particular, the instability criterion (4.15) indicates that while the $n=1$ mode is always unstable, the second Fourier mode $n=2$ changes from unstable to stable as the density ratio D passes through the value $D=6$ . This can be seen in Figure 1(a), where $- \mathcal {C}_2 \approx 0$ for the second mode $n=2$ at smaller source radii $r_S < 0.6$ (although for larger values of $r_S$ , this mode too evidently becomes unstable). For smaller source radii, all the higher modes $n \geq 3$ have $- \mathcal {C}_2 < 0$ and so are expected to be stable. Therefore, the analysis of this section suggests that a one-sided jet would be formed for this case $D=6$ when the source radius $r_S$ is sufficiently small, as found by Forbes [10].

Figure 1 Stability curves for the first five spherical modes $n = 1 , \ldots , 5$ . The constant $- \mathcal {C}_n$ obtained from (5.3) is plotted against the source radius $r_S$ for (a) zero-gravity source with $D = 6$ and (b) density ratio $D = 1.05$ and Froude number $F = 1$ .

An interesting feature is present in Figure 1(a) in a relatively narrow band of values of source radius $r_S$ near one. Although the first mode $n=1$ is the only unstable one for smaller radius $r_S$ , the higher modes become dominant in the approximate interval $0.9 < r_S <1$ . As a result, at $r_S = 1$ , the growth rate increases with the mode order n. Therefore, at low source radius $r_S$ , the eventual outflow shape is predicted to be a one-sided configuration, but for $r_S \approx 1$ , the outflow morphology would be expected to reflect the initial interface shape. Although the results in Figure 1 correspond to purely inviscid outflows, a somewhat similar phenomenon was observed by Terrones and Heberling [26] in the viscosity-dominated case. They found that a similar stability exchange occurs as a parameter corresponding to a type of Reynolds number was made to increase, so that the outflow would be dominated by progressively higher Fourier modes as this occurred.

Figure 1(b) shows stability curves for Froude number $F = 1$ and density ratio $D = 1.05$ . In this case, all the spherical Fourier modes are predicted to be unstable and the degree of instability increases with the mode number n. This is true regardless of the source radius $r_S$ and for this case, the shape of the outflow as time develops would most likely reflect the original interface shape at $t = 0$ .

6 Complex distributed source

So far, we have considered only simple spherical outflow, either as a singularity (2.1), similar to the previous studies of Forbes [10, 11], or as a source (5.1) of finite radius $r_S$ . However, a distributed source over the finite region $r < r_S$ provides an opportunity to consider more complex outflow patterns and these are addressed briefly here.

In addition to the mass-producing radial outflow (2.1), we have also introduced higher spherical modes into the source behaviour to study their effect on outflow morphologies. In spherical polar coordinates $( r , \phi , \theta )$ , the fluid velocity vector $\mathbf {q}_S$ generated by the source can be expressed as

$$ \begin{align*} \mathbf{q}_S = u_S^{\textrm{sph}} \mathbf{e}_r + w_S^{\textrm{sph}} \mathbf{e}_{\phi} + v_S^{\textrm{sph}} \mathbf{e}_{\theta}, \end{align*} $$

in which the unit vector $\mathbf {e}_r = \mathbf {r} / r$ points radially outward, as in (2.1). The remaining two unit vectors $\mathbf {e}_{\phi }$ and $\mathbf {e}_{\theta }$ point in the direction of the positive declination angle $\phi $ down from the z-axis and positively in the azimuthal $\theta $ -direction, respectively. The simple radial outflow (2.4) is now generalized to become

(6.1) $$ \begin{align} u_S^{\textrm{sph}} = \begin{cases} \displaystyle \frac{f_{\ell}}{4\pi r^2} - \varepsilon_n \frac {(n+1)}{r^{n+2}} P_n^{\ell} ( \cos\phi ) \cos ( \ell\theta ) & \mbox{for} \ r> R_S \\[8pt] \displaystyle \frac {f_{\ell} r^2}{4\pi r_S^4} - \varepsilon_n \frac {(n+1) r^2}{r_S^{n+4}} P_n^{\ell} ( \cos\phi ) \cos ( \ell\theta ) & \mbox{for} \ r < R_S \end{cases} \end{align} $$

in the radial direction,

$$ \begin{align*} w_S^{\textrm{sph}} = \begin{cases} \displaystyle - \frac {\varepsilon_n}{r^{n+2}} ( P_n^{\ell} )^{\prime} ( \cos\phi ) \, \sin\phi \cos ( \ell\theta ) \quad \mbox{for} \ r> R_S\\[8pt] \displaystyle - \frac {\varepsilon_n r^2}{r_S^{n+4}} ( P_n^{\ell} )^{\prime} ( \cos\phi ) \, \sin\phi \cos ( \ell\theta ) \quad \mbox{for} \ r < R_S \end{cases} \end{align*} $$

in the direction of the declination angle $\phi $ , and

(6.2) $$ \begin{align} v_S^{\textrm{sph}} = \begin{cases} \displaystyle - \varepsilon_n \frac {\ell}{r^{n+2} \sin\phi} P_n^{\ell} ( \cos\phi ) \sin ( \ell\theta ) \quad \mbox{for} \ r> R_S\\[8pt] \displaystyle - \varepsilon_n \frac {\ell r^2}{r_S^{n+4} \sin\phi} P_n^{\ell} ( \cos\phi ) \sin ( \ell\theta ) \quad \mbox{for} \ r < R_S \end{cases} \end{align} $$

azimuthally. The Cartesian velocity components $( u_S , v_S , w_S )$ generated by this complex source (6.1)–(6.2) are then found from

(6.3) $$ \begin{align} \begin{aligned} u_S & = u_S^{\textrm{sph}} \sin\phi \cos\theta + w_S^{\textrm{sph}} \cos\phi \cos\theta - v_S^{\textrm{sph}} \sin\theta , \\ v_S & = u_S^{\textrm{sph}} \sin\phi \sin\theta + w_S^{\textrm{sph}} \cos\phi \sin\theta + v_S^{\textrm{sph}} \cos\theta , \\ w_S & = u_S^{\textrm{sph}} \cos\phi - w_S^{\textrm{sph}} \sin\phi. \end{aligned} \end{align} $$

The divergence of the velocity field can be calculated directly from its polar form (6.1)–(6.2) using Batchelor [5, page 601]. We obtain

(6.4) $$ \begin{align} \begin{cases} \displaystyle \textrm{div} \mathbf{q}_S = 0 \quad \mbox{for} \ r> r_S \\[8pt] \displaystyle \textrm{div} \mathbf{q}_S = \frac {f_{\ell} r}{\pi r_S^4} - \varepsilon_n \frac {(n+1)(n+4) r}{r_S^{n+4}} P_n^{\ell} ( \cos\phi ) \cos ( \ell\theta ) \quad \mbox{for} \ r < r_S. \end{cases} \end{align} $$

The divergence is zero outside the source region, $r> r_S$ , as required, since no mass is generated there. Furthermore, when $\textrm {div} \mathbf {q}_S$ in (6.4) is integrated over the source volume $r < r_S$ , the total flux of the source is calculated to be $f_{\ell }$ as required. Thus the extra higher-order terms in the complex source produce no net additional flux.

In the results to be discussed in Section 7, we will focus on the bi-polar case $n=2$ . Indeed, this seems likely to be of the most interest in astrophysical applications, due to the additional presence of a magnetic field. In the purely axi-symmetric case $P_2^0 ( \cos \phi ) = (1/2) ( 3\cos ^2 \phi - 1 )$ , the three equations in (6.3), combined with the spherical forms (6.1)–(6.2), show that the appropriate Cartesian velocity components in this case are

(6.5) $$ \begin{align} \begin{cases} \displaystyle u_S = \frac {x}{r} \bigg[ \frac {f_{\ell}}{4\pi r^2} - \frac {3\varepsilon_2}{2 r^4} \bigg( 5 \frac {z^2}{r^2} - 1 \bigg) \bigg] \quad \mbox{for} \ r> r_S\\[.4cm] \displaystyle u_S = \frac {x}{r} \bigg[ \frac {f_{\ell} r^2}{4\pi r_S^4} - \frac {3\varepsilon_2}{2 r_S^6} ( 5 z^2 - r^2 ) \bigg] \quad \mbox{for} \ r < r_S , \end{cases} \end{align} $$

for the component along the x-axis. The y-component of the source velocity is the same, except that the variable x is replaced with y. The component on the z-axis is

(6.6) $$ \begin{align} \begin{cases} \displaystyle w_S = \frac {z}{r} \bigg[ \frac {f_{\ell}}{4\pi r^2} - \frac {3\varepsilon_2}{2 r^4} \bigg( 5 \frac {z^2}{r^2} - 3 \bigg) \bigg] & \mbox{for} \ r> r_S\\[.4cm] \displaystyle w_S = \frac {z}{r} \bigg[ \frac {f_{\ell} r^2}{4\pi r_S^4} - \frac {3\varepsilon_2}{2 r_S^6} ( 5 z^2 - 3 r^2 ) \bigg] & \mbox{for} \ r < r_S. \end{cases} \end{align} $$

By contrast, the most severely nonaxi-symmetric bi-polar outflow is obtained with $P_2^2 ( \cos \phi ) = 3 \sin ^2 \phi $ . For this case, the Cartesian velocity component directed along the x-axis is

(6.7) $$ \begin{align} \begin{cases} \displaystyle u_S = \frac {x}{r} \bigg[ \frac {f_{\ell}}{4\pi r^2} + \frac {6\varepsilon_2}{r^4} + 15 \varepsilon_2 \frac {( y^2 - x^2 )}{r^6} \bigg] & \mbox{for} \ r> r_S\\[.4cm] \displaystyle u_S = \frac {x}{r} \bigg[ \frac {f_{\ell} r^2}{4\pi r_S^4} + 6\varepsilon_2 \frac {r^2}{r_S^6} + 15 \varepsilon_2 \frac {( y^2 - x^2 )}{r_S^6} \bigg] & \mbox{for} \ r < r_S. \end{cases} \end{align} $$

The y-directed velocity component is similar, except that the $x/r$ term to the left of the large parentheses is replaced by $y/r$ . In addition, the sign of the second term is changed from positive to negative. The z-component of the source velocity takes the slightly simpler form

(6.8) $$ \begin{align} \begin{cases} \displaystyle w_S = \frac {z}{r} \bigg[ \frac {f_{\ell}}{4\pi r^2} + 15 \varepsilon_2 \frac {( y^2 - x^2 )}{r^6} \bigg] \quad \mbox{for} \ r> r_S\\[.4cm] \displaystyle w_S = \frac {z}{r} \bigg[ \frac {f_{\ell} r^2}{4\pi r_S^4} + 15 \varepsilon_2 \frac {( y^2 - x^2 )}{r_S^6} \bigg] \quad \mbox{for} \ r < r_S. \end{cases} \end{align} $$

For completeness, we also give the divergence

(6.9) $$ \begin{align} \begin{cases} \displaystyle \textrm{div} \mathbf{q}_S = 0 \quad \mbox{for} \ r> r_S\\[8pt] \displaystyle \textrm{div} \mathbf{q}_S = \frac {f_{\ell} r}{\pi r_S^4} + 54 \varepsilon_2 \frac {r}{r_S^6} \frac {( y^2 - x^2 )}{r^2} \quad \mbox{for} \ r < r_S , \end{cases} \end{align} $$

since this quantity is needed in the nonlinear terms in (3.9).

It is not possible to carry out a full stability analysis for this more complex flow, as was possible in Section 4 for pure radial source flow, because the unperturbed interface in the inviscid case is no longer the simple spherical shape (4.5). However, if the amplitude $\varepsilon _n$ is regarded as a small parameter, then it is again possible to linearize about the purely radial source flow with its spherical interface, precisely as was done in (4.6). The inviscid potentials are then sought in the forms

$$ \begin{align*} \Phi^L_1 ( r , \phi , \theta , t ) & = \bigg[ P^L_1 (t) r^n + \frac {1}{r^{n+1}} \bigg] P_n^m ( \cos\phi )\cos (m\theta ), \quad r < R_0 (t), \\[4pt] \Phi^L_2 ( r , \phi , \theta , t ) & = P^L_2 (t) r^{- (n + 1)} P_n^m ( \cos\phi ) \cos (m\theta ) , \quad r> R_0 (t) , \end{align*} $$

generalizing the previous expression (4.7), and the nonspherical perturbation to the interface is

$$ \begin{align*} R^L ( \phi, \theta , t ) = A^L (t) P_n^m ( \cos\phi ) \cos (m\theta ) , \end{align*} $$

exactly as in (4.8). The analysis proceeds as in Section 4 and after some algebra, the equation for the linearized interface amplitude function $A^L (t)$ , replacing (4.12), is

$$ \begin{align*} R_0^2 \frac {d^2 A^L}{d R_0^2} + R_0 \frac {d A^L}{d R_0} + ( a_n - b_n R_0^3 ) A^L = c_n R_0^{1-n}. \end{align*} $$

The two constants $a_n$ and $b_n$ are as defined in (4.13) and the additional constant on the right-hand side is

$$ \begin{align*} c_n = \bigg( \frac {4\pi}{f_{\ell}} \bigg) \frac {(n+1)^2 (2n + 1)}{Dn + n + 1}. \end{align*} $$

This is just an inhomogeneous version of the previous amplitude equation (4.12), which is to be expected since it has been necessary to treat the nonspherical contributions to the flow as small perturbations to the pure spherical outflow, exactly as in Section 4, and so the same stability conditions are obtained here also.

7 Results

The spectral technique described in Section 3, and the modified form of it in Section 6, has been coded in a parallelized form of the Fortran programming language and run on a purpose-built high-performance computer using graphics-card based architecture. Specifically, the computations were performed using 64-bit floats on an Nvidia GP100 graphics card using Nvidia’s NVFortran compiler, on Ubuntu 18.04. The classic fourth-order Runge–Kutta method was used to march the solution forward in time (Atkinson [4, page 371]). At regular intervals, points on the “interface” were identified by (3.20) and saved to file. These points were then plotted using Matlab’s “pcshow” function. Since these numerical results are based on the Boussinesq approximation, we have taken the density ratio to be $D = 1.05$ consistent with that theory. The inverse Reynolds number $1/R_e$ and the density diffusion parameter $\sigma $ have both been given values in the interval $10^{-4}$ to $10^{-3}$ to control small-wavelength instabilities in the numerical results and the overall outflow morphologies are not affected by these choices.

7.1 Spherical distributed source

We begin this presentation of results by investigating perturbations to the spherical distributed source in Sections 2 and 3. Figure 2 shows the interfacial zone in which the initial velocity profile (3.4) had been perturbed at the first Fourier–Legendre mode $n = 1$ . The disturbance was axisymmetric, with $m = 0$ , and had initial amplitude $\epsilon = 0.025$ . The density ratio is $D = 1.05$ , with Froude number $F = 1$ , inverse Reynolds number $1 / R_e = 4 \times 10^{-3}$ and diffusion constant taken to be $\sigma = 10^{-4}$ . The radius of the distributed source is $r_S = 0.6$ and its total integrated flux is $f_{\ell } = 1$ . The outflow is displayed at the four different (dimensionless) times $t = 4$ , $8$ , $12$ and $16$ , and these solutions were obtained using $M = N = P = 24$ Fourier modes in each of the three physical coordinates, giving a total of $13,824$ modes overall and $120$ grid points in each coordinate variable. The size of the computational region is described by the parameter $L = 3$ .

Figure 2 Outflow from a spherical distributed source ( $f_{\ell } = 1$ ), perturbed at the first axisymmetric Legendre mode $n=1$ , $m=0$ with amplitude $\epsilon = 0.025$ . The interface is shown (left to right) at the four times $t = 4$ , $8$ , $12$ , $16$ .

The initial $n=1$ mode disturbance has maintained its original shape, to a large extent, as time progresses, although nonlinear effects have clearly become more important at later times. Thus, by time $t = 16$ in the last frame on the right of Figure 2, the portion of the outflow at the top of the diagram, for $z> 0$ , has formed an overturning mushroom-shaped cap. In addition, a weaker secondary outflow has formed at the bottom of the diagram, for $z < 0$ , and this also begins to form a small overturning region.

In contrast to Figure 2, we show in Figure 3 a solution with no outflow from the source region, $\,f_{\ell } = 0$ . The source radius in this instance is $r_S = 0.4$ and there is a weaker gravity field for which $1/F^2 = 0.4$ . To obtain these images at sufficient accuracy over this time interval, the number of Fourier modes was increased to $M = N = P = 64$ , with $320$ grid points in each coordinate, and $L = 3$ as the computational boundary.

Figure 3 Evolution of a disturbed spherical distributed region, with no source present ( $f_{\ell } = 0)$ , perturbed at the first axisymmetric Legendre mode $n=1$ , $m=0$ with amplitude $\epsilon = 0.1$ . The interface is shown (left to right) at the four times $t = 10$ , $20$ , $30$ , $40$ .

Since there is no outflow, a source is not present in the results displayed in Figure 3, and so this sequence of diagrams represents the evolution, under a central gravity field, of a bubble with velocity perturbed at the first mode with amplitude $\epsilon = 0.1$ . In the first diagram for $t = 10$ at the left of this figure, the interfacial zone is almost spherical, but at the later times shown, an overturning cap develops and moves outward up the positive z-axis, while continuing to overturn. This is qualitatively similar to the interface shape in Figure 2, except that there is no flow down the negative z-axis, since a source is not present. A remnant of the initial spherical ball centred at the origin persists at these later times, even as the cap moves further away. From the computational viewpoint, it is infeasible to continue to much later times, but it is likely that the spherical cap will eventually detach from the portion remaining near the origin, similar to the “lazy plume” results obtained by Allwright et al. [2].

The source is reintroduced at the origin in Figure 4, so that $f_{\ell } = 1$ and the choice of parameters is the same as in Figure 2. However, the initial condition in this case has been altered to be an axisymmetric velocity perturbation at the second spherical mode, with $n=2$ and $m=0$ and initial amplitude $\epsilon = 0.025$ . In this case, the outflow is dominated by its initial bipolar velocity distribution, which develops and becomes more pronounced as time progresses. By the last time $t = 16$ shown, both outflow jets have overturned to form axisymmetric mushroom-shaped plumes.

Figure 4 Outflow from a spherical distributed source ( $f_{\ell } = 1$ ), with initial velocity perturbed at the second axisymmetric Legendre mode $n=2$ , $m=0$ with amplitude $\epsilon = 0.025$ . The interface is shown (left to right) at the four times $t = 4$ , $8$ , $12$ , $16$ .

Figure 5 again shows an outflow in which there is no source present, $f_{\ell } = 0$ , so that the morphology that develops is entirely due to the initial disturbance, which in this case is a perturbation at the second axisymmetric mode $n = 2$ , $m = 0$ . The other parameters in the solution are identical to those used in Figure 3. As time progresses, the initially spherical interface forms bipolar jets that move in opposite directions along the z-axis, eventually forming overhanging plumes. As with Figure 3, the situation depicted in Figure 5 is analogous to the “lazy plume” computed by Allwright et al. [2], and it is possible that the two caps may eventually detach from the portion that remains near the origin and move separately in opposite directions on the z-axis.

Figure 5 Evolution of a disturbed spherical distributed region, with no source present ( $f_{\ell } = 0)$ , perturbed at the second axisymmetric Legendre mode $n=2$ , $m=0$ with amplitude $\epsilon = 0.1$ . The interface is shown (left to right) at the four times $t = 10$ , $20$ , $30$ , $40$ .

An axisymmetric perturbation at the third Fourier mode $n=3$ , with initial amplitude $\epsilon = 0.025$ , is explored in Figure 6, with nett outflow $f_{\ell } = 1$ from the source region $r < r_S = 0.6$ . This result was obtained using the same parameter values as in Figures 2 and 4. In this case where $n=3$ , $m=0$ , the initial disturbance involves a small flow on the positive z-axis and the negative z-axis as well as a circular outflow on the $xy$ -plane. As time progresses, each of these features is maintained in the evolving interface shape, as is evident in Figure 6. The circular waistband around the $xy$ -plane continues to expand outwards, and there are jets formed up both the positive and negative z-axes. Nevertheless, it is interesting to observe that there is a strong asymmetry in these two jets, with a smaller one on the negative z-axis but a much larger jet up the positive z-axis, that also overturns to form a mushroom-shaped structure. Although the solution was started just with a perturbation to the third mode $n=3$ , nonlinearity is responsible for coupling of all the Fourier–Legendre modes and this evidently allows the $n=1$ mode to grow rapidly with time t, thus contributing to the strong asymmetry about the $xy$ -plane.

Figure 6 Outflow from a spherical distributed source ( $f_{\ell } = 1)$ , perturbed at the third axisymmetric Legendre mode $n=3$ , $m=0$ with amplitude $\epsilon = 0.025$ . The interface is shown (left to right) at the four times $t = 4$ , $8$ , $12$ , $16$ .

7.2 Outflow from a complex source region

We now consider the complex source in Section 6. For these solutions, we choose not to perturb the initial flow field and so $\epsilon = 0$ . Rather, the interface shape that evolves with the passage of time will be driven by the additional higher-order terms in (6.1)–(6.2) with their forcing amplitude $\varepsilon _n$ , and with particular focus on the bipolar modes in (6.5)–(6.9) with forcing amplitudes $\varepsilon _2$ . For the results to follow, we have chosen Reynolds number $R_e = 4 \times 10^{-3}$ , diffusion constant $\sigma = 10^{-4}$ and Froude number $F = 1$ . The density ratio is $D = 1.05$ and the outflow flux is $f_{\ell } = 1$ . In each of these solutions, we have used $M = N = P = 64$ Fourier modes, a total of $262,144$ modes overall, with $320$ mesh points over each of the three coordinate variables. The computational window was defined by the constant $L = 2.5$ in each case.

Figure 7 shows the evolution of the interface obtained with the complex source in (6.5), (6.6) with forcing at the axisymmetric $(n,m) = (2,0)$ mode. In this case, the parameter $\varepsilon _2$ is negative, with value $\varepsilon _2 = - 0.06$ , and solutions are shown at the four times $t = 2$ , $4$ , $6$ , $8$ . The initial interface shape is spherical, but as time progresses, it develops the dumbbell-like shape shown with two lobes centred along the z-axis. These continue to grow, since there is nett outflow $f_{\ell } = 1$ from the source region $r < r_S = 0.6$ , and the two lobes move further apart. We have not been able to follow the development of this solution for times much greater than approximately $t = 8$ , which is the last image of the interface on the right of this figure, primarily due to the accumulation of round-off errors in the enormous volume of computer arithmetic undertaken by this code, as has been observed in other fluid-mechanical contexts by Walters and Forbes [27]. These errors are responsible for the small fluctuations on and about the central $xy$ -plane that may be visible in the last image for $t = 8$ on the right-hand side of Figure 7.

Figure 7 Outflow from the complex distributed source ( $f_{\ell } = 1$ ) described by (6.5), (6.6) for the second axisymmetric Legendre mode $n=2$ , $m=0$ , with negative amplitude $\varepsilon _2 = - 0.06$ . The interface is shown (left to right) at the four times $t = 2$ , $4$ , $6$ , $8$ .

The same situation is depicted in Figure 8 as for the previous diagram in Figure 7, except that the amplitude parameter is now positive, so that $\varepsilon _2 = 0.06$ . The numerical computations in this case are more difficult to perform accurately than for the corresponding situation in Figure 7 with $\varepsilon _2$ negative, and so the four results presented are for dimensionless times $t = 1.2$ , $2.4$ , $3.6$ , $4.8$ . For the last two times $t = 3.6$ and $4.8$ shown in Figure 8, small numerically generated oscillations are produced in the density perturbation function $\overline {\rho }$ near the z-axis and outside the source radius $r_S$ , where the expression (3.3) is most sensitive to errors in the Fourier coefficients $C_{m n p} (t)$ . These are responsible for the small perturbations visible in the last two times shown, but do not affect the interface shapes shown. The outflow morphology in this case is very different from that in Figure 7, since now the interfacial region forms a flatter torus-like object, pulling in slightly near the z-axis and lying parallel to the $xy$ -plane.

Figure 8 Outflow from the complex distributed source ( $f_{\ell } = 1$ ) described by (6.5), (6.6) for the second axisymmetric Legendre mode $n=2$ , $m=0$ , with positive amplitude $\varepsilon _2 = 0.06$ . The interface is shown (left to right) at the four times $t = 1.2$ , $2.4$ , $3.6$ , $4.8$ .

It seems at first surprising that a simple change in the sign of the amplitude parameter $\varepsilon _2$ could lead to such a large qualitative difference in the shape of the interface between Figures 7 and 8. An explanation can be obtained, however, by considering the radial velocity component $u_S^{\textrm {sph}}$ of the source in spherical polar coordinates, given in (6.1). For axisymmetric bipolar outflow, $n=2$ and $\ell = 0$ , this becomes

$$ \begin{align*} u_S^{\textrm{sph}} =\frac {r^2}{4\pi r_S^4} \bigg[ f_{\ell} - \varepsilon_2 \frac {3\pi}{r_S^2} ( 1 + 3\cos (2\phi ) ) \bigg] \end{align*} $$

inside the source region $r < r_S$ . When $\varepsilon _2$ is negative, and for the parameter values used here, this velocity component is positive on the z-axis where either $\phi = 0$ or $\phi = \pi $ , and negative on the $xy$ -plane where $\phi = \pi / 2$ . Thus, in Figure 7, the fluid flows outwards up and down the z-axis and inwards along the centre-plane $z = 0$ , giving rise precisely to the two axisymmetric lobes encountered in that sequence of diagrams. However, for the positive value $\varepsilon = 0.06$ used to generate Figure 8, the radial velocity component $u_S^{\textrm {sph}}$ is negative on the z-axis and positive on the centre-plane $z = 0$ . The fluid flows inwardly along the z-axis and outwardly near the $xy$ -plane, in an axisymmetric pattern, thus creating the flat platelet-shaped interfaces seen in Figure 8, with the pinched-in portions near the z-axis.

The outflow shown in Figure 9 corresponds to a source forced at the $(2,2)$ spherical mode with negative amplitude parameter $\varepsilon _2 = - 0.06$ in the source-flow velocity components (6.7) and (6.8). Solutions are shown at the four times $t = 2$ , $4$ , $6$ , $8$ and it is evident that the interface shape evolves from an initial sphere to a structure with two lobes along the x-axis. In this solution, minor numerical inaccuracies in the Fourier coefficients $C_{m n p} (t)$ , due to the enormous volume of computer arithmetic involved in this solution, again give rise to localized errors in the evaluation of the series (3.3) for the density perturbation function, in this instance near the y-axis, and these are picked up by our plotting procedure (3.20) and appear in Figure 9 as disturbances near that axis. Although irritating, they have no effect on the accuracy of the interface shape on the two lobes, however. Unlike the axisymmetric outflows considered in Figures 7 and 8, changing the sign of the amplitude parameter in Figure 9 does not cause major changes to the outflow morphology on this case; rather, the shapes depicted in Figure 9 are merely rotated by $90$ degrees about the z-axis.

Figure 9 Outflow from the complex distributed source ( $f_{\ell } = 1$ ) described by (6.7), (6.8) for the second Legendre mode $n=2$ , $m=2$ , with negative amplitude $\varepsilon _2 = - 0.06$ . The interface is shown (left to right) at the four times $t = 2$ , $4$ , $6$ , $8$ .

8 Conclusion

When a source of fluid is introduced into an ambient fluid of different density, an interfacial zone forms between the fluids; as time progresses, the interfacial zone is driven outwards and can form complex shapes. This paper is concerned with understanding the conditions that affect these outflow morphologies. In pursuit of this aim, we have developed a novel approach to the description of the source region, representing it as a sphere of finite radius $r_S$ within which a background velocity is prescribed such that the divergence $\textrm {div}\, \mathbf {q}$ of the velocity field is nonzero and creates mass within that region only. Furthermore, the integral of the divergence over the spherical source region gives the required volume flux $f_{\ell }$ . We use the bi-streamfunction approach of Forbes and Brideson [13] and Walters and Forbes [27] in a Boussinesq model of the fluids to formulate an efficient spectral method that solves for the density and the fluid velocity vector; nevertheless, in unsteady three-dimensional flow, this method still requires enormous computational resources, and we use a high-performance computer with high-level parallelizing capabilities based around the use of graphics-processing units.

Various linearized stability analyses have been investigated in this paper for two ideal inviscid fluids separated by an infinitessimally thin interface. In particular, we have attempted in Section 5 to estimate the effect that a source of finite radius $r_S$ might have on the eventual shape of the outflow region, although we were only able to do this for early times. In addition, the boundary condition (5.1), crucial in this analysis, is nevertheless not the same as that implicit to the full numerical scheme in Section 3, and so results such as those in Figures 1 are at best only a guide to the likely behaviour of the full system. There is clearly considerable scope for future developments in this analysis.

Numerical solutions have been presented here using our spectral algorithm that can accommodate fully three-dimensional time-dependent outflow from the distributed source. A major aim of this present work has been to determine to what extent the outflow, from the source of finite radius $r_S$ , retains its memory of initial conditions or whether it ultimately achieves some different long-term configuration. In particular, Forbes [10] had argued that many flows from a spherical point source should eventually produce an interface with a one-sided jet, in which the first spherical model $n = 1$ ultimately dominates, although Forbes [11] and Forbes and Brideson [12] later showed how flow conditions in the surrounding heavier fluid could alter that conclusion. The effect of the ambient medium upon outflow conditions has been considered both theoretically and experimentally in an astrophysical situation in a recent article by Schulreich and Breitschwerdt [24].

The results presented here suggest that the relationship between the outflow shape and the initial configuration are possibly somewhat more nuanced than Forbes [10] indicated. When a small perturbation is made at the first spherical mode, a strongly one-sided outflow does indeed develop, as indicated in Figure 2. After some time has elapsed, the single outflow jet also overturns and forms a rotationally symmetric mushroom shape. In some other instances, a one-sided outflow does also evidently become dominant, such as for the mode- $3$ example illustrated in Figure 6, in which an asymmetric shape evolves, with a strong outflow at the top of the diagram and a far weaker one at the bottom. However, in some other instances illustrated in Section 7, the situation is much less clear, with apparently symmetric bipolar lobes forming up and down the z-axis, as in Figure 7, or along some other axis, as illustrated in Figure 9.

The question still remains, however, to what extent an apparently symmetric outflow (illustrated in Figure 7, for example) might nevertheless eventually collapse onto a one-sided configuration if its initial condition were not simply a pure mode- $2$ perturbation. Given the enormous computing resources needed to generate these results, this question cannot be answered comprehensively here, but we have attempted to address this in Figure 10. Here, a gravity-free solution has been computed $(1/F^2 =~0)$ with the slightly smaller value $r_S = 0.22$ , where the approximate stability analysis of Section 5 would suggest that only the $n = 1$ mode would in fact be unstable, as in Figure 1. In Figure 10, the initial condition consisted of a velocity perturbation at the second Fourier-Legendre mode $n=2$ with amplitude $\epsilon = 0.025$ , with an additional much weaker disturbance at the first mode $n=1$ having amplitude $\epsilon = 0.0025$ . The results in Figure 10 develop very differently from those in Figure 7; for example, now, although the second mode is clearly dominant at the first time $t = 2.5$ presented in Figure 10, at the next time $t = 5$ shown, a significant asymmetry has developed with a stronger outflow at the top of the diagram than at the bottom. At the third time $t = 7.5$ indicated, the outflow morphology is evidently now dominated by the first mode and this pattern continues in the final time $t = 10$ presented. In addition, the outflow jets have not yet developed overturning structures at their tips, so that they appear more similar to the outflows reported by Forbes [10]. Thus, while the Forbes suggestion that one-sided jets eventually dominate may be true for complex, multi-mode initial conditions, it is not yet clear whether this is always the case or if so, how long it might take for the one-sided jet to become apparent. There is clearly more to be understood on this intriguing question.

Figure 10 Outflow with infinite Froude number ( $1/F^2 = 0$ ), perturbed primarily at the second axisymmetric Legendre mode $n=2$ , $m=0$ , but with a small additional contribution at the first mode. The interface is shown (left to right) at the four times $t = 2.5$ , $5$ , $7.5$ , $10$ .