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 [Reference Forbes and Brideson13] and Walters and Forbes [Reference Walters and Forbes27], 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 [Reference Forbes10], 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. [Reference Lin, Gerya, Tackley, Yuen and Golabek17] 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 [Reference Abramowitz and Stegun1, 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 [Reference Forbes10, Reference Forbes11] 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 [Reference Epstein8]. 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. [Reference Peng, Zabusky and Zhang22]), 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 [Reference Forbes10]. 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 [Reference Forbes10] 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 [Reference Abramowitz and Stegun1, 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 [Reference Forbes10] 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 [Reference Forbes10], 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 [Reference Forbes10].

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 [Reference Terrones and Heberling26] 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 [Reference Forbes10, Reference Forbes11], 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 [Reference Batchelor5, 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.