6.2.6 Evaluation of the vorticity

As was discussed in section 4.1, in the redistribution method the diffusion of vorticity is achieved by changing the strengths $\Gamma ^n_i$ of the vortices. Unlike some other methods, the continuous vorticity field

$$\omega(\vec{x},t)\;=\;\sum_{i}\,\Gamma_i\, \phi_\delta(\vec{x}-\vec{x}_i(t)) \quad \ ,$$ (6.10)

is not differentiated or even evaluated during diffusion. The equations 4.8 governing diffusion are completely independent of the ``vortex core shape" $\phi_\delta$ in the representation of the vorticity field above. Thus the vortex core $\phi_\delta$ does not appear in the actual computation. To be precise, while the diffusion process is independent of the vortex core, the implementation of convection of section 6.1 still requires one. Our choice (6.2), (6.3) for the core shape $\phi_\delta$ has a relatively low second order accuracy [106], but it is everywhere positive. The reason we chose this function is that unsteady separating flows at large Reynolds numbers involve short scale vorticity features that can be lost by large core sizes. This makes a small core size desirable regardless of the order of accuracy of the core. For a small core size, a relatively low order of accuracy can be sufficient; and other considerations may be more important. In particular, positive second order cores such as the one we used have the advantage that they cannot introduce false vorticity of opposite sign.

Another place we use the vortex core is when we need to evaluate pointwise vorticity values for output purposes. In this case the considerations for the choice of the most desirable core are somewhat different. For maximum visual smoothness, a large core is desirable, since a large core gives the greatest reduction in short wave errors. These short wave errors are further also much more pronounced in the vorticity field than in the velocity evaluation during convection. On the other hand, there is much less risk that a larger core would smooth small features, since the actual solution is now known, and the effect of core size can be determined experimentally without repeating any of the computation. Further, even if there would be some loss of information about the shorter wave lengths in the output for the vorticity, this loss does not affect the further computation. Such considerations suggest the use of a second core different from the one used for convection to do the output. In fact, there is no good reason why the two cores would need to be the same. So, for evaluation of the vorticity, we choose a second core with a high order of accuracy, since these can be larger for a given accuracy. Thus, for the computations of flows in free space in sections 7.1 and 7.2, we choose our second smoothing function to be a relatively large but infinite-order core

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

in which

$$\phi(\vec{x}) = \frac{1}{2\,\pi} \; \frac{J_1(\vert\vec{x}\vert)}{\vert\vec{x}\vert} \quad \ ,$$ (6.12)

where $J_1$ is the Bessel function of the first kind of order one. This smoothing function was first proposed by Leonard [126].

As shown in section 7.2, the results are not very sensitive to the precise choice of either of our core shapes. In most computations in this thesis, as in [201], a core size $\delta$ of $2.4\sqrt{\nu\,\Delta t}$ to $3.5\sqrt{\nu\,\Delta t}$ was used, where $\nu$ is the kinematic viscosity and $\Delta t$ is the time step.

Near boundaries, the determination of the pointwise vorticity field runs into difficulties, since the vorticity field beyond the boundary is unknown. In order to compute the correct vorticity near a boundary, the vorticity field must be extrapolated into the boundary in some way.

For our computations of flows around circular cylinders of unit radius in chapter 8, we extrapolated the vorticity into the cylinder by mirroring every vortex into the cylinder to a radial position that is the inverse of the radial position of the original vortex. Further, we changed the strengths of the thus mirrored vortices so that a constant vorticity level outside the cylinder would be extrapolated into a constant vorticity level within the cylinder. It is easily seen that this requires that the strength of each mirrored vortex is reduced by a factor equal to the radial position of the vortex to the power four.

Using this procedure, vorticity fields that are about constant near the wall can be evaluated without difficulty. However, if there are appreciable gradients, they will affect the evaluation of the vorticity. This is evident in our computations of the flow about a cylinder in rotational oscillations in figure 8.2. As expected, our procedure produces a local average of the vorticity near the wall, rather than a pointwise vorticity value at the wall.

To remove such errors would require that the vorticity field is linearly extrapolated into the cylinder, rather than as a constant. This would require some additional coding effort, but does not seem particularly difficult from a fundamental point of view. For a cylinder in steady rotation, figure 8.1 the vorticity flux vanishes at the wall and the error in our simple procedure is much smaller.

One additional modification needed near solid walls concerns the smoothing function. The smoothing function (6.11) decays too slowly at large distances to be useful for flows with solid boundaries; for these flows the vorticity can only be extrapolated a small distance into the boundary. Instead, for the flows over circular cylinders in chapter 8 we chose $\phi$ to be a fourth-order Gaussian of form [17],

$$\phi(\vec{x}) \; = \; \frac{1}{\pi}\;\left( 2\, e^{-\vert\ve... ...\;-\; \frac{1}{2}\,e^{-2\vert\vec{x}\vert^2} \right)\quad \ .$$ (6.13)

In these computations a core size $\delta=3.5\sqrt{\nu\,\Delta t}$ was used.

Figure 6.4: Smoothing functions in (a) Fourier space and (b) Physical space. Broken lines are nonconvergent smoothing. Solid lines are modified smoothing.

A final smoothing function was used to evaluate the vorticity in three-dimensional flows in free space in section 7.3. For two-dimensional flows in free space in sections 7.1 and 7.2, Leonard's infinite order smoothing function (6.12) above gives excellent results. In three dimensions, the equivalent infinite order smoothing function is

$$\phi(\vec{x}) = \frac{1}{(2\,\pi)^{3/2}} \; \frac{J_{\frac{3}{2}}(\vert\vec{x}\vert)}{\vert\vec{x}\vert^{3/2}} \quad \ ,$$ (6.14)

where $J_{\frac{3}{2}}$ is the fractional order Bessel function of the first kind of order $\frac{3}{2}$. However, this function decays too slowly at large distances to be useful: even constant vorticity fields cannot be represented by it. The reason for the slow decay of the smoothing function (6.14) is its Fourier transform, which is shown as a dotted line in figure 6.4(a). It is unity for wave number $k$ less than one and vanishes for larger $k$. The resulting discontinuity causes the slow decay of the smoothing function itself. To obtain faster decay, we must smooth the discontinuity. We choose:

$$\widehat{\phi}(k)= \frac{1}{e^{2(k^2-k^{-2})} + 1} \ ,$$ (6.15)

where $k$ is the wave number in the radial direction. This produces an infinite order smoothing function that also decays exponentially at large distance. Since this core cannot readily be evaluated in physical space, we approximated it by a spline interpolant (figure 6.4b).

To apply the numerical implementation discussed so far to flows over solid walls we need to address the numerical handling of the boundary condition on solid walls; we will describe that in the next section.