2.1. Similarity model of a pure plume

The fluid in a plume moves under gravity and is driven by a density difference between the plume fluid and its surrounding environment. The plume body exhibits significant mean vertical velocity $w$ and reduced gravity $g'$. The reduced gravity $g'$ is defined as $g' = {g(\rho _a-\rho )}/{\rho _0}$, where $\rho _a$ and $\rho$ denote the density of the ambient fluid and the plume, respectively, $\rho _0$ is some reference density and $g$ is the acceleration due to gravity. The radial profiles of these time-averaged quantities are approximately Gaussian, and may be represented as

(2.1)\begin{gather} \tilde{w} = C_w{F_0}^{1/3} z^{{-}1/3}\exp{({-}r^2/\tilde{b}^2)}, \end{gather}

(2.2)\begin{gather}\widetilde{g'} = C_{g'}{F_0}^{2/3} z^{{-}5/3}\exp{({-}r^2/\tilde{b}^2)}, \end{gather}

where $\tilde {w}(r,z)$ is the vertical velocity, $\widetilde {g'}(r,z)$ is the reduced gravity, $\tilde {b}(z)$ is the $e$-folding length, $r$ is the radial or cross-plume coordinate, $z$ is the vertical coordinate, and $C_w$ and $C_{g'}$ are constants. In the self-similar regime, $\tilde {b} \propto z$. A small difference in the $e$-folding length scale between the two profiles is often observed, but also often assumed negligible so that the approximation of a single length scale is made, as above.

Such buoyant plumes are typically modelled with integral models. The classical model of Morton, Taylor & Turner (Reference Morton, Taylor and Turner1956) comprises a set of three coupled ordinary differential equations, which are based on three main assumptions: self-similarity, Boussinesq approximation and the constant entrainment coefficient $\alpha$. In an unstratified environment the equations for volume flux $Q$, momentum flux $M$ and buoyancy flux $F$ are

(2.3a–c)\begin{equation} \frac{\text{d}Q}{\text{d}z}=2 \alpha M^{1/2} , \quad \frac{\text{d}M}{\text{d}z} = 2 \frac{QF}{M} , \quad \frac{\text{d}F}{\text{d}z} = 0. \end{equation}

The ‘top-hat’ formulation replaces the peak values of $w$ and $g'$ and the Gaussian $e$-folding length scale with horizontally averaged values of $w$ and $g'$ and an associated mean radius $b$. The fluxes can then be expressed as

(2.4a–c)\begin{equation} Q = b^{2}w, \quad M=b^{2}w^{2}, \quad F = b^{2}wg'. \end{equation}

From here on, we use $w(z)$ and $g'(z)$ to denote the top-hat variables, replacing the Gaussian profiles of $\tilde {w}(r,z)$ and $\widetilde {g'}(r,z)$, respectively.

In the formulation of these equations it is assumed that the plume originates from a point source. In actuality, plumes will begin from a source with some finite radius. To account for the finite nature of the source, one can extrapolate the sides of the plume back to some point origin, as shown in figure 1(a). This virtual origin correction is made by taking the height as $z+z_v$ (when calculating the plume radius by the similarity solution (2.5), for instance), where $z_v$ is the height of the finite source above the virtual point source. At larger heights when $z \gg z_v$, the virtual origin correction will be negligible.

Figure 1. The configuration of two plumes with vertical offset $\zeta$ and horizontal source separation $\chi$ shown in the $(a)$ $xz$ plane and $(b)$ $xy$ plane. The ‘lower’ plume is on the left and the ‘upper’ plume is on the right. The origin of coordinates $O$ is positioned along the centreline of the lower plume at the height of the upper plume virtual origin such that the virtual origin of the lower plume is at $(0,0,-\zeta )$, and that of the upper is at $(\chi,0, 0)$. The shaded circular region shows the cross-sectional area of the lower plume at $z=0$.

Under the assumption that the ambient fluid is of uniform density ($\rho _0 \equiv \rho _a$), Morton et al. (Reference Morton, Taylor and Turner1956) give a similarity solution for an axisymmetric pure plume as

(2.5)\begin{gather} b = \frac{6\alpha}{5}z, \end{gather}

(2.6)\begin{gather}w = \frac{5}{6\alpha}\left(\frac{9}{10}\alpha F\right)^{1/3}z^{{-}1/3} , \end{gather}

(2.7)\begin{gather}g' = \frac{5F}{6\alpha}\left(\frac{9}{10}\alpha F\right)^{{-}1/3} z^{{-}5/3}. \end{gather}

2.2. The sink model of entrainment

A line sink of strength $-m(z)$ positioned at the origin has complex potential

(2.8)\begin{equation} \varOmega ={-}\frac{m}{2{\rm \pi}}\ln Z, \end{equation}

where $Z = x+{\rm i}y = r\text {e}^{\text {i}\theta }$. Since the flow induced by plume entrainment is horizontal and irrotational, Kaye & Linden (Reference Kaye and Linden2004) used line sinks to represent the flow exterior of two interacting plumes.

The principle of superposition implies that the complex potentials of individual line sinks can be summed to give a total complex potential representing the entrainment field of an equivalent configuration of interacting plumes. The real part of the resulting complex potential $\varOmega = \phi + \text {i}\psi$ corresponds to the velocity potential.

As described previously by Rooney (Reference Rooney2015) and Rooney (Reference Rooney2016), contours of equal velocity potential can be used to approximate the mean plume boundaries at any height. Outside the plume boundaries, the streamfunction from the complex potential yields the combined horizontal entrainment field of multiple plumes. The vertical evolution of the system is controlled by the buoyancy of the plumes and the entrainment along the perimeter of the plume system. The distortion and merging of the plume boundaries is mirrored by that of the velocity-potential contours (e.g. Rooney Reference Rooney2016, figure 7).

3.1. A modified circle theorem

The circle theorem of Milne-Thomson (Reference Milne-Thomson1940) expresses how an irrotational, two-dimensional flow of incompressible and invisid fluid is altered upon the insertion of a circular cylinder into the flow. The circle of radius $R$ obtained from a cross-section of the cylinder is denoted $C$ and has the equation $\lvert Z \rvert = R$. The theorem applies when the complex potential of the initial flow $f(Z)$ has no rigid boundaries and no singularities in the region $\lvert Z \rvert \leqslant R$ for $R \in \mathbb {R}$. Placing the circle into this flow changes the complex potential to

(3.1)\begin{equation} \varOmega(Z) = f(Z) + \overline{f\left({R^2}/{\bar{Z}}\right)}, \end{equation}

where the bar denotes the complex conjugate. Since $R^2/\bar {Z}$ is the reflection of $Z$ in $C$ and $\lvert Z \rvert > R$ implies $\lvert R^2/\bar {Z}\rvert < R$, it follows that the additional, perturbing term has no singularities outside of $C$. Furthermore, the circle is a streamline. This can easily be illustrated by considering $Z=R\text {e}^{\text {i}\theta }$ such that the complex potential becomes

(3.2)\begin{equation} \varOmega(R\text{e}^{\text{i}\theta}) = f(R\text{e}^{\text{i}\theta}) + \overline{f\left(\frac{R^2}{R\text{e}^{-\text{i}\theta}}\right)} = f(R\text{e}^{\text{i}\theta}) + \overline{f(R\text{e}^{\text{i}\theta})}, \end{equation}

and, hence, is entirely real.

To instead identify the circle as a velocity-potential contour, the theorem must be modified to ensure a purely imaginary potential on $C$. Similar to the original theorem, we require the perturbing term to have all its singularities within the circle. The complex potential

(3.3)\begin{equation} \varOmega(Z) = f(Z) - \overline{f\left({R^2}/{\bar{Z}}\right)} \end{equation}

satisfies this condition, and is entirely imaginary on $C$.

As a first example, we consider uniform flow around a circle of radius $R$ centred at the origin. Taking $f(Z) = Z$, the complex potential given by the original circle theorem is

(3.4)\begin{align} \varOmega &= f(Z) + \overline{f(R^2/\bar{Z})} \nonumber\\ &= r\text{e}^{\text{i}\theta} + \frac{R^2}{r}\text{e}^{-\text{i}\theta} \end{align}

using $Z = r\text {e}^{\text {i}\theta }$. The imaginary part of the complex potential corresponds to the streamfunction $\psi$, where

(3.5)\begin{equation} \psi(r, \theta) = \sin\theta\left( r - R^2/r\right). \end{equation}

A streamline is a line everywhere tangent to the local fluid velocity and, hence, $\psi$ remains constant along streamlines. The streamlines corresponding to the streamfunction in (3.5) are plotted in figure 2(a). As expected, the boundary $\lvert Z \rvert = R$ corresponds to a level set of the streamfunction as the potential is entirely real on the circle (imaginary part is constant).

Figure 2. Streamlines for uniform flow around a circle of unit radius calculated with the Milne–Thomson circle theorem and the adapted theorem are shown in (a,b), respectively.

Applying the adapted circle theorem, the streamfunction becomes

(3.6)\begin{equation} \psi(r, \theta) = \sin\theta\left( r + R^2/r\right), \end{equation}

and the corresponding streamlines are plotted in figure 2(b). The streamlines demonstrate that there is now flow both into and out of the circle perpendicular to its boundary such that the total net flow is zero due to symmetry. Since streamlines are orthogonal to velocity-potential contours, this flow pattern highlights the circular velocity-potential contour at $\lvert Z \rvert = R$ that forms as the complex potential is entirely imaginary on the circle.

3.2. Vertically offset sources

A line (or point) sink may be used to represent a plume from a point source. To represent sources with a vertical as well as a horizontal offset, a model is developed to represent a two-plume system like that shown in figure 1(b). This requires two additions to the domain containing the line sink. The first is the insertion of a circular velocity-potential contour as the boundary of the lower plume, using the adapted circle theorem. This aligns the streamlines normal to the boundary, as in figure 2(b). The second is to represent the entrainment of the lower plume via a second sink, concentric with the circle. The additional sink being concentric means that the streamline alignment at the circle boundary does not change due to its presence. However, there is now a net inflow across the boundary, representing entrainment of the lower plume. The magnitude of the inflow is determined by the strength of the second sink.

Coelho & Hunt (Reference Coelho and Hunt1989) discussed the related topic of the flow field external to a jet in a crossflow, and described that as being similar to potential flow around a cylinder with suction. This method of modelling the lower plume may be described likewise.

To insert a circular potential contour, we can apply the adapted theorem to the aforementioned cross-section at $z=0$. The complex potential for a line sink of strength $-m(z)$ at distance $\chi$ along the $x$ axis is given as

(3.7)\begin{equation} f(Z) ={-}\frac{m}{2{\rm \pi}}\ln(Z-\chi). \end{equation}

The adapted circle theorem requires that $R<\chi$ such that the singularity of $f(Z)$ lies outside of the circle. Under this assumption, the total complex potential is

(3.8)\begin{equation} \varOmega(Z) ={-}\frac{m}{2{\rm \pi}}(\ln(Z-\chi) - \ln(R^2/Z-\chi)), \end{equation}

using $\overline {\ln (Z)} = \ln \bar {Z}$. Non-dimensionalisation by the horizontal source separation $\chi$ yields

(3.9)\begin{equation} \varOmega(Z^{*}) ={-}\frac{m}{2{\rm \pi}}\ln\left(\frac{Z^{*}(Z^{*}-1)}{\gamma^2-Z^{*}}\right), \end{equation}

where $Z^{*} = Z/\chi$ and $\gamma =R/\chi$ is a dimensionless constant representing the ratio of the lower plume radius at $z=0$ to the horizontal source separation. For $Z \in \mathbb {C}$ with $Z = x + \text {i}y$,

(3.10)\begin{equation} \ln Z = \frac{1}{2}\ln{(x^2+y^2)}+\text{i}\arctan\left(\frac{y}{x}\right). \end{equation}

Hence, the complex potential in (3.9) can be written as

(3.11)\begin{equation} \varOmega(Z^{*}) ={-}\frac{m}{2{\rm \pi}}\ln\left(\frac{x^{*}(x^{*}-1)-y^{*2}+\text{i}\left(y^{*}(x^{*}-1)+x^{*}y^{*}\right)}{\gamma^2-x^{*}-\text{i}y^{*}}\right). \end{equation}

The streamfunction is obtained from the imaginary part as

(3.12)\begin{equation} \psi(x^{*},y^{*}) ={-}\frac{m}{2{\rm \pi}}\arctan\left(\frac{y^{*}\left(\gamma^2(2x^{*}-1)-(x^{*2}+y^{*2})\right)}{\gamma^2(x^{*}(x^{*}-1)-y^{*2})-x^{*}(x^{*2}+y^{*2})+x^{*2}+y^{*2}}\right), \end{equation}

where asterisks are used to denote non-dimensional quantities.

The corresponding streamlines plotted in figure 3(a) demonstrate that flow leaves the circle in the direction of the upper plume source. To increase the flow into the circle and, hence, create a circular sink, an additional sink is required at the centre of the circle to represent the entrainment of the lower plume. We label the sinks such that the original sink at a distance $\chi$ along the $x$ axis has strength $-m_1(z)$ whilst the further sink at the origin has strength $-m_2(z)$.

Figure 3. Streamlines for sink flow around a circle calculated by the adapted circle theorem for (a) no additional sink at the origin, and (b) an additional sink of equal strength placed at the origin. The black circles represent the radii of the plumes given by (2.5) at height $z^{*}=0$. At this height, the lower plume has radius $R$ and the upper plume is a point source (not to scale). The stagnation point in (b) is marked with a red cross.

Upon addition of the second sink, the complex potentials in (3.8) and (3.9) become

(3.13)\begin{equation} \varOmega(Z) ={-}\frac{m_1}{2{\rm \pi}}\left[\ln(Z-\chi) - \ln(R^2/Z-\chi)\right] - \frac{m_2}{2{\rm \pi}}\ln Z \end{equation}

and

(3.14)\begin{equation} \varOmega(Z^{*}) ={-}\frac{m_2}{2{\rm \pi}}\ln\Bigg[\left(\frac{(Z^{*}-1)Z^{*}} {\gamma^2-Z^{*}}\right)^MZ^{*}\chi\Bigg], \end{equation}

respectively, where $M=m_1/m_2$ is the ratio of sink strengths. The equation for the streamfunction changes similarly as

(3.15) \begin{align} \psi(x^{*},y^{*}) &={-}\frac{m_1}{2{\rm \pi}}\arctan\left(\frac{y^{*}(\gamma^2(2x^{*}-1)-(x^{*2}+y^{*2}))}{\gamma^2(x^{*}(x^{*}-1)-y^{*2})-x^{*}(x^{*2}+y^{*2})+x^{*2}+y^{*2}}\right) \nonumber\\ &\quad - \frac{m_2}{2{\rm \pi}} \arctan\left({\frac{y^{*}}{x^{*}}}\right). \end{align}

An example with equal sink strengths ($M=1$) is plotted in figure 3(b). The additional sink at the origin means that a greater proportion of the flow is now directed into the lower plume. As a result, a stagnation point forms between the two plumes and its position is controlled by the ratio of sink strengths $M$.

As mentioned above, Coelho & Hunt (Reference Coelho and Hunt1989) also pointed out the similarity between an entraining jet or plume and a cylinder with suction at the surface. The method presented above (applying the modified circle theorem plus an additional sink to represent a circular sink of known strength) has applications for modelling other situations. One example would be the insertion of a porous cylinder into a two-dimensional flow for the purpose of extracting fluid. For example, a sink concentric with the circular boundary in figure 2(b) would have the effect of drawing some of the ambient flow into the cylinder, altering the flow field accordingly. If the additional sink were replaced by a source, that would alternatively represent fluid extrusion through a porous cylinder. Sucking or blowing at the boundary of cylindrical obstacles is one method of modifying fluid–structure interaction to reduce drag or to increase flow stability for instance (Fransson, Konieczny & Alfredsson Reference Fransson, Konieczny and Alfredsson2004; Chen et al. Reference Chen, Huang, Chen, Yu and Gao2022). It is possible that a suitable coordinate transformation may allow application of this method using the circle theorem to represent a sink of non-circular geometry (for example, this would be straightforward for an elliptical sink).

3.3. Entrainment in offset sources

To calculate the sink strengths at any given height, we take the sink strength to be equal to the entrainment across any velocity-potential contour, which is shown by the entrainment-flux integral calculations of Rooney (Reference Rooney2016). Since the entrainment constant $\alpha$ is selected as in Morton et al. (Reference Morton, Taylor and Turner1956), the entrainment flux into the plume at any particular height is equivalent to the sink strength and is defined by (Kaye & Linden Reference Kaye and Linden2004)

(3.16)\begin{equation} m = 2{\rm \pi}\alpha b w. \end{equation}

Up until height $z=0$, the lower plume is assumed to be unaffected by the upper plume source. Therefore, it is axisymmetric and the total entrainment into its circular area at this height can be calculated simply with (3.16). The similarity solution of Morton et al. (Reference Morton, Taylor and Turner1956) implies that

(3.17)\begin{gather} b \sim \alpha z \end{gather}

(3.18)\begin{gather}w \sim \alpha^{{-}2/3}F^{1/3}z^{{-}1/3}, \end{gather}

where $F$ is the buoyancy flux of the plume.

When a small height above the upper plume point source is considered, both plumes can be assumed to remain approximately axisymmetric (Li & Flynn Reference Li and Flynn2021). As the height above the plume source increases, the plumes will distort so that there is no single value of radius. However, in the region under consideration, this distortion is initially small (as will be evident from RANS results in later sections) and the plumes will be assumed to have a circular cross-section, the radius of which can be calculated by the similarity solution (2.5).

Using the subscripts $1$ and $2$ to correspond to the values of the upper and lower plumes, respectively, the entrainment of the upper plume at a height $z$ above the upper plume is given as

(3.19)\begin{equation} m_1 \sim \alpha^{4/3}F_1^{1/3} z^{2/3}. \end{equation}

For a vertical offset of $\zeta$, the entrainment of the lower plume at this height is

(3.20)\begin{equation} m_2 \sim \alpha^{4/3}F_2^{1/3}(z + \zeta)^{2/3}, \end{equation}

and, hence, the ratio of entrainment is

(3.21)\begin{equation} M = \left(\frac{F_1}{F_2}\right)^{1/3}\left(\frac{z }{z + \zeta }\right)^{2/3}. \end{equation}

This indicates that $M$ is zero at $z=0$, i.e. at the level of the point source. Additionally, in an unstratified environment it is assumed that the buoyancy flux does not vary with height and, hence, the ratio $F_1/F_2$ at all heights $z$ is equal to the initial ratio of buoyancy fluxes for the two plumes, so $M \rightarrow (F_1/F_2)^{1/3}$ for large $z$.

3.4. Flow speed on contours

From the complex potential in (3.14), the velocity field may be obtained as $\text {d}{\varOmega }/\text {d}Z = U-\text {i}V$. For vertically offset plumes, this is

(3.22)\begin{align} \frac{\text{d}{\varOmega}}{\text{d}{Z}} &= \frac{1}{\chi}\frac{\text{d}{\varOmega}}{\text{d}{Z^{*}}} \end{align}

(3.23)\begin{align} &={-}\frac{m_2}{2{\rm \pi}\chi} \frac{(M+1)Z^{*2}-(\gamma^2(2M+1)+1)Z^{*}+(M+1)\gamma^2}{(Z^{*}-\gamma^2)(Z^{*}-1)Z^{*}}, \end{align}

so that

(3.24)\begin{equation} U ={-}\frac{m_2}{2{\rm \pi}\chi}\frac{N}{D} \text{,} \end{equation}

where

(3.25)\begin{align} N &= (M+1)[(\gamma^{2}+x^{*2}-y^{*2})(1-x^{*})(\gamma^{2}-x^{*})x^{*}] \nonumber\\ &\quad -x^{*2}(\gamma^{2}-x^{*})(1-x^{*})(\gamma^{2}(2M+1)+1)\nonumber\\ &\quad +y^{*2}(M+1)(\gamma^{2}+1-3x^{*})(\gamma^{2}+x^{*2}-y^{*2})\nonumber\\ &\quad -x^{*}y^{*2}(\gamma^{2}+1-3x^{*}) (\gamma^{2}(2M+1)+1)+ 2x^{*2}y^{*2}(M+1)(3x^{*}-2\gamma^{2}-2)\nonumber\\ &\quad +2x^{*}y^{*2}(M+1)(\gamma^{2}-y^{*2})-x^{*}y^{*2}(3x^{*}-2\gamma^{2}-2)(\gamma^{2}(2M+1)+1) \nonumber\\ &\quad -y^{*2}(\gamma^{2}(2M+1)+1)(\gamma^{2}-y^{*2}), \end{align}

(3.26)\begin{align} D &= x^{*2}\gamma^4 - 2x^{*3}\gamma^2(\gamma^2+1)+x^{*4}(1+4\gamma^2+\gamma^4)-2x^{*5}(\gamma^2+1)+x^{*6} \nonumber\\ &\quad+y^{*4}(\gamma^4+1)+ x^{*2}y^{*2}(3(x^{*2}+y^{*2})+2\gamma^2(\gamma^2+1)+2(\gamma^2+1)-4x^{*}(1+\gamma^2)) \nonumber\\ &\quad-2x^{*}y^{*4}(\gamma^2+1)-2\gamma^2x^{*}y^{*2}(\gamma^2+1)+y^{*2}(y^{*4}+\gamma^4). \end{align}

To find the position of the stagnation point between the two plumes at height $z=0$, we note that the velocity component $V$ goes to zero along the line $y^{*}=0$ that connects the two plumes. Hence, the flow speed along this line is entirely in the $x$ direction and can be found by evaluating

(3.27)\begin{equation} U|_{y^{*}=0} = \frac{(M+1)(\gamma^2+x^{*2})-x^{*}(\gamma^2(2M+1)+1)}{x^{*}(x^{*}-1)(x^{*}-\gamma^2)}. \end{equation}

To find the non-dimensional distance of the stagnation point from the lower plume centreline $x_{s}^{*}$, $U_{y^{*}=0}$ is set equal to zero and the equation is solved for $x^{*}$ to obtain

(3.28)\begin{equation} x^{*}_{s} = \frac{\gamma^2(2M+1)+1\pm \sqrt{(\gamma^2(2M+1)+1)^2 - 4\gamma^2(M+1)^2}}{2(M+1)}. \end{equation}

Physically, the stagnation point must lie between the two plumes, otherwise one plume would be forced to behave as a source rather than a sink. Additionally, the stagnation point is valid only when it lies outside of the radii of the plumes. The solution given by the negative square root of (3.28) lies inside the radius of the lower plume and can therefore be disregarded. If the stagnation point given by the positive square root also lies within the lower plume radius, then plume detrainment occurs. The limiting case is defined where the stagnation point lies at $x^{*}_{s} = \gamma$, i.e. on the radius of the circular sink at the point closest to the upper plume source. This value can be substituted into the positive square root case of (3.28) to find the corresponding critical ratio of sink strengths $M_c$ as

(3.29)\begin{equation} M_c = \frac{1-\gamma}{2\gamma}. \end{equation}

The value $M_c$ describes the maximum ratio of sink strengths that can be reached before detrainment of the lower plume occurs and the flow field becomes unrealistic. If $M \leqslant M_c$, the sink at the origin is sufficiently strong to overcome detrainment and the lower plume continues to act as a sink at all points on its boundary. In this case, a valid stagnation point forms and its position is given by taking the positive square root in (3.28). Letting $(M+1)^{-1} = \phi \leqslant 1$, then (3.28) may be rewritten as

(3.30)\begin{equation} x^{*}_{s} = \frac{\phi}{2} + \gamma^2 \left(1 - \frac{\phi}{2} \right) \pm \sqrt{ \left( \frac{\phi}{2} + \gamma^2 \left(1 - \frac{\phi}{2} \right) \right)^2 - \gamma^2 }, \end{equation}

which helps demonstrate that $x^{*}_{s} \rightarrow 1$ as $M \rightarrow 0$.

When $M>M_c$, there is no stagnation point and the flow from the lower plume is directed into the upper plume such that the lower plume will eventually be absorbed by the upper plume. This is more likely to occur if the upper plume is ‘stronger’ (has a greater buoyancy flux) than the lower plume. In this case, the lower plume cannot maintain its own structure. The model therefore indicates the conditions under which the lower plume loses its integrity, and detrains into the upper plume. It should be noted that, at and beyond this point, the model described here will no longer represent the plume flow structure, since the flow has undergone a qualitative change. Streamlines for examples of different sink strength ratios $M$ are shown in figure 4.

Figure 4. Streamlines calculated by the adapted circle theorem with an additional sink positioned at the origin and $M_c = 2$. The ratio of sink strengths are (a) $M = 1$, (b) $M = 2$ and (c) $M = 3$. The stagnation point is marked with a red cross in (a,b), but does not exist outside of the plume boundaries in (c). Therefore, (c) depicts a flow field that would be unrealistic for merging offset plumes.

4.1. The governing equations

We used ANSYS Fluent Academic Research Mechanical and CFD Release 2020 R2 (ANSYS Inc., Canonsburg, PA, USA) to perform RANS simulations. The steady RANS equations for mass and momentum conservation in an incompressible fluid, solved in a three-dimensional Cartesian coordinate system are

(4.1)\begin{gather} \frac{\partial (\rho u_i)}{\partial x_i} = 0, \end{gather}

(4.2)\begin{gather}\frac{\partial (\rho u_i u_j)}{\partial x_j} ={-}\frac{\partial p}{\partial x_i} + \frac{\partial}{\partial x_j}\left(\tau_{ij} - \rho \overline{u^{\prime}_{i} u^{\prime}_{j}} \right) + (\rho - \rho_a) g_i, \end{gather}

where $u_i$ are the mean velocity components, $u_i^\prime$ are the fluctuating velocity components, $p$ is pressure, $\rho$ is density, $\rho _a$ is ambient density and $g_i = (0,0,-g)$ in the $x$, $y$ and $z$ directions, where $g$ is the magnitude of acceleration due to gravity. The Reynolds stresses $\rho \overline {u^{\prime }_{i} u^{\prime }_{j}}$ are modelled using turbulence models to close the system of equations. The deviatoric stress tensor $\tau _{ij}$ is defined as

(4.3)\begin{equation} \tau_{ij} = \mu\left(\frac{\partial u_j}{\partial x_i} + \frac{\partial u_i}{\partial x_j}\right) - \frac{2}{3}\mu\frac{\partial u_k}{\partial x_k}\delta_{ij}, \end{equation}

where subscript $i,j,k= 1,2\text { and }3$ and $\mu$ is dynamic viscosity.

Using the Boussinesq approximation, we assume small temperature differences and consider only density variations in the buoyancy term in the momentum equation

(4.4)\begin{equation} (\rho - \rho_a)g \approx{-}\rho_a\beta(T-T_a)g, \end{equation}

where $\beta$ is the coefficient of thermal expansion, $T$ is a reference temperature and $T_a$ is the constant ambient temperature. Therefore, the density $\rho$ can be expressed as

(4.5)\begin{equation} \rho = \rho_a \left (1-\beta(T-T_a) \right ) . \end{equation}

To close the system of equations, we use a $k$-$\epsilon$ turbulence model that is known to be robust in modelling a range of turbulent flows. In particular, we use the renormalisation group (RNG) $k$-$\epsilon$ turbulence model, based on its success in modelling a single plume (Cook & Lomas Reference Cook and Lomas1998) and multiple plumes (e.g. Mokhtarzadeh-Dehghan et al. Reference Mokhtarzadeh-Dehghan, König and Robins2006; Durrani et al. Reference Durrani, Cook, McGuirk and Kaye2011). The turbulent kinetic energy $k$ and dissipation rate $\epsilon$ transport equations for the RNG $k$-$\epsilon$ turbulence model are given as

(4.6)\begin{gather} \frac{\partial (\rho k u_i)}{\partial x_i} =\frac{\partial }{\partial x_j} \left(\alpha_k \mu_{t} \frac{\partial k }{\partial x_j} \right) +G_k +G_b - \rho \epsilon, \end{gather}

(4.7)\begin{gather}\frac{\partial (\rho \epsilon u_i)}{\partial x_i} =\frac{\partial }{\partial x_j} \left(\alpha_\epsilon \mu_{t} \frac{\partial \epsilon }{\partial x_j} \right) + C_{1\epsilon} \frac{\epsilon}{k}(G_k+ C_{3 \epsilon} G_b)- C_{2 \epsilon}^* \frac{\epsilon^2}{k} , \end{gather}

with turbulent viscosity $\mu _t = \rho C_\mu k^2/\epsilon$. The inverse effective Prandtl numbers $\alpha _k$ and $\alpha _\epsilon$ correspond to the turbulent kinetic energy $k$ and dissipation rate $\epsilon$, respectively, and their values are obtained from the RNG theory by Orszag et al. (Reference Orszag, Yakhot, Flannery, Boysan, Choudhury, Maruzewski and Patel1993).

Here $G_k$ is the generation of turbulent kinetic energy $k$ due to the mean velocity gradients and is defined as

(4.8)\begin{equation} G_k = \mu_t(2S_{ij}S_{ij})^{1/2} , \end{equation}

where the strain rate tensor $S_{ij}$ is

(4.9)\begin{equation} S_{ij} = \frac{1}{2} \left( \frac{\partial{u_i}}{\partial{x_j}} + \frac{\partial{u_j}}{\partial{x_i}} \right) . \end{equation}

The effect of buoyancy on turbulence was established to be significant by Nam & Bill (Reference Nam and Bill1993). Therefore, the generation of turbulent kinetic energy due to buoyancy $G_b$ in (4.6) and (4.7) is modelled as

(4.10)\begin{equation} G_b = \beta g \frac{\mu_t}{\textit{Pr}_t} \frac{\partial T}{\partial x_i} , \end{equation}

with turbulent Prandtl number for energy $\textit {Pr}_t = 0.85$.

The values for the model constants are $C_{1\epsilon }=1.42$, $C_{2\epsilon }^*=2.0$ and $C_\mu =0.0845$. The constant that controls the effect of buoyancy in (4.7) is modelled according to Henkes, Van Der Vlugt & Hoogendoorn (Reference Henkes, Van Der Vlugt and Hoogendoorn1991) and is defined as

(4.11)\begin{equation} C_{3\epsilon} = \tanh{| v/u |}, \end{equation}

where $v$ is the vertical velocity component parallel to the gravitational vector and $u$ is the horizontal velocity component perpendicular to the gravitational vector.

Finally, the energy transport equation is modelled as

(4.12)\begin{equation} \frac{\partial (\rho u_i h)}{\partial x_i} =\frac{\partial }{\partial x_j} \left[\left( \kappa + \frac{c_p \mu_t}{{\textit{Pr}}_t}\right) \frac{\partial T}{\partial x_j} \right] , \end{equation}

where the thermal conductivity $\kappa =0.0248\ \text {W}\ (\text {m} \text {K})^{-1}$, specific heat capacity at constant pressure $c_p=1006.43 \text {J}\ (\text {kg}\ \text {K})^{-1}$ and $h$ is the specific enthalpy.

4.2. Computational domain and boundary conditions

A schematic of the computational domain and mesh for two vertically offset plume sources is shown in figure 5. The mesh shown is indicative of the Cartesian mesh used, but is much more coarse. Additionally, the actual source is much smaller compared with the domain, but was enlarged to make its geometry and the boundary conditions visible on the figure. To reduce the number of grids and computational complexity, we use an $xz$ symmetry plane at $y=0$ and only model half of the plumes since the problem is symmetric about $y=0$. The two plumes are in close proximity to each other, therefore allowing for plume interaction in the domain. Both plumes have the same source length scale $d=0.02$ m, and the lower plume source is located at the bottom boundary $z= -\zeta$, where $\zeta$ is the vertical offset. To model the higher plume, we extrude a prism of dimensions $d \times d/2 \times \zeta$ from the base domain, which is located at a horizontal distance $\chi$ from the lower plume. The size of the domain measured in multiples of the plume source length scale is $30d\times 15d \times 50 d$ in the $x$, $y$ and $z$ directions, respectively. To prevent large-scale recirculation and instabilities, the size of the domain is relatively large compared with the length scale of the plume source. The horizontal extent of the domain in the $x$ direction of $30d$ is comparable to the set-up for two interacting plumes used by Durrani et al. (Reference Durrani, Cook, McGuirk and Kaye2011), where their domain is $40d$ in the $x$ and $y$ directions. Our domain has a ratio of $z_{max} / x_{max} =1.6$, which is similar to Pham, Plourde & Doan (Reference Pham, Plourde and Doan2007). We also ensured that the domain is high enough to allow the plumes to merge, thus preventing the influence of the top boundary conditions on the near-field region.

Figure 5. Computational domain and mesh with $xz$ symmetry plane at $y=0$. The boundary conditions are colour coded. Grey, pressure inlet at ambient temperature for the bottom and side planes; red, pressure inlet at source temperature for plume sources; blue, pressure outlet for top plane; green, no-slip wall for plume offset. (Note: this is not to scale, the grid cells and plume sources have been scaled up for better visualisation)

Although the theory presented above assumes a point source, we can only model the plumes with a finite source. Therefore, we compare the results from both methods by estimating the virtual point source $z_v$ of our plumes with finite sources. We model the plumes with square sources instead of a circular source because of the advantages of using a structured grid and the challenges associated with meshing in the case of a higher plume with a circular source inside a Cartesian domain. A structured grid provides a high level of quality and control for a mesh independence study, faster convergence because the elements are aligned to the flow, and lesser computational time due to less memory for storing mesh structure. Most importantly, if the source size is relatively small compared with the domain (as is the case for our set-up), the effect of source shape on the steady behaviour of an axisymmetric plume can be neglected because the cross-section quickly transitions to a circular cross-section, as observed by Marjanovic, Taub & Balachandar (Reference Marjanovic, Taub and Balachandar2017). In our simulations the cross-section becomes circular over heights of less than $3d$. Moreover, square sources in a Cartesian rectangular domain have been successfully used to model single axisymmetric plumes (e.g. Mokhtarzadeh-Dehghan et al. Reference Mokhtarzadeh-Dehghan, König and Robins2006; Yan Reference Yan2007; Marjanovic et al. Reference Marjanovic, Taub and Balachandar2017), and multiple axisymmetric plumes (e.g. Durrani et al. Reference Durrani, Cook, McGuirk and Kaye2011; Lou et al. Reference Lou, He, Jiang and Han2019).

The computational mesh is a structured grid with quadrilateral (two-dimensional) and hexahedral (three-dimensional) elements. A constant grid spacing $\varDelta$ in the $x$ and $y$ directions and a stretched grid clustered at the sources with a growth rate of $1.02$ (i.e. $2\,\%$ increase in neighbouring cells) in the $z$ direction was implemented to capture the high turbulent kinetic energy that exists off-axis and away from the source (Pham et al. Reference Pham, Plourde and Doan2007). Table 1 provides the number of grid points $N_x$, $N_y$ and $N_z$ in the $x$, $y$ and $z$ directions for three different meshes used for the mesh independence study, which will be presented later. The details of the mesh used for all simulations are given under mesh M1 in table 1, where $\varDelta =0.1d$ in the $x$ and $y$ directions was found to be sufficient.

Table 1. Number of grid points $N$ in $x, y$ and $z$ directions for three different meshes.

To closely model ‘interacting plumes’ in an unrestricted domain, we impose a Dirichlet boundary condition on the side, top and bottom boundaries, allowing for ambient fluid to enter at temperature $T_a= 300\ {\rm K}$ and leave the domain. This was implemented using a pressure inlet boundary condition on the sides and bottom boundaries and imposing a pressure outlet condition on the top boundary. The sources were modelled using a pressure inlet condition at $T_0 = 800\ {\rm K}$, where $T_0>T_a$. Transition to turbulence occurs at about Grashof number $Gr\approx 10^9$ (Bill & Gebhart Reference Bill and Gebhart1975), where $Gr$ is defined as

(4.13)\begin{equation} Gr = { \left ( \frac{\mu_a}{\rho _a} \right ) } ^2 g z^3 \frac{(T_0-T_a)}{T_a}, \end{equation}

and $\mu _a$ is the ambient dynamic viscosity.

Therefore, $T_0$ is chosen such that the plume is turbulent relatively quickly at a height $z/\chi \approx 1$ in the computational domain. The gauge pressure for all boundaries was set to zero. A no-slip condition was imposed on the side walls of the higher plume source with length $\zeta$ for representing vertically offset plumes. The horizontal separation $\chi$ is such that the lower plume is unaffected by the no-slip side wall of the higher plume i.e. $z<0$, consistent with the assumption made in our analytical model.

4.3. Numerical method

The transport equations are discretised using a finite-volume-based method. A second-order upwind scheme was used for all convective terms. The least square method was used to compute the scalar values at the cell faces, velocity derivatives and secondary diffusion. We used the SIMPLE-consistent (SIMPLEC) scheme to solve the pressure–velocity coupling. A multigrid algorithm that computes corrections on a series of coarse grid levels was implemented to accelerate convergence. At the start of each simulation run, the entire domain was initialised to a temperature $T=T_a + \delta {T}$ to improve convergence, where $\delta {T}$ is a small, uniform perturbation. The ratio $\delta {T}/T_a \approx 0.016$; however, the results were observed to be independent with respect to changes in $\delta T$. All simulations ran until the convergence criteria is satisfied, such that the residual of the conservation and transport equations is less than $10^{-5}$.

Since we only specify gauge pressure at the source and the source does not exhibit self-similarity, we calculate the dimensional source volume $Q_0$, momentum $M_0$ and buoyancy $F_0$ fluxes by summing the flux of each grid cell at the square source, where the fluxes for each grid cell are defined as

(4.14a–c)\begin{equation} Q_0= {\varDelta}^2 w, \quad M_0 = {\varDelta}^2w^{2}, \quad F_0 = {\varDelta}^2wg' , \end{equation}

and $\varDelta =0.1d$ is the constant grid spacing.

Assuming the plume is axisymmetric and self-similar, we determine the dimensional fluxes at different heights above the source by integrating as follows:

(4.15)\begin{equation} \left.\begin{gathered} Q(z) = 2{\rm \pi}\int_{0}^{ \tilde{b}(z)} \tilde{w}(r,z)r \, \text{d}r, \quad M(z) = 2{\rm \pi}\int_{0}^{ \tilde{b}(z)} {\tilde{w}(r,z)}^2r \, \text{d}r, \\ F(z) = 2{\rm \pi}\int_{0}^{ \tilde{b}(z)} \tilde{w}(r,z) \widetilde{g'}(r,z)r \, \text{d}r . \end{gathered}\right\} \end{equation}

Here $r$ is the radial coordinate relative to the local plume centre, and $\tilde {w}(r,z)$ and $\widetilde {g'}(r,z)$ correspond to the vertical velocity and buoyancy profiles, respectively.

The boundary conditions allow ambient fluid into the domain that leads to a component of vertical flow within the domain. Therefore, we only integrate within the plume by defining a plume radius $\tilde {b}(z)$ to accurately quantify the fluxes. According to Cook & Lomas (Reference Cook and Lomas1998), the Gaussian plume radius $\tilde {b}(z)$ is estimated as the radial distance from the plume axis to the point at which the vertical velocity $\tilde {w}$ has fallen to $1/e$ of its axial value. We evaluate the accuracy of this approximation by comparing the volume flux $Q$ scaling with theory (Morton et al. Reference Morton, Taylor and Turner1956) in § 5.2. In the case of two interacting plumes, the flux for each independent plume is summed since we only compare results in the region where the plumes can be referred to as independent, as described by Cenedese & Linden (Reference Cenedese and Linden2014). Finally, to compare numerical results with theory, we use the approximation proposed by Cook (Reference Cook1998) that links the Gaussian profile to the top-hat profile

(4.16a,b)\begin{equation} w = \tilde{w}_c/2, \quad b = \sqrt{2} \tilde{b}, \end{equation}

where $\tilde {w}_c$ is the vertical centreline (peak) velocity of the Gaussian profile, and $w$ and $b$ are the top-hat velocity and top-hat radius, respectively.