8.2.9 Performance compared with previous methods

In the previous subsections 8.2.2 through 8.2.6 we established that our results are invariably comparable to the most ones accurate available. In this subsection we will now try to argue the actual superiority of our method.

As was discussed in chapter 1, our method was designed to be a replacement of the random walk method. Like that method it has the capability to use independent points, a capability that none of the other methods have. As we mentioned there, this capability is crucial to deal effectively with the strong distortions of computational points caused by the strong convection during unsteady separation process [241]. Our primary objective was to maintain these advantages of the random walk method, but to do something about its large and random numerical errors.

The key question is therefore whether we succeeded in doing so, in particular whether the small amount of additional work we introduce is worth it. Based on the results in this section, we can answer with a resounding yes. We need merely point to our own results for $\Delta t=0.04$ and $\epsilon _\Gamma =10^{-5}$ in figures 8.33 and 8.34 as compared to the random walk results in figure 8.23. These two computations use the same number of vortices, yet our new results are dramatically better. And so are other measures, such as the drag (see figures 8.28 and 8.29).

Note that this comparison is fair: we used the same fast velocity summation scheme as the random walk computations, we used the same handling of the wall boundary condition, the same evaluation of the drag, and so on. The only essential difference is that we replaced the random walk with the equally powerful redistribution method. (The random walk results of Cheer [42] use only 900 vortices and we would not expect more than qualitative features for this little resolution).

Wang [245] used the random walk method to study the control of flow separation from aircraft wings. He complained that the large and random errors in the random walk method make it difficult to study the effect of the various parameters. Indeed, it should be clear from mere comparison of the two random walk computations (a) and (b) in figure 8.23 that the random errors in this method must be large enough to effectively mask the effect of changing the parameters. Yet this is the exact same computation; only the random numbers changed! Our new method makes it possible for the first time to use vortex methods with independent numerical points for real applications.

In addition to this, the results of subsections 8.2.2 through 8.2.6 show an impressive performance of our method compared not just to the random walk method, but also to other much more limited numerical methods. This is truly astonishing, since the geometry of a circular cylinder, at a Reynolds number as low as $Re=9,500$, is obviously much too simple to do justice to the capabilities of our mesh-free method.

Take the Particle Strength Exchange vortex methods. The computations based on those methods shown in this chapter initially locate the vortices on a very smooth, regular mesh. Next, they take only a few computational steps before they create a new mesh of vortices. It may be expected that those computations could gain tremendous accuracy benefits from having at all times a very smooth, ordered vortex structure. After all, they give up a significant amount of flexibility for it.

In fact, the results in subsections 8.2.2 through 8.2.6 show that their results, in particular those reported by Koumoutsakos & Leonard [117] do rank among the very best. Yet they use much more vortices than we to achieve this. Figures 8.33, 8.34, and 8.36 show that our method produces converged results at $\Delta t=0.02$ and $\epsilon _\Gamma =10^{-5}$. This computation requires as little as 60,000 vortices at time $t=3.0$. In contrast, Koumoutsakos & Leonard [117] used 350,000 vortices at the same time. The fact that the particle scheme uses little computational time for diffusion above the time for convection is irrelevant for such disparities in the number of vortices.

Even with the large number of vortices used by Koumoutsakos & Leonard by [117], it appears to us that their resolution could be improved. For example, compare their vorticity field at time $t=3.0$ in figure 8.22 with our converged vorticity field. Although we attempted to get the colors identical, it still appears that some of the finest flow features are simply thicker in the particle results. This would suggest that there is still some computational smoothing of the smallest scales.

We also note that the drag coefficient predicted by Koumoutsakos & Leonard [117] in figure 8.28 has a slight ``dip" near time $t=2.0$, similar to, but not the same as, the results of Anderson & Reider [3]. Such a dip does not seem consistent with our convergence studies, such as in figure 8.29. It also does not agree with the spectral element computation of Kruse & Fischer [120], figure 8.28, using over a million computational points.

One possible theoretical reason why the particle methods, also including Fishelov's method [78], may experience more difficulties with resolving the shortest scales will be given in section 9.1. Yet, the fact remains that the particle methods did give excellent results for this flow. There is another, possibly more important difficulty faced by those methods. It has to do with reliability.

The work of Koumoutsakos & Leonard [117] was the culmination of a fairly long research effort at Caltech into particle methods starting with the work of Pépin [170]. There has been considerable variation in the predictions of the particle method for the flow over the circular cylinder during that time. For example, we mention the drag coefficients in figure 8.28. Even the results found in Koumoutsakos's thesis [119] show significantly different vorticity fields and drag coefficients from those presented in Koumoutsakos & Leonard [117]. The difficulty is probably that the particle method has a considerable number of uncertainties to deal with. For example, we mentioned above and in subsection 8.2.1 that the particle computations very frequently create fresh vortex distributions. This raises difficult questions such as how to distribute the new set of vortices, how to interpolate the strengths of the old set of vortices onto the new set, and what is the effect of the errors involved in each of these steps.

Further, the particle discretization of diffusion has a longer numerical region of influence than our redistribution method. Such regions of influence make it much harder to handle solid wall boundary conditions. In fact, the handling of boundary conditions for the Particle Strength Exchange scheme has received considerable attention [119,148,170] and the final procedure developed by Koumoutsakos, Leonard & Pépin [118] is considerably more complicated than our simple approach of section 6.3. Such more complicated procedures will pose additional difficulties for more complicated bodies than circular cylinders.

There are other uncertainties in the particle methods. For example, these methods discretize diffusion using some form of approximate integral operator. This operator requires a chosen discretization function, called the kernel [72]; it is not obvious what shape of this kernel to choose, or its general size, yet such decisions may have a major influence on the final results. Shiels [208] at Caltech is currently trying to clarify some of those issues.

Our method does not face such difficulties. We have never experienced the need to improve our results by attempting to adjust computational parameters; our results turned out correct the first time. In fact, we took most of the few computational decisions that our method requires directly from the random walk scheme. This robustness of our scheme is a matter of reliability: not for all flows is there as much reference data as for the circular cylinder.

The impressive performance of our method compared to other vortex methods raises the question how well it stands up to hybrid and finite difference methods. Being mesh-based, these methods are much more restrictive than ours, but this should not be a great concern for the circular cylinder. In addition, these methods have seen many decades of development compared to our new procedure, and can therefore make use of such advanced concepts as nonlinear viscosity, higher order schemes, fast Fourier transforms in the tangential direction, and so on. It is therefore truly remarkable to see that our method outperforms these established procedures.

One obvious example is the hybrid vortex method results of Chang & Chern [40]. Their computed drag coefficient, figure 8.28, can simply not be supported by any of the available data. The reason is easily guessed: their computation uses as little as $256\times200$ computational points. These authors do use radial mesh stretching, which produces seemingly adequate radial resolution right at the cylinder. However, a simple check shows that using their radial stretching, most of the radial resolution has already been lost in the middle of the boundary layer. Chang & Chern [40] report on checking the radial resolution by changing the resolution at the wall by ``a factor close to one"; they did not find a ``substantial difference" in doing so. Since multiplying the radial resolution by a factor equal to one will not show any difference at all, this test is not adequate. Chang & Chern [40] wrote that they have extended their work to Reynolds numbers up to $10^6$, and still obtained very stable numerical results. It is our opinion that this highly unlikely stability is caused by excessive artificial diffusion due to inadequate numerical resolution.

The finite difference computations of Wu, Wu, Ma & Wu [249] are in excellent agreement with our own, as shown in figures 8.15 and 8.28. However, a large number of mesh points was needed to achieve this. For most of their results, they used a $512\times 301$ mesh that extends outward towards 5 times the cylinder radius. Yet, their own results show that their values for the drag and velocity profiles using this mesh are not converged. In fact, they would show significant deviations from the correct results in figures 8.15 and 8.28. Hence instead, in the figures we used their results from a mesh that according to the listed mesh spacing has roughly $512\times900$ mesh points, with an exponential stretching in the radial direction. With this mesh, their results for the relatively simple cylinder flow do agree very well, with only some minor variations in the drag coefficient near the plateau.

Note however, that Wu, Wu, Ma, & Wu reach the conclusion that the numerical results will not converge as the grid refines, limiting the application of their method for very high resolution computations. Clearly, even despite the simple configuration and the long development time of finite difference procedures, there remains much uncertainty.

The computation of Anderson & Reider [3] also produced very good results (figures 8.15 and 8.28), save for a minor ``dip" in the drag coefficient that is not supported in this form by other computations. Yet, their numerical procedures are somewhat limited even for a finite difference procedure. In particular, they only actually compute a strip of half a radius thickness around the cylinder. Outside, they use the analytical potential flow Fourier series solution in polar coordinates. They refer to this as a domain decomposition procedure. Yet, in this form it does seem quite limited in applicability. Further constraints arise from their higher order procedures. They used as many as $2048\times256$ mesh points in their thin strip. (We note that Anderson & Reider lost their actual data. The data presented in figures 8.15 and 8.28 were measured from Wu, Wu, Ma, & Wu's paper [249]).

An interesting earlier procedure by Anderson & Reider [4] used the boundary layer approximation. This could be useful to reduce computational effort or to simplify boundary conditions for certain flows. However, some caution is needed here, since it has been well established that the unsteady boundary layer equations develop ``Van Dommelen & Shen" singularities. Current evidence indicates that viscous-inviscid interaction will not eliminate that singularity, in particular, see Peridier, Walker & Smith [171], Cowley, Van Dommelen & Lam [68], and Cassel, Smith & Walker [37].

The finite difference computation of Hakizumwami [105] works quite well at Reynolds number $Re=3,000$, figure 8.12, but at $Re=9,500$ the need for additional resolution starts to show up. His grid was only $300\times128$. He takes advantage of the simple geometry of the problem by using a Fourier-Galerkin method in the tangential direction. However, in this problem the smallest scales are clearly smaller than he can resolve.

The computations of Loc & Bouard [134] at $101\times301$ points and Loc [133] at $61\times61$ mesh points would need more resolution, as well as true radial mesh stretching. In addition, we do not quite understand how the initial data are generated by Loc & Bouard [134]; the procedure described in that reference does not seem to make sense. In any case, the evolution of the drag at $Re=550$ in figure 8.25 is simply unreconcilable with the analytical solutions of figure 8.24.

The conclusion we reach from all this is that the accuracy of our method is outstanding compared to the much more developed and much more limited finite difference methods. In fact, the apparent convergence of our data leads us to believe that our results may actually be better than all the finite difference computations. For example, using half a million mesh points in a small strip, Anderson & Reider [3] still experience a dip in the drag that is not supported by the other computations, figure 8.28. Koumoutsakos & Leonard have a dip in their drag, but it is a different dip. Yet, our drag is converged using as little as 60,000 vortices, figure 8.36(b).

Also using about half a million mesh points, the results of Wu, Wu, Ma & Wu [249] do not show the dip in the drag curve, but their results for the drag plateau in figure 8.28 seem somewhat too high compared to ours. We tend to believe our own results based on our convergence studies, figure 8.29, as well as based on the preliminary results of Kruse & Fischer [120] in figure 8.28. The latter computation is a spectral element scheme with over a million mesh points (ignoring symmetry). In addition, we see at least one possible explanation why the finite difference scheme may compute a plateau that is somewhat too high. Figure 8.36(a) shows that ignoring the exponentially small velocities above the boundary layer can result in an increase in the drag plateau. Now the typical length scale of those velocities decreases going radially outward, while the mesh of Wu, Wu, Ma & Wu [249], as of most other finite difference schemes, expands going outward. Hence resolution difficulties are conceivable.

This then brings us to the final question how well our results compare to the spectral element computations. A spectral scheme can be considered the very opposite of our method: while in our scheme the computational points (vortices) are completely independent points, simply put a spectral method is based on making rigid associations between each point and many of its neighbors. These associations allow spectral accuracy (or at least a very high order of accuracy) to be achieved, but at the price of a loss in flexibility. The spectral element method tries to reduce this loss in flexibility, greatly increasing the range of practical applications, but clearly it can never achieve the flexibility of our totally independent points. The methods are in this sense indeed opposite.

Yet, although our method and the spectral element method are not really competitors, we are gratified to see how good the agreement is between the two. Only a very slight difference in the level of the drag plateau in figure 8.28 is observed. Even this difference might still be due to the spectral method; the spectral method could conceivably still have a slight difficulty with the exponentially small velocities above the boundary layer. The spectral element sizes increase away from the wall, while, as mentioned, the length scales of the exponentially small velocities decrease. However, note that the spectral computations are only preliminary results kindly made available to us by the authors. We will need to know their final results before drawing any final conclusion.