4.1 Mathematical formulation

Figure 4.1: Redistribution of the circulation of a vortex $\Gamma ^n_i$.

The purpose of the vorticity redistribution method is to simulate the diffusion of each vortex during a time-step. As sketched in figure 4.1, this is done by distributing fractions of the circulation $\Gamma_i^n$ of each vortex $i$ to its neighboring vortices. The question is how to select the neighboring vortices and the fractions so that the correct diffusion is approximated. We will answer that question in the following discussion.

First, vortices will be considered to be within the neighborhood of a given vortex $i$ if they are within a predetermined distance from that vortex. We take this distance to be of the order of the typical diffusion distance during a time-step. To be precise, the typical diffusion distance $h_v$ will be defined as

$$h_v \equiv \sqrt{\nu \Delta t}\ ,$$ (4.1)

where $\Delta t$ is the size of the time step and $\nu$ the coefficient of kinematic viscosity. A vortex $j$ is part of the neighborhood of vortex $i$ if

$$\vert\vec x_j - \vec x_i\vert \le R\,h_v \quad \ ,$$ (4.2)

where $R$ is a chosen constant. We used a maximum distance $\sqrt{12} h_v$ in all computations presented in this thesis. Some guidelines for choosing this distance will be given in subsection 6.2.3.

The diffusion of a vortex $i$ will be approximated by moving fractions of its circulation towards the other vortices within this neighborhood. We will indicate the fraction moved from vortex $i$ to a vortex $j$ by $f^n_{ij}$, where $n$ indicates the time level. Implementation of the redistribution method is in principle merely a matter of determining fractions $f^n_{ij}$ that approximate the correct diffusion over a time-step accurately and stably.

Yet, we choose not to identify the vortex strengths with any particular smooth interpolated vorticity distribution. The reason is that due to straining effects, the vortex locations can become very irregular. In the absence of a continuous vorticity field, the question arises how a meaningful representation of the diffusion process can still be achieved. Regardless of the interpretation of that representation, the redistribution method changes a vorticity distribution at time level $n$

$$\omega^n = \sum_i \Gamma_i^n \phi_\delta(\vec x -\vec x_i)$$ (4.3)

into

$$\omega^{n+1} = \sum_i \sum_j f^n_{ij} \Gamma_i^n \phi_\delta(\vec x-\vec x_j) \ .$$ (4.4)

at the next time level $n+1$. We would like this change to approximate the true diffusion over the time-step in some way. Our approach will be to demand that all finite wave numbers of the Fourier transform are correctly damped. This is similar to a weak formulation in which Fourier modes are used as weighting functions.

The Fourier transform of the new vorticity distribution is

$$\widehat\omega^{n+1}= \widehat\phi(k \delta) \sum_i \Gamma_i... ...sum_j f^n_{ij} e^{-{\rm i}\vec k\cdot (\vec x_j- \vec x_i)}\ .$$ (4.5)

This is to be compared with the Fourier transform of the exactly diffused vorticity:

$$\widehat\omega^{n+1}_e = \widehat\phi(k \delta) \sum_i \Gamma_i^n e^{-{\rm i}\vec k\cdot\vec x_i} e^{-k^2\nu\Delta t}\ .$$ (4.6)

The two Fourier transforms cannot be equal for all values of $k$ using only a finite number of vortices. However, within the neighborhood of vortex $i$, the distance $\vert\vec x_j- \vec x_i\vert$ is a small quantity of order $O(\sqrt{\Delta t})$, compare (4.1) and (4.2). This makes it possible to approximate the trailing exponentials in the two Fourier transforms by a truncated Taylor series. It does turn out to be possible to equate the Fourier transforms using these truncated Taylor series. The detailed derivation is given in Appendix A.

The resulting equations are the redistribution equations we were looking for. They involve scaled relative vortex positions defined as

$$\vec \xi_{ij} \equiv \frac{\vec x_j - \vec x_i}{h_v}\ ,$$ (4.7)

that are bounded by the neighborhood radius; $\xi_{ij}\le R$.

In terms of the scaled coordinates, the final redistribution equations are

$$O(1): \sum_j f^n_{ij} = 1 \quad ;$$ (4.8)

$$O(\Delta t)^{1/2}: \sum_j f^n_{ij} {\xi_1}_{ij} = 0 \quad ; \quad \sum_j f^n_{ij} {\xi_2}_{ij} = 0 \quad ;$$ (4.9)

$$O(\Delta t): \sum_j f^n_{ij} {\xi_1^2}_{ij} = 2 \quad ; \quad \sum_j f^n_{ij} {\xi_1}_{ij} {\xi_2}_{ij} = 0 \quad ; \quad$$

$$\sum_j f^n_{ij} {\xi_2^2}_{ij} = 2 \quad ;$$ (4.10)

$$O(\Delta t)^{3/2}: \sum_j f^n_{ij} {\xi_1^3}_{ij} = 0 \quad ; \quad \sum_j f^n_{ij} {\xi_1^2}_{ij} {\xi_2}_{ij} = 0 \quad ;$$

$$\sum_j f^n_{ij} {\xi_1}_{ij} {\xi_2^2}_{ij} = 0 \quad ; \quad \sum_j f^n_{ij} {\xi_2^3}_{ij} = 0 \quad ;$$ (4.11)

$$O(\Delta t)^{m/2}: \hbox{ Higher-order moment equations, } m = 4, \ldots, M+1 \ . \$$ (4.12)

From these equations, the redistribution fractions $f^n_{ij}$ are to be found.

Consistency requires that the numerical solution approximates the $O(\Delta t)$ diffusive changes in the exact solution: the redistribution fractions must at least satisfy (4.8) through (4.10). This results in a truncation error of order $O(h_v)$. Subsequent equations, (4.11), (4.12), can be included to achieve a higher order of accuracy $O(h_v^M)$.

Thus in principle, the accuracy can be increased arbitrarily, although for the Navier-Stokes equations the splitting error also has to be considered. The conditioning of the above system of equations also needs to be taken into account in a practical application. The equations could be recast in terms of orthogonal polynomials such as Legendre polynomials to improve the conditioning. On the other hand, the conditioning of the system may not be very important; the requirement is not to find a particular solution for the fractions $f^n_{ij}$, but to satisfy the equations accurately. In the numerical results in this paper, we simply solved (4.8) through (4.10) in the form shown.

The redistribution equations are similar to the equations obtained when a Taylor series expansion of the exact solution is substituted into a finite difference formula, or to the moment conditions in the particle methods. In fact, consistency of a finite difference scheme requires the same agreement for finite wave numbers; see Strikwerda ([219] (10.1.3)), for example. For uniform point spacing and redistribution fractions, the redistribution method is equivalent to an explicit finite difference scheme. The redistribution equations do not involve the smoothing function $\phi_\delta$ in (4.3). This allows us to choose this function after the actual computation has already been completed.

In implementing the redistribution scheme, it is important to realize that not all solutions $f^n_{ij}$ to (4.8) and following will lead to a convergent approximation. For example, a consistent but unstable explicit finite difference scheme would satisfy the equations. Some form of stability condition needs to be imposed; following Van Dommelen [235], we will demand that all fractions are positive:

$$f^n_{ij} \geq 0 \ .$$ (4.13)

This ensures that the $l_1$ norm $\Gamma^n = \sum_i \vert\Gamma_i^n\vert$ of the circulation cannot grow.

In the next two sections we will further justify the above conditions using physical and mathematical arguments. However, the truly relevant questions are clearly whether the equations are solvable, whether they can be solved using only a finite number of neighboring points within a finite scaled distance $R$, and whether the numerical solution approaches the exact solution with the expected rate of convergence. In the following chapters 5 and 6 we will prove that the answer to all these questions is affirmative for the linear Stokes equations. To verify that our method also works for the nonlinear Navier-Stokes equations, we will present example computations with nontrivial convection effects in chapters 7 and 8.