6.1 Numerical implementation of convection

As discussed in section 3.2, the convection is modeled by moving the vortices according to the equation

$$\frac{d}{dt} \vec{x}_i(t) = \sum_{j}\,\Gamma_j \, \vec{K}_\delta(\vec{x}_i(t);\vec{x}_j(t)) \quad \ ,$$ (6.1)

where $\Gamma_j$ is the circulation of vortex $j$, $\vec{x}_i$ is the position of the vortex $i$ at time $t$, and $\Gamma_j \vec{K}_\delta(\vec{x};\vec{x}_j)$ the velocity field of the vortex $j$ alone.

To integrate (6.1), we used the fourth order Runge-Kutta time stepping scheme proposed by Blum [27].

Unfortunately, the cost of computing the velocity of all the vortices using the summation (6.1) is proportional to $O(N^2)$, where $N$ is the number of vortices. Such computational effort would be unrealistic for flows where significant small-scale motion requires a fine vortex spacing, in other words, large $N$.

To solve this dilemma, `fast' algorithms were developed by, among others, Greengard and Rokhlin [97] and Carrier et al. [36], and independently by Van Dommelen and Rundensteiner [233,240]. For the computations in this work the latter of these three schemes was used; it seems to have been the first available scheme that was solution adaptive [240], but it can be noticeably slower than the first two schemes. For all these schemes, the required amount of work is roughly proportional to $N$. More recent variations have been proposed by a number of authors, such as [2,5,10,74,75].

As explained in section 3.1, the velocity field $\vec{K}_\delta$ of a vortex is obtained by convolving the velocity field $\vec{K}$ of a point vortex with a smoothing function $\phi_\delta(\vec{x})$ which is usually taken to be of the form

$$\phi_\delta(\vec{x}) \;=\;\frac{1}{\delta^2}\: \phi\left( \frac{\vec{x}}{\delta} \right) \quad \ ,$$ (6.2)

with $\phi(\vec{x})$ an axisymmetric function that integrates to unity. Hence, to obtain $\vec{K}_\delta$ in (6.1) a smoothing function $\phi(\vec{x})$ must be chosen. Based on procedures for the random walk computations [213,231,237], a low-order algebraic smoothing function was chosen of form

$$\phi(\vec{x}) = \frac{1}{\pi}\;\frac{1}{ \left(1+\vert\vec{x}\vert^2\right)^2} \quad \ .$$ (6.3)

This smoothing function falls in a class discussed by Hald [106] for the Euler equations. A relatively small core size $\delta = \sqrt{0.5\nu\Delta t}$ was used, as is common in practical applications; see Goodman, Hou, & Lowengrub [91].