The motion of buoyant point vortices

TL;DR

This paper extends the Hamiltonian framework to include buoyancy effects in point vortices, systematically analyzing their dynamics, stability, and chaotic behavior in various vortex configurations.

Contribution

It introduces a novel formulation for buoyant vortices, explores their interactions, and identifies conditions for bounded or unbounded motion and chaos.

Findings

Two-vortex systems can be bounded or unbounded with a clear criterion.

Maximum height achievable by buoyant vortex pairs is theoretically established.

Chaotic motion is demonstrated in the buoyant three-vortex system.

Abstract

A general formulation is presented for studying the motion of buoyant vortices in a homogeneous ambient fluid. It extends the well-known Hamiltonian framework for interacting homogeneous point vortices to include buoyancy effects acting on the vortices. This is then used to systematically examine the buoyant 1-, 2-, and 3-vortex problems. In doing so we find that 2 buoyant vortices may either evolve as a pair in bounded circular orbits, or as two independent unbounded vortices that drift apart, and a criteria is found to distinguish these cases. Special attention is given to the buoyant vortex couple, consisting of two vortices of equal and opposite circulation, and equal buoyancy anomaly. We show that a theoretical maximum height is generally possible for the rise (or fall) of such couples against buoyancy forces. Finally, the possibility and onset of chaotic motions in the buoyant 3-vortex problem is addressed. In contrast to the homogeneous 3-vortex problem, the buoyant vortex system shows evidence that chaos is present. We also demonstrate the chaotic advection of tracer parcels arising from the flow field induced by just 2 buoyant vortices.

Equations

$$\frac{\partial}{\partial t}\int_{V}\mathbf{u}(\textbf{\em x},t)\mathrm{d}V=\int_{V}\bigg{(}\mathbf{u}\times\bm{\omega}+\textbf{\em F}\bigg{)}\mathrm{d}V-\int_{\partial V}\bigg{(}\frac{p}{\rho_{a}}+\frac{1}{2}|\textbf{u}|^{2}\bigg{)}\hat{\mathbf{n}}\,\mathrm{d}A,$$

$$\text{emF}\equiv -\frac{g}{{\rho }_a}\rho \hat{\mathbf{j}}.$$

$$\frac{\mathrm{d}}{\mathrm{d}t}\mathbf{I}^{(v)}=\int_{V^{(v)}}\bigg{(}\mathbf{u}_{\mathrm{ext}}\times\bm{\omega}+\textbf{\em F}\bigg{)}\mathrm{d}V.$$

$${\mathbf{I}}^{(v)}=\frac{1}{2}{\int }_{V^{(v)}}\text{emx}\times \mathbf{\omega }\mathrm{d}V.$$

$${\mathbf{L}}^{(v)}=-\frac{1}{2}{\int }_{V^{(v)}}{\Vert \text{emx}\Vert }^2\mathbf{\omega }\mathrm{d}V,$$

$$\frac{\mathrm{d}}{\mathrm{d}t}\mathbf{L}^{(v)}=\int_{V^{(v)}}\textbf{\em x}\times\bigg{(}\mathbf{u}_{\mathrm{ext}}\times\bm{\omega}+\textbf{\em F}\bigg{)}\mathrm{d}V.$$

$${\mathbf{u}}_{\mathrm{ext}}(\text{emx})=\frac{1}{2(n-1)\pi }{\int }_V\frac{{\mathbf{\omega }}_{\mathrm{ext}}({\text{emx}}^{\prime })\times (\text{emx}-{\text{emx}}^{\prime })}{{\Vert \text{emx}-{\text{emx}}^{\prime }\Vert }^n}\mathrm{d}{V}^{\prime },$$

$${\mathbf{I}}^{(v)}=-{\int }_{A^{(v)}}\omega \text{emx}\times \hat{\mathbf{k}}\mathrm{d}A\mathrm{and}{\mathbf{L}}^{(v)}=\frac{1}{2}{\int }_{A^{(v)}}{\Vert \text{emx}\Vert }^2\omega \hat{\mathbf{k}}\mathrm{d}A.$$

$$\omega (x,y)=\sum _{i=1}^N{\Gamma }_i\delta (\text{emx}-{\text{emx}}_i).$$

$$\text{emF}(x,y)=\sum _{i=1}^N{\tilde{g}}_i\delta (\text{emx}-{\text{emx}}_i)\hat{\mathbf{j}}.$$

$${\tilde{g}}_i\equiv {\int }_{A_i^{(v)}}\text{emF}\cdot \hat{\mathbf{j}}\mathrm{d}A=-{\int }_{A_i^{(v)}}\frac{\rho (x,y)}{{\rho }_a}g\mathrm{d}A.$$

$$\mathcal{H}=-\sum _{i,j}^N\frac{{\Gamma }_i{\Gamma }_j}{4\pi }\ln [(x_i-x_j{)}^2+(y_i-y_j{)}^2]-\sum _{i=1}^N{\tilde{g}}_i{y}_i,$$

$${\Gamma }_i\frac{\mathrm{d}{x}_i}{\mathrm{d}t}=-\frac{\partial \mathcal{H}}{\partial y_i}\Rightarrow \frac{\mathrm{d}}{\mathrm{d}t}{I}_y^{(i)}=\sum _{i,j}^N\frac{{\Gamma }_i{\Gamma }_j}{2\pi }\frac{y_i-y_j}{(x_i-x_j{)}^2+(y_i-y_j{)}^2}+{\tilde{g}}_i,$$

$${\Gamma }_i\frac{\mathrm{d}{y}_i}{\mathrm{d}t}=\frac{\partial \mathcal{H}}{\partial x_i}\Rightarrow -\frac{\mathrm{d}}{\mathrm{d}t}{I}_x^{(i)}=-\sum _{i,j}^N\frac{{\Gamma }_i{\Gamma }_j}{2\pi }\frac{x_i-x_j}{(x_i-x_j{)}^2+(y_i-y_j{)}^2},$$

$$\sum _{i=1}^N\frac{\partial \mathcal{H}}{\partial x_i}=0\Rightarrow -\frac{\mathrm{d}}{\mathrm{d}t}{I}_x\equiv \sum _{i=1}^N{\Gamma }_i\frac{\mathrm{d}{y}_i}{\mathrm{d}t}=0,$$

$$-\sum _{i=1}^N\frac{\partial \mathcal{H}}{\partial y_i}=\sum _{i=1}^N{\tilde{g}}_i\Rightarrow \frac{\mathrm{d}}{\mathrm{d}t}{I}_y\equiv \sum _{i=1}^N{\Gamma }_i\frac{\mathrm{d}{x}_i}{\mathrm{d}t}=\sum _{i=1}^N{\tilde{g}}_i.$$

$$\sum _i\frac{\partial \mathcal{H}}{\partial {\theta }_i}=-\sum _{i=1}^N{\tilde{g}}_i{x}_i\Rightarrow \frac{\mathrm{d}}{\mathrm{d}t}{L}^2\equiv \frac{1}{2}\frac{\mathrm{d}}{\mathrm{d}t}\sum _{i=1}^N{\Gamma }_i(x_i^2+y_i^2)=-\sum _{i=1}^N{\tilde{g}}_i{x}_i.$$

$$\dot{x}=\tilde{g}/\Gamma \mathrm{and}\dot{y}=0,$$

$$\mu \ddot{x}+\Gamma \dot{y}=0\mathrm{and}\mu \ddot{y}-\Gamma \dot{x}+\tilde{g}=0,$$

$$\dot{x}(t)=C_1\cos (\alpha t)+C_2\sin (\alpha t)+\tilde{g}/\Gamma \mathrm{and}\dot{y}(t)=C_1\sin (\alpha t)-C_2\cos (\alpha t),$$

$$\Pi \equiv \frac{{\mu }^2{\tilde{g}}^4}{{\Gamma }^8},$$

$${\Pi }^{1/4}\sim \frac{g^{\prime }a}{V_{\max }^2},$$

$$\text{emr}\equiv (\zeta ,\eta )=(x_1-x_2,y_1-y_2),$$

$$\textbf{\em R}\equiv(X,Y)=\Big{(}\frac{\Gamma_{1}x_{1}+\Gamma_{2}x_{2}}{\Gamma},\frac{\Gamma_{1}y_{1}+\Gamma_{2}y_{2}}{\Gamma}\Big{)},$$

$${\text{emx}}_1=\text{emR}+\frac{{\Gamma }_2}{\Gamma }\text{emr}\mathrm{and}{\text{emx}}_2=\text{emR}-\frac{{\Gamma }_1}{\Gamma }\text{emr}.$$

$$\mathcal{H}=-\frac{{\Gamma }_1{\Gamma }_2}{4\pi }\ln ({\zeta }^2+{\eta }^2)-\frac{{\tilde{g}}_1{\Gamma }_2-{\tilde{g}}_2{\Gamma }_1}{\Gamma }\eta -\tilde{g}Y,$$

$$\frac{{\Gamma }_1{\Gamma }_2}{\Gamma }\dot{\zeta }=-\frac{\partial \mathcal{H}}{\partial \eta },\frac{{\Gamma }_1{\Gamma }_2}{\Gamma }\dot{\eta }=\frac{\partial \mathcal{H}}{\partial \zeta },$$

$$\Gamma \dot{X}=-\frac{\partial \mathcal{H}}{\partial Y},\Gamma \dot{Y}=\frac{\partial \mathcal{H}}{\partial X},$$

$$\frac{\mathrm{d}}{\mathrm{d}t}\text{emR}=(\tilde{g}/\Gamma ,0),$$

$$\mathcal{H}_{\ast}=-\frac{\Gamma_{1}\Gamma_{2}}{\Gamma^{2}}\Big{[}\frac{1}{4\pi}\mathrm{ln}(\zeta^{2}_{\ast}+\eta^{2}_{\ast})+D\eta_{\ast}\Big{]}-Y_{\ast},$$

$${\dot{\zeta }}_{\ast }=\frac{1}{2\pi }\frac{{\eta }_{\ast }}{{\zeta }_{\ast }^2+{\eta }_{\ast }^2}+D\mathrm{and}{\dot{\eta }}_{\ast }=-\frac{1}{2\pi }\frac{{\zeta }_{\ast }}{{\zeta }_{\ast }^2+{\eta }_{\ast }^2}.$$

$${\zeta }_{\ast }^2+{\eta }_{\ast }^2=C^2\exp (-4\pi D{\eta }_{\ast })$$

Full text

††preprint: APS/123-QED

I Introduction

Ever since the classical paper of Helmholtz, [16], vortex dynamics has been used as a tool to understand the behaviour of complex fluid flows. In the words of Aref, [1]: ‘…the evolution of vorticity, and thus the motions of vortices, are essential ingredients of virtually any real flow. Hence vortex dynamics is of profound practical importance’. Based on this principle, in the present paper we consider a general model for the interaction of singular ‘point’ vortices that are subject to buoyant forces in a homogeneous ambient fluid. We find that the presence of buoyancy leads to richer dynamics than in the homogeneous vortex problem, and displays some non-intuitive solutions, including the more rapid appearance of chaos compared to homogeneous flows as the number of vortices is increased. The present work generalises previous studies to account for any number of buoyant vortices, and in any initial configuration. This allows for a number of generalisations of the behaviour of buoyant vortex interactions. It also has the advantage of fitting into the well-studied Hamiltonian framework.

The first study that considered two interacting buoyant point vortices was that of Turner, [25]. Motivated by the release of effluent from chimneys, and buoyant ‘thermals’ found in convection flows, he investigated the motion of two vortices of opposite circulation and equal buoyancy anomaly relative to the ambient fluid. The motion of this buoyant vortex couple, was found by Turner, [25] for the special configuration of a purely vertical trajectory in the direction of the buoyancy force (i.e., for a buoyant couple rising, or dense couple falling). In this case he made the non-intuitive finding that increasing the buoyancy difference leads to a decrease in the rate of rise of the buoyant couple.

Since the original paper of Turner, [25], there have been a number of studies that have extended his analysis, nearly always focussing on the important case of two vortices with opposite circulation [e.g., 20], and often in more complicated ambient environments with shear and/or stratification [e.g., 21, 15]. In a series of papers by Arendt [6, 7, 8] the motion of buoyant vortices in polytropic stratified fluids, where the pressure and density are related through a power law, have been systematically explored. Many of the results in Arendt, [8] will be found to have analogues in the present investigation where we focus on buoyant vortices in a homogeneous ambient. Recently, Ravichandran et al., [20] examined the collision and collapse of such buoyant vortex couples, and conducted simulations using many vortex patches to demonstrate the presence of this mechanism in a random field of vortices.

Similar to the general framework established by Arendt, [8], we formulate the motion of many buoyant point vortices of arbitrary buoyancy anomaly and circulation in a homogeneous fluid, into a single general framework. Our motivation for such a formulation arises from the study of instability development in stratified shear layers, which has been found to lead to the formation of concentrated patches of buoyancy and vorticity [e.g., 11, 23]. An idealised model for the interaction of these buoyant vortex structures, such as formulated here, is therefore of interest in order to understand the evolution of stratified shear layers, with implications for oceanic and atmospheric mixing [14].

The paper is organised as follows. In the next section we derive a series of conservation equations for the motion of buoyant vortices, and demonstrate that they are captured by a Hamiltonian framework. Given this general formulation of the laws governing the motion of buoyant point vortices, we then go on to explore the solutions by increasing the number of vortices from $N=1$, to $N=3$, in sections III to V. Special cases of interest, such as the buoyant vortex couple, and an exploration of chaotic motion in the buoyant 3-vortex problem, will be examined along the way. A summary and conclusions follows in the final section.

II General formulation

II.1 Conservation laws

In this subsection II.1, we present a brief background of the formulation of vortex conservation laws for Boussinesq flows [for details about homogeneous flows see 22].

Consider an inviscid, non-diffusive, Boussinesq fluid inside a simply connected domain of fixed volume $V$ (which could be infinite). The background state, with density $\rho_{a}$, is assumed to be irrotational and in hydrostatic balance. Denoting perturbation velocity, vorticity, pressure and density respectively by $\mathbf{u}$, $\bm{\omega}(\equiv\nabla\times\textbf{u})$, $p$ and $\rho$, the perturbed Boussinesq Euler equations on integration yield

[TABLE]

where

[TABLE]

Here $g$ denotes gravity, $\hat{\mathbf{j}}$ the unit vector in the vertical $y$ direction, and $\hat{\mathbf{n}}$ the unit vector normal to the bounding surface $\partial V$. We will focus on the dynamics of an isolated vortex patch (i.e., vorticity distribution with a compact support) of volume $V^{(v)}\in V$. The vortex patch may have a density different from $\rho_{a}$, hence yielding the body force (or buoyancy) term F. The vector quantity $\mathbf{u}\times\bm{\omega}$ is referred to by different names by different authors, such as the vortex force [22], Magnus force [20], or the lift force [19], however, we will use the latter term in this study.

By taking $\textbf{\em x}\times(\nabla\times)$ of the integrands in Eq.  (1) and making use of various vector calculus identities, detailed in Saffman, [22], we arrive at

[TABLE]

This expresses the conservation of linear impulse of the vortex patch, defined as

[TABLE]

As discussed in [22] and [13], this result relies on the fact that both $\mathbf{u}$ and $p$ due to an isolated vortex decay as $\mathcal{O}(r^{-3})$, where $r$ is the distance from the vortex. Therefore, the surface integrals which appear in this process vanish. In Eq.  (2), the quantity $\mathbf{u}_{\mathrm{ext}}$ represents the induced velocity within $V^{(v)}$ due to all vorticity external to it. Similarly, defining the angular impulse of an isolated vortex as

[TABLE]

we obtain the conservation of angular impulse

[TABLE]

If multiple such vortex patches are present (and are far away from the domain boundary) in the irrotational background flow field, $\mathbf{u}_{\mathrm{ext}}$ on a given patch would arise from the action of the other vortex patches present in the system. It can be obtained using the Biot-Savart law, which inverts the vorticity field $\bm{\omega}_{\mathrm{ext}}$ existing due to the other patches:

[TABLE]

where the number of spatial dimensions is $n=2$ or $3$. From Eq.  (2) and Eq.  (3), we recover the well-known result that in the absence of any external velocity field and applied force, both linear and angular impulses of a homogeneous vortex patch (i.e., with the same density as the background) are constants of motion.

In this paper, we will be focusing on two-dimensional ($x$–$y$) systems. The vorticity vector is then given by $\bm{\omega}=-\omega\mathbf{\hat{k}}$, where $\omega(x,y)=\partial u/\partial y-\partial v/\partial x$, and $(u,v)$ are fluid velocities in the $(x,y)$ directions. The linear and angular impulses of an isolated vortex patch of area $A^{(v)}$ are given by [22]

[TABLE]

The missing $1/2$ in $\mathbf{I}$ for two-dimensional flows is a consequence of the vortex lines not being closed. When working with two-dimensional flows we also consider area integrals to be expressed per unit length.

The conservation of linear and angular impulse essentially represent rewritten forms for the linear and angular momentum that are more practical for vortex dynamics studies. This will become apparent in the next subsection when we use them to formulate point vortex motions.

II.2 Equations governing buoyant point vortex dynamics

Consider $N$ vortices, denoted by the subscript $i=1,...,N$, that move across the unbounded $x$–$y$ plane in time, $t$. We assume that the diameter of the vortices is much smaller than the distance between them, and so we can take them to be represented by ‘point vortices’ at the locations $\textbf{\em x}_{i}=(x_{i},y_{i})$. The vorticity field is zero everywhere except at the locations of the vortices, where it has delta-function behaviour, i.e.,

[TABLE]

Here $\Gamma_{i}$ indicates the circulation strength of vortex $i$, with a $\Gamma_{i}>0$ indicating clockwise circulation.

In addition, we allow each vortex to be composed of fluid with a density anomaly from the constant ambient density $\rho_{a}$. Each vortex will then experience a buoyant point force, and the field of these forces can be written as

[TABLE]

F is expressed as a force per unit mass (per unit length in two dimensions), so that the ‘reduced gravity’ of each vortex, $\tilde{g}_{i}$, is defined by the difference in density with the ambient $\rho(x,y)$, and the vortex cross sectional area $A_{i}^{(v)}$, as

[TABLE]

Hence $\tilde{g}_{i}>0$ implies that the vortex is (positively) buoyant. This definition of the reduced gravity, $\tilde{g}_{i}$, differs from the usual definition by the area factor. Since we are assuming the fluid to be Boussinesq, the inertia of the vortices resulting from their difference in density from the ambient fluid is neglected. Note that in taking our vortex to be concentrated at a point and shrinking the area to zero, we ensure that the net buoyant force on the vortex patch is fixed. This limit might appear to violate the Boussinesq approximation, as the density difference with the ambient should become infinite, however, it must be kept in mind that these limits are only mathematical tools used when making a point approximation to a finite-sized vortex. Another option would be to take $g\to\infty$, as in Aref and Tryggvason, [3], which would avoid this issue. We observe that the integral terms in Eqs. (2)–(3) preserve the conservation of linear and angular impulse in this point vortex representation [e.g., 25].

As in the case of a system of homogeneous vortices [17], in the system just described it is possible to use a Hamiltonian framework to find the vortex trajectories. The Hamiltonian $\mathcal{H}$, is given by the energy

[TABLE]

being composed of both a ‘kinetic energy of interaction’ (given by the first term on the right hand side), and a potential energy (last term), and is a constant of motion 111In this Hamiltonian system the generalised coordinates, and momenta take the form of $q_{i}=\Gamma_{i}y_{i}$, and $p_{i}=x_{i}$, respectively, where $dp_{i}/dt=-\partial\mathcal{H}/\partial q_{i}$ and $dq_{i}/dt=\partial\mathcal{H}/\partial p_{i}$.. Note that the sum in the kinetic energy term is not carried out over the singular terms $i=j$. Substituting Eqs.  (6)–(7) in the linear impulse conservation equation Eq.  (2), and using the definitions in Eqs.  (4, 5), we obtain the following Hamiltonian dynamical system of $2N$ dimensions:

[TABLE]

where $\mathbf{I}^{(i)}=\left(I_{x}^{(i)},I_{y}^{(i)}\right)$ represents the linear impulse of vortex $i$. It is important to note that the summation terms in Eq.  (10) arise from the externally induced lift force $\int_{V^{(v)}}\mathbf{u}_{\mathrm{ext}}\times\bm{\omega}\mathrm{d}V$ of the linear impulse conservation Eq.  (2).

It is also possible to find global conservation laws for all buoyant vortices in a straight forward manner directly from the Hamiltonian description. In this regard we follow the analyses of [18], who deduced the conservation laws for homogeneous point vortices,

[TABLE]

Also, if we transform into polar coordinates, with $(x_{i},y_{i})=(r_{i}\cos\theta_{i},r_{i}\sin\theta_{i})$, it is possible to derive the conservation of the global angular impulse:

[TABLE]

This shows that when a global budget (including all vortices) is considered, the lift forces can be regarded as internal forces, and will cancel to produce no net effect.

The Hamiltonian in Eq. (9) has correctly reproduced the conservation laws of linear and angular impulse, derived in a general form in the previous subsection, for a buoyant point vortex system. This can be verified by substitution of Eqs. (6)–(7) into Eqs. (2)–(3). These conservation laws, and their alterations under the force of buoyancy from the homogeneous vortex case, will be critical for determining the nature of the solutions in what follows. Note that the Hamiltonian itself is also a conserved quantity, corresponding to energy conservation. These conservation laws and Hamiltonian structure also have an analogue in the configuration studied by Arendt, [8].

III A single buoyant vortex

For the case of a single buoyant vortex (whose parameters are not denoted with a subscript) the Hamiltonian equations of motion in Eqs.  (10a)–(10b) become

[TABLE]

where the dot represents ordinary differentiation with respect to $t$. This shows that the vortex propagates in the horizontal direction with speed $U=\tilde{g}/\Gamma$, and merely expresses the conservation of linear impulse for a single vortex. Note that for a vanishing buoyancy, this solution reduces to the homogeneous result of a stationary vortex. The non-intuitive result that a strictly vertical buoyancy force results in the horizontal motion of the vortex can be understood as follows. The Hamiltonian formulation reveals that total energy $\mathcal{H}$ is a constant of motion. Since there is no kinetic energy of interaction between vortices, the total potential energy must remain constant, giving $\dot{y}=0$. However, in order to keep the vortex at the same vertical level, there must be an additional force that is needed to balance the vertical buoyancy force. As long as the vortex is moving horizontally, this force is provided by the aerodynamic lift (or vortex force) that is present on any body possessing a circulation in a cross flow [19]. In terms of force per unit length, the lift force is $\rho_{a}\Gamma U$, and the buoyancy force is $\rho_{a}\tilde{g}$, which leads directly to the resulting translation velocity. The appearance of an equivalent horizontal drift of a vortex once subject to a body force is also derived by Saffman, [22], and discussed by Arendt, 1993b [7]. In Appendix A, we numerically show that a buoyant point vortex (modeled as a Gaussian vortex patch) indeed travels horizontally with a constant speed.

One may also be tempted to infer that the horizontal motion of the buoyant vortex arises due to the absence of an inertial term (from making the Boussinesq approximation) that can lead to a vertical acceleration of the vortex under the action of gravity. This can easily be added in an ad hoc manner, and the equations of motion written as

[TABLE]

with $\mu$ a measure of the mass of the vortex. A similar balance was considered in Ravichandran et al., [20]. However, their model is based on potential flow around a circular cylinder and it is not clear how it may be constructed from the equations of motion applied to vortex motions (i.e. proceeding similar to that outlined in section II). Proceeding nonetheless, a general solution for the velocities $(\dot{x},\dot{y})$ of this system can be found as

[TABLE]

where $C_{1},C_{2}$ are constants that depend on the initial conditions. This solution describes a horizontal drifting motion of the vortex at the speed $\tilde{g}/\Gamma$ and an oscillation at the frequency $\alpha\equiv\Gamma/\mu$.

The vortex motion is comprised of a balance between inertia and the buoyancy and lift forces. Note that if inertia is neglected, the lift-buoyancy balance leads to the constant horizontal drift of $\dot{x}=\tilde{g}/\Gamma$. This arises from the fact that the lift force is always perpendicular to the vortex motion, and is the only force that can balance the vertical buoyancy force. The oscillation arises due to the exchange between the kinetic and potential energy components of the total energy, which must be conserved. The addition of inertia to the system, therefore does not significantly change the fundamental force balance, or the vortex trajectory. This neglect of the vortex inertia will be appropriate as long as the dimensionless number

[TABLE]

which quantifies the relative importance of the inertial term, is small. If we approximate the mass of the vortex by its area, $A^{(v)}$, assume that the vortex is circular with radius $a$, and approximate the circulation by a solid body rotation with maximum velocity $V_{\mathrm{max}}$, we can write

[TABLE]

with $g^{\prime}\equiv\Delta\rho g/\rho_{a}$ the standard reduced gravity. This shows that $\Pi^{1/4}$ has the form of a Richardson number, or $\Pi^{-1/8}$ a Froude number.

IV Two interacting buoyant vortices

IV.1 Solution

Using the Hamiltonian system Eqs.  (10a)–(10b) with $N=2$ vortices, now leads to a non-zero kinetic energy of interaction term, and considerably richer dynamics. In analogy with the two-body problem of classical physics, it is possible to make considerable analytical progress in solving this system by introducing new variables given by the vortex separation

[TABLE]

and the centre of vorticity

[TABLE]

as illustrated in Fig 1.

Here we have defined $\Gamma\equiv\Gamma_{1}+\Gamma_{2}$, and the original coordinates can be recovered by

[TABLE]

With the new coordinates the Hamiltonian takes the form

[TABLE]

where $\tilde{g}\equiv\tilde{g}_{1}+\tilde{g}_{2}$. The original $2N$-dimensional Hamiltonian system has been reduced to a $(2N-2)$-dimensional Hamiltonian system, and we can rephrase the whole problem as the following set of equations:

[TABLE]

for R. Notice that in this set of equations the motion of the centre of vorticity R is decoupled from the vortex separation r. In particular, the centre of vorticity moves with the constant horizontal drift as if it were a single vortex with circulation $\Gamma$, and reduced gravity $\tilde{g}$. This may be seen from Eq.  (22b), i.e.,

[TABLE]

giving the same result as Eq.  (13), which can easily be integrated. Then we are left with a system of two first-order ODEs for the components of r.

At this point it is useful to non-dimensionalize the equations using $\Gamma^{2}/\tilde{g}$ as a length scale, and $\tilde{g}/\Gamma$ as a velocity scale. The corresponding form for the Hamiltonian becomes

[TABLE]

where all variables with an asterisk are non-dimensional, and we have defined the dimensionless number $D\equiv(\Gamma/\tilde{g})(\tilde{g}_{1}/\Gamma_{1}-\tilde{g}_{2}/\Gamma_{2})$. The system of equations describing the trajectories is then written as

[TABLE]

Note that for $D=0$ the equations describing two homogeneous point vortices are recovered. We can gain considerable insight into the solution for the vortex trajectories by using the fact that the total energy is conserved. Since $\mathcal{H}$ is a constant of the motion, given by its initial value $\mathcal{H}(t=0)\equiv E$, it allows us to write

[TABLE]

with

[TABLE]

a dimensionless constant given by the initial conditions. Note that from Eq.  (23) we can always define our coordinate system so that $Y=0$ for all $t$ (this is just a choice of datum for the potential energy).

Equation (26) defines a set of curves that describe the vortex separation distance over the entire evolution of the vortex interaction, and are shown in Fig 2(a).

Similar to the results of Arendt, [8] for the case of a stratified polytropic fluid, we can generally split the solution into two different categories: (i) those that remain bounded (i.e., within a finite distance) to the centre of vorticity, and (ii) those that are unbounded. For a given value of $D$, there is a critical value of the dimensionless parameter $C$, given by $C_{\mathrm{cr}}=(2\pi eD)^{-1}$ and shown in Fig 2(b) (which can be thought of as a measure of the initial distance between the vortices, or equivalently as the initial energy), above which there are unbounded orbits (Fig 2, white area). The separatrix, represented by the thick dark contour in Fig 2(a), separates the bounded and unbounded regions. This separatrix appears as a homoclinic orbit, with its hyperbolic fixed point located at

[TABLE]

and has eigenvalues $\pm 2\pi D^{2}$ and eigenvectors $(\mp 1,1)$. Hence the eigenvalues and eigenvectors are independent of the sign of $D$, leading to a clockwise circulation (in the bounded region) for both $D>0$ and $D<0$. The far-field flow, however, depends on the sign of $D$ as per Eq.  (25). We also note in passing that $(0,0)$ is not a fixed point, but on the contrary it corresponds to a point of infinite relative velocities, as can be seen from the coloured arrows in Fig 2(a). This singularity is not realistic as we have assumed that the separation of the vortices is much larger than their diameter.

If the initial energy of the vortices is large enough, i.e., $C$ is large enough, then unbounded orbits will result, and this is more likely to occur for large $D$. Note that for the singular case of homogeneous point vortices $D=0$, and we recover the result that all vortex trajectories are bounded. The dependence of the orbits on the parameter $D$ can be interpreted more easily if we write it in terms of the individual vortex velocities assuming there is no interaction between them, i.e., $U_{i}\equiv\tilde{g}_{i}/\Gamma_{i}$. Then we have

[TABLE]

where $U_{\mathrm{cv}}\equiv\tilde{g}/\Gamma$ is the velocity of the centre of vorticity. As we would expect intuitively, when the difference in the individual vortex velocities is large compared to that of the centre of vorticity, they are able to ‘escape’ from a bounded orbit.

IV.2 Special cases

IV.2.1 Bounded and unbounded orbits - an example

It is helpful to examine an example of the transition between bounded and unbounded vortex trajectories. Fixing $D=(4\pi)^{-1}$ and choosing $C=2e^{-1}\pm\epsilon$ will produce a bounded trajectory for $C<C_{\mathrm{cr}}=2e^{-1}$, and an unbounded trajectory for $C>C_{\mathrm{cr}}$. In practise, we choose $\Gamma_{1}=2\Gamma_{2}$, and can compute the trajectories relative to R.

The results are shown in Fig 3. In the unbounded case, the vortices approach a constant horizontal trajectory in opposite directions, as the distance between them increases, and thus the interaction weakens. Recall that for a vanishing interaction each vortex will approach a horizontal trajectory at the speed $\tilde{g}_{i}/\Gamma_{i}$. Only an extremely small change in the initial energy (and therefore position) of the vortices, expressed by $C$, alters the orbits of the two from bounded to unbounded (Fig 3b). This diverging of the orbits occurs as the two are farthest from one another, and the interaction is weakest. This sensitive dependence on the initial condition in the two vortex problem is only present at this point in parameter space close to the critical separation between bounded and unbounded orbits.

IV.2.2 The vortex couple

In the original analysis of Turner, [25], he considered the special case of the strictly vertical propagation of a buoyant vortex couple, i.e., vortices with $\Gamma_{1}=-\Gamma_{2}$. This vortex configuration often arises in idealized models of convection ‘thermals’, and in starting jets. In these applications, and in Turner’s analysis, it is assumed that the couple is formed from an initial buoyant cloud, or jet, and therefore it is also the case that $\tilde{g}_{1}=\tilde{g}_{2}\equiv\tilde{g}_{0}$. In this subsection we extend Turner’s analysis and examine the motion of buoyant vortex couples in general.

To solve for the trajectory of the vortex couple we cannot use the formulation in the previous section since $\Gamma=0$; instead we will use the conservation laws directly. The conservation of linear impulse Eq.  (10) reduces to

[TABLE]

where $\Gamma_{0}\equiv\Gamma_{1}=-\Gamma_{2}$. This equation states that in order for a vortex couple to conserve linear impulse the vertical distance between vortices remains constant in time, while the horizontal distance changes linearly with time. The total energy for the couple can be written

[TABLE]

where we have defined the mean vertical position of the couple as $\bar{y}\equiv(y_{1}+y_{2})/2$. Differentiating this equation with respect to $t$, and combining it with Eq.  (29) allows us to solve for the evolution of the couple.

For large times, we can effectively ignore the initial separation of the vortices and, by Eq.  (29), write

[TABLE]

This can be combined with Eq.  (30) to find the long-time evolution of the elevation of the couple as,

[TABLE]

This relation demonstrates a number of noteworthy results. First, buoyant couples will rise, and dense couples will fall, inline with our intuition. Second, they will do so at a decreasing rate, proportional to $t^{-1}$. Last, it is interesting that as the magnitude of the buoyant force (measured by $\tilde{g}_{0}$) is increased, the rate of elevation change decreases. This is a generalisation of Turner’s result that less buoyant thermals will rise faster. The reason for this is clearly that the rate of separation increases with increasing $|\tilde{g}_{0}|$, due to the dominant balance for large times between the lift and buoyancy forces on the individual vortices. This results from increasing horizontal motion of the vortices required to increase the lift in response to increased buoyancy.

For short times, however, the couples can propagate in a direction that is opposite to the buoyant forces. An example in Fig 4 shows such an orientation and trajectory. In an arbitrary initial orientation it can be shown that the rate of elevation change of the couple is

[TABLE]

where $|\textbf{\em r}|^{2}=\zeta(t)^{2}+\eta_{0}^{2}$ with $\eta_{0}$ the initial vertical separation (which is constant in time). This shows that it is the vertical advection speed of each vortex on the other that is responsible for the mean rate of elevation change (but with the separation governed by the conservation of linear impulse). Using the conservation of $\mathcal{H}$ allows us to derive a maximum vertical displacement of the couple against buoyancy forces of

[TABLE]

where $|\textbf{\em r}_{0}|$ is the initial separation of the couple. This can alternatively be expressed in terms of the angle, $\theta$, that a line connecting the two vortices makes with the horizontal,

[TABLE]

This relationship for the dimensionless elevation change, $\Delta\bar{y}\tilde{g}_{0}/\Gamma_{0}^{2}$, is plotted in Fig.  4(b).

Note that the special case treated by Turner (with $\eta_{0}=0$ and $\theta=0$) can result in a singularity in finite time when the couple is directed opposite to the buoyant force, as it results in the collision of the vortices and an infinite vertical velocity of the couple. This result is clearly unrealistic as our assumption that the diameter of the vortices is much smaller than the distance between them breaks down. This case has been examined by Ravichandran et al., [20] who found that the collision results in significant deformation of the vorticity field into complicated structures, leading to the collapse of the individual vortex cores. The process of vortex couple collapse is also likely to happen for small, but non-zero $\theta$, as the distance between the vortices decreases during the elevation change of the couple (Fig.  4a). This will likely lead to a decrease in the elevation change predicted by the curve in Fig.  4(b).

IV.2.3 A ‘pseudo-homogeneous’ case

In this case, we note that trajectories of the buoyant vortices can occur as if they were homogeneous when $D=0$. This occurs when $\tilde{g}_{1}\Gamma_{2}=\tilde{g}_{2}\Gamma_{1}$, or similarly if the individual vortex speeds are equal, i.e., $U_{1}=U_{2}$. The solution for $|\textbf{\em r}|$ is

[TABLE]

which corresponds exactly to the equations for the homogeneous two-vortex problem with $C$ the initial distance between the vortices. The only difference between the trajectories of these buoyant vortices and two similar homogeneous vortices is the constant motion of the centre of vorticity at the speed $U_{\mathrm{cv}}=\tilde{g}/\Gamma$. The vortices are also rotating about the centre of vorticity at the frequency of $\Gamma/2\pi C^{2}$, independent of buoyancy. Rather than plotting the trajectories in the $xy$-plane, it is also possible to get an accurate qualitative feel for the trajectories by using the $\zeta\eta$-plane. This can be seen by writing

[TABLE]

so that the individual vortex trajectories, in a frame of reference moving with the centre of vorticity, are just a scaled version of the curve in the $\zeta\eta$-plane.

V Three buoyant vortices

A general property of Hamiltonian systems is the appearance of chaotic motions as the number of degrees of freedom is increased [24]. In the case of homogeneous point vortices, it is well known that chaotic motions appear once the number of vortices (equal to the number of degrees of freedom of the Hamiltonian system) $N\geqslant 4$ [2, 1]. The key to integrating the equations, and thereby eliminating the possibility of chaotic trajectories, is to determine $N$ independent integrals (constants) of motion, $F_{1},...,F_{N}$, where $F_{j}=\mathrm{constant}$ for all vortex positions and all $t$. These integrals of the motion are usually formed from the conservation laws that we found in Eq.  (10) and Eq.  (12), as well as the Hamiltonian (energy) itself. For example, in the homogeneous vortex case, three integrals of motion can be identified as $\mathcal{H}$, $I_{x}^{2}+I_{y}^{2}$, and $L^{2}$ [2]. These integrals can be used to effectively reduce the number of degrees of freedom so that the system may be solved (in theory) by an integration 222In Hamiltonian systems, each integral allows the elimination of a pair of conjugate variables (e.g., $\Gamma_{i}y_{i},\,x_{i}$ in our case). Thus if we have $N$ integrals of motion for a system of $2N$ Hamiltonian equations (i.e. $N$ degrees for freedom, and $N$ vortices) we can, therefore, reduce the problem to quadratures. Such systems are said to be Liouville integrable.. Although techniques exist to derive integrals of motion, such as Lie algebra and symmetries of the Legrangian [e.g., 12, 10, 9], one often relies on intuition [see 24, p. 39].

For the global conservation laws of the buoyant vortex system derived in Eqs.  (11)-(12) we can see that the addition of buoyancy forces has resulted in non-constant angular impulse $L$, and vertical linear impulse, $I_{y}$. We are therefore left with only two conserved quantities $\mathcal{H}$, the energy, and $I_{x}$, the horizontal linear impulse. Hence, the buoyant $3$-vortex problem appears to be non-integrable, as was also the case found for the stratified polytropic fluid system of Arendt, [8]. The purpose of this section is therefore to answer the question: Can three interacting buoyant vortices exhibit chaotic motion?

To address this question, we begin by showing in Fig.  5 an example of the evolution of three buoyant vortices, and comparing to the homogeneous case (which is known to be non-chaotic). In both cases, the circulation of the vortices is identical with $\Gamma_{i}=\{1,2,3\}$, but the buoyant 3-vortex system has $\tilde{g}_{i}=\{-0.2,0,0.2\}$ (Fig.  5b,d,f). Note that since $\sum_{i}\tilde{g}_{i}=0$, by Eq.  (11b), there is no net movement of the centre of vorticity. It can be seen from Fig.  5 that, whereas the homogeneous vortices trace out regular symmetric patterns, the buoyant vortices have irregular asymmetric trajectories. This is an indication of the presence of chaos in the buoyant 3-vortex problem, as we might expect from the loss of an independent integral of motion.

To further address the presence of chaos in the buoyant 3-vortex problem, we have numerically calculated the Lyapunov exponents. A defining feature of chaos, is that nearby initial conditions will produce diverging solutions after finite times. This feature is quantified through the maximum Lyapunov exponent defined as

[TABLE]

It characterises the exponential rate of separation of infinitesimally close trajectories whose initial separation is $\delta\mathbf{X}_{0}$. An autonomous non-linear dynamical system with $\lambda_{\mathrm{max}}>0$ is non-integrable, implying chaos is a possibility. We have numerically computed the Lyapunov exponents up to $t=100,000$ using the procedure outlined in Wolf et al., [26]. In Fig.  6, $\lambda_{\mathrm{max}}$ is shown for both three homogeneous vortices (black curve) and three buoyant vortices (red curve), corresponding to the case examined in Fig.  5. In the absence of reduced gravity, i.e., for the homogeneous 3-vortex case, $\lambda_{\mathrm{max}}\rightarrow 0$ with increasing $t$, as expected. However, the presence of buoyancy in the problem leads to a positive $\lambda_{\mathrm{max}}$, indicating that two nearby trajectories will exponentially diverge, further supporting the finding that chaos is present.

A point worth mentioning in the context of differences between the homogeneous and buoyant 3-vortex problems is the concept of chaotic advection – simple time-dependent flows can cause chaotic motion of tracer parcels. The basic idea stems from the fact that the motion of a point vortex with vanishing circulation and zero buoyancy anomaly will behave as a tracer particle, being advected by the fields induced by the other vortices present. When such a ‘tracer vortex’ is included in a field of three homogeneous point vortices, although the unsteady flow field produced by the three point vortices is integrable, the motion of the tracer vortex is not [5]. This result echos a similar finding for the 3-body problem with one of the point masses being negligible, and is commonly called the restricted 3-body problem. This concept can be easily extended to understand chaotic advection due to buoyant vortices, with an example of the restricted buoyant 3-vortex problem shown in Fig.  7. The flow field produced by two buoyant vortices is integrable (Fig.  7a) with repeated closed orbits, however, the motion of neutrally buoyant tracers in this flow field is not. The chaotic advection of the tracer particle is demonstrated in Fig.  7(b) to consist of irregular non-repeating orbits indicative of chaos.

VI Summary and conclusions

In this paper we have formulated a general Hamiltonian framework for the investigation of buoyant point vortices in a homogeneous ambient fluid. This extends previous work in this area that looked only at special cases of buoyant vortices of relevance for particular applications. It reveals analogous results to the general Hamiltonian framework developed for vortices in a stratified polytropic ambient fluid by Arendt, [8]. We have provided a systematic study of the 1- and 2-vortex problems, and discovered that the evolution of two buoyant vortices can be split into two different solution types: those in which the vortices remain bounded (i.e., within a finite distance) to each other, and those in which the vortices drift apart in time. The boundary between these two cases was also derived. As the number of vortices is increased to 3, the solutions become complex, and irregular chaotic motions result. This feature of the buoyant 3-vortex problem is in contrast to the homogeneous 3-vortex problem, and arises due to the loss of symmetry in the Hamiltonian due to the buoyancy force, thus reducing the number of integral invariants. We also investigate the possibility of chaotic advection of tracer parcels arising from the simple, time dependent flow field induced by two buoyant vortices. Future work could include a deeper investigation of chaotic motion in the 3-buoyant vortex problem, including a mapping of chaotic regions throughout parameter space, as well as investigating the statistical mechanics of many buoyant vortices.

Acknowledgements.

We would like to acknowledge funding from the Helmholtz Foundation under the PACES II programme, and the Alexander von Humboldt Association for Humboldt research fellowship to A. Guha carried out at the Helmholtz-Zentrum Geesthacht.

Appendix A Numerical validation of the constant horizontal drift of a buoyant vortex patch

We numerically solve the non-dimensional perturbation vorticity and scalar (density) conservation equations in two dimensions ($x$-$y$):

[TABLE]

where $D/Dt$ denotes the material derivative, $\mathrm{Re}$, $\mathrm{Fr}$ and $\mathrm{Pr}$ respectively denote the Reynolds, Froude and Prandtl numbers. The above set of equations are solved using a pseudospectral code in a $2\pi\times 2\pi$ domain with doubly periodic boundary conditions. The number of mesh-points used are $(N_{x},N_{y})=(512,512)$; the Crank–Nicolson method is used for time integration. The following parameters are assumed: $\mathrm{Re}=10^{3}$, $\mathrm{Fr}=10$ and $\mathrm{Pr}=1$. We simulate a generalized Gaussian vortex centered at ($\pi,\pi$) having the following vorticity and density structure at $t=0$:

[TABLE]

where $\beta=6$. While $\beta=1$ results in a standard Gaussian distribution, higher values of $\beta$ allow a flatter peak, resulting in a vortex patch of nearly constant vorticity and density (i.e. approximating a ‘buoyant’ Rankine vortex). We show the initial and final states of the vortex patch in Figs.  8(a) and 8(b), respectively. In Fig.  8(b), we find that the vortex has diffused, as expected in a viscous, diffusive flow. Nevertheless, the vortex indeed travels along a horizontal trajectory with a constant speed (Fig.  8c). This speed can be derived based on Eq. (13) for a finite sized buoyant vortex patch giving $U=g^{\prime}_{0}/\Omega_{0}$, where $\Omega_{0}$ is the near-constant initial vorticity of the patch and $g^{\prime}_{0}$ the buoyancy anomaly (standard reduced gravity). Non-dimensionalising this speed gives $U_{\ast}=\mathrm{Fr}^{-2}=0.01$, corresponding to our chosen numerical setup. This provides an excellent match to the results of the translation speed observed in the simulation of $U_{\ast}=0.0099$.