JOURNAL OF THEORETICAL

AND APPLIED MECHANICS

42, 2, pp. 223-238, Warsaw 2004

RANDOM VORTEX METHOD FOR THREE DIMENSIONAL
FLOWS. PART II: NUMERICAL IMPLEMENTATION AND

SAMPLE RESULTS 1

Andrzej Styczek

Piotr Duszyński

Marta Poćwierz

Jacek Szumbarski

Institute of Aeronautics and Applied Mechanics, Warsaw University of Technology

e-mail: jasz@meil.pw.edu.pl

Continuation of the first part of the paper focuses on implementation is-
sues, computational efficiencyandpresentationof test calculations. Sam-
ple results have been obtained for a viscous flow past a spherical body
immersed in an uniform stream.The solution algorithm for the external
Neumannboundaryproblem, the construction of the vortexparticles ne-
ar the material boundary and the evaluation of the stretching effect are
described in some details. The problem of design of efficient algorithms
for induced velocity computation is discussed briefly. The presented re-
sults include patterns of instantaneous velocity and vorticity field aswell
as selected streamlines showing the complexity of flow near the aft part
of the body and inside the aerodynamic wake.

Key words: vortex methods, vortex stretching, Fokker-Planck-
Kolmogorow equation, Itô stochastic differential equations

1. Introduction

In the first part of the paper Styczek et al. (2004), the Lagrangian vortex
method for 3D viscous flows has been formulated. Here, certain implemen-
tation issues are discussed and sample results are presented. The focus is on
the investigation of the algorithm features rather that on solving a particular

1
Part I: Mathematical background has been published in Journal of Theoretical and

AppliedMechanics, 42, 1, 2004



224 A.Styczek et al.

physical problem. Our test case is the flow past a spherical body, which is
geometrically simple but still nontrivial due to appearance of complex vortex
patterns of the wake behind the body. Simple geometry of the flow domain
facilitates an efficient numerical solution to the external Neumman boundary
problem for the Laplace equation, which is conceptually the most complex
element of the proposedmethod. An other essential element – solving the Ito
stochastic differential equations for vortex trajectories of the particles – is not
a difficult task. Since the diffusion coefficient is constant, these equations are
in fact of the Stratonowicz type, and as such can be integrated using standard
numerical procedures (see Gardiner, 1990). The crucial problem is a compu-
tational rather that theoretical one: how to calculate quickly the induced part
of the velocity field.
The sample results presented in this work include patterns of instantane-

ous velocity and vorticity fields, instantaneous streamlines and vortex spatial
distribution of the particles. The Reynolds number assumed in the numerical
simulations is the range where the wake behind the sphere becomes essen-
tially unstable, loses symmetry and the small-scale motion inside it develops
complicated, apparently chaotic behavior.
The obtained results show good andweak points of the proposedmethod.

The sources of difficulties and measures to avoid them are discussed in some
extent.

2. The numerical method for the external Neumann boundary
problem

In order to calculate auxiliary potentials ϕU, ϕO and ϕki, needed for the
velocity determination (Eqs (4.7), (6.5), (6.6) and (6.8) of the Part I), the
external Neumann boundary value problem for the Laplace equation has to
be solved.
The boundary condition for each of these problems can be written in the

following form
∂ϕ

∂n
=
∂ϕ

∂r

∣

∣

∣

r=R
= f(θ,γ) (2.1)

In the above, ϕ denotes a harmonic function in the exterior domain
|r| ≡ r > R, and the function f is given and depends on the particular
case considered. All boundary problems are well posed, meaning that

∮

f(θ,γ)cosθ dθdγ=0 (2.2)



Random vortex method for 3D flows. Part II 225

The symbols θ and γ denote the spherical coordinates on the sphere: the
longitudinal and azimuthal angles, respectively.
The solution can be expressed in the well-known form of the infinite series

ϕ=−R
∞
∑

n=1

Yn(θ,γ)

n+1

(R

r

)n+1
(2.3)

In the above, Yn(θ,γ) is a spherical function of the nth order. With the use
of boundary conditions (2.1), one gets

f(θ,γ)=
∞
∑

n=1

Yn(θ,γ)=
∞
∑

n=1

n+1
∑

m=0

Cmn P
m
n (cosθ)e

inγ (2.4)

The symbol Pmn denotes associated Legendre polynomials, i.e.

Pmn (x)= (1−x2)
m

2

dmPn(x)

dxm

where the Legendre polynomials Pn(x) are defined as

Pn(x)=
1

n!2n
dn

dxn
(x2−1)n

The expansion coefficients Cmn can be found by means of the standard
harmonic analysis on the sphere (Chorlton, 1969).

3. Computing the trajectories of the vortex particles

Asdescribed inPart I (see Section 5,Eq. (5.5)), the trajectory of the center
of a vortex particle is described by the stochastic differential equation

dri =vidt+
√
2ν dWi (3.1)

The increment of the Wiener process dWi can be defined as

dWi =
√
dt ξi

where ξi is the randomvectorwhose components are statistically independent
variables with the standardGaussian distribution N(0,1).
The velocity at the particle center vi is computed by summing up the

induction and other components accordingly to formula (4.7) from Part I.



226 A.Styczek et al.

We also already mentioned that stochastic differential equations (3.1) can be
integrated using standard numerical methods (the viscosity of a fluid is con-
stant). The trajectories of the vortex particles are traced using themulti-step
integration scheme of the 3rd order, except for newly created particles, whose
trajectories are initially calculated using a low-order single-step scheme. The
important problem is the choice of initial conditions, i.e. the initial locations
of the newly created particles. When these points are projected orthogonally
onto the sphere, a surface grid is obtained. Putting things the opposite way,
each surface grid on the sphere canbeused to generate a set of initial positions
for the vortex particles. It seems, however, that some types of the grid (like
the geographical one or the grid based of the polyhedral division) may not
be appropriate, because they would introduce particular symmetries, which
might lead to artifacts in the vorticity and velocity fields. Therefore, a grid
generation algorithm based on the maximization of the following measure

D=
∑

i,j

|rsrfi −r
srf
j | (3.2)

has been implemented. In the above, r
srf
i −r

srf
j is the vector linking ith and

jth points of the surface grid on the sphere.
The set of surface points maximizing the measure D can be determined

iteratively. In each iteration, the surface points are moved along the vectors

wi =
∑

j

ri−rj
|ri−rj|

ε

(ε denotes a small parameter introduced in order to restrict the magnitude
of the displacements of the surface points), and then projected orthogonally
back onto the sphere. Note that the vectors wi are parallel to the gradient
of the measure D, hence the described procedure leads to magnification of
this measure. Numerical computations show that the mesh generation algo-
rithm works correctly and efficiently. The iterations are terminated until the
modification of D becomes negligible. Since the obtained locations are not
changed during flow simulation, the numerical cost of the above algorithm is
not relevant.
Having the grid points on the spherical surface, one can consider the set

of lines emanating from these points and pointing outside in the normal direc-
tions. New vortex particles are initially located on these lines and the distance
from the surface is such that the tangent component of the induced velocity is
maximal. Note that the radius of the particle σ (see formula (3.8) fromPart I)
and the optimal distance are connected accordingly to the form of the core
function f(r).



Random vortex method for 3D flows. Part II 227

4. Computation of the stretching effect

Theeffect of thevortex stretching is describedbyEq. (5.7) given inStyczek
et al. (2004)

∂tω=(ω ·∇)v≡
d

dλ
v(t,r+λω)

∣

∣

∣

λ=0

Using the representation of the vorticity field for a single vortex particle
we obtain (see Section 6 of Part I)

dΩk
dt
=
3

2f(0)

d

dλ
v
(

t,r+λω(t,rk)
)

∣

∣

∣

λ=0
(4.1)

The same result can be derived the otherway. Indeed, integration of the stret-
ching equationwithinaballwith the radius r0, followedbycalculations similar
to those in Section 3 in Styczek et al. (2004), yields

∫

∂tωk(t,ξ) d3r=4πΩ
′

k(t)
[

r0
∫

0

ξ2f(ξ) dξ−
r0
∫

0

F(ξ)

ξ
dξ
]

−

−4π
3
Ω′k

[

r0
∫

0

ξ2f(ξ) dξ−3
r0
∫

0

F(ξ)

ξ
dξ
]

=

∫

S(r0)

(nω)v(t,r) dS(r0)

Using the following relation between the functions f(r) and F(r)

F(r)=

r
∫

0

ξ2f(ξ) dξ

some terms in the above expression cancel out, and finally we get

Ω
′

k(t)=

∫

S(r0)

(nω)v(t,r) dS(r0)

8π
3

r0
∫

0
ξ2f(ξ) dξ

(4.2)

Assuming r0 very small, the following approximation is admissible

Ω
′

k(t)≈
(ω∇)v

2
r3
0

r0
∫

0
ξ2f(ξ) dξ

=
(ω∇)v
2
3
f(0)

The velocity vector v contains the contribution induced by the kth particle.
We will derive closed formulae describing the part of the vortex stretching



228 A.Styczek et al.

due to the interaction between the jth and kth particles. The sum of these
contributions over the indices i=1,2,3, . . . defines the total stretching effect
of the vorticity-induced component of the velocity field on the kth particle.
The remaining part of the stretching is due to the potential component and
can be computed by numerical differentiation in the direction of the vorticity
vector in the particle centre. In order to derive necessary formulas, we use a
local coordinate systemwith the origin at the particle centre. Next, we have

s(r,j,k) =
d

dλ

Ωj × (r+λΩk)
|r+λΩk|3

F(|r+λΩk|)
∣

∣

∣

λ=0
=

(4.3)

=
[Ωj ×Ωk
r3

−
3(Ωkr)(Ωj ×r)

r5

]

F(r)+
(Ωkr)(Ωj ×r)

r2
f(r)

The problem arises when the interacting particles have their centers in the
same point. If r → 0, the second and the third terms in (4.3) cancel each
other since

lim
r→0

F(r)

r3
= lim
r→0

1

r3

r
∫

0

ξ2f(ξ) dξ=
1

3
f(0)

The first term in (4.3) will vanish if only f(0) = 0. The reason for such an
assumption is the avoidance of the pathological self-stretching effect.
Summing up all contributions defined by (4.3), we get

S(k)=
∑

j

s(rk−rj,j,k) (4.4)

whichdetermines thevorticity-inducedpart of thederivative Ω′k(t).The stret-
ching effect due to thepotential part is calculatednumericallyusing the central
finite difference formula of the secondorder in thedirection of theunaryvector
Ωk/Ωk.

5. Computation of the vortex interaction (induced velocity)

Usability of the method of vortex particles in 3D is mostly determined by
its computational efficiency. Like in any particle method, the numerical cost
of induced velocity evaluation is proportional to the squared number of the
vortex particles. The stretching determination is also expensive since it requ-
ires additional calculation of the velocity field at two points near each particle
center. Another time-consuming part of the algorithm is the determination



Random vortex method for 3D flows. Part II 229

of the right-hand side vector of the linear system for the charges of the new
particles. Although its cost increases only linearly with the number of partic-
les, there is a large computational ”overhead” due to complicated expressions
involved. Other elements of the algorithm give rather negligible contributions
to the overall CPU effort.

The most crucial for computational effectiveness is the reduction of time
for the induced velocity evaluation. It seems that the simplest idea is to intro-
duce an interpolation grid.Themutual inductionof close particles is computed
from exact formulas on each-to-each basis. The induction of distant particles
is calculated indirectly: first the velocity at the grid vertices is determined and
next it is interpolated from the vertices to the particle centers. The major
drawback of such approach is that simple interpolation rules cannot be di-
vergent free and vice versa. Besides, this method works efficiently only if the
spatial distribution of the particles is sufficiently uniform, which is rarely a
case, especially for wake flows behind immersed bodies.

Other, much more sophisticated approaches to interacting particles have
been proposed since late ’80s. Themost popular ones are:

• methods based on the usage of multipole expansions, the most eminent
example being the famous Greengard-Rokhlin’s Fast Multipole Method
(FMM) (see Pringle (1994) for comprehensive survey of this approach),

• methods using tree structures, for example the Bernes-Hut algorithm
(see Blelloch and Narlikar, 1997).

Both algorithms have their pros and cons. The BH algorithm is easier to im-
plement but it offers worse performance, especially when higher accuracy is
demanded Blelloch and Narlikar, 1997. With a moderate number of particles
(say, 50 thousands), the numerical cost related to error reduction by one order
of magnitude in the BHmethod is nearly the same as the cost of error reduc-
tion by three orders in FMM (Gibbon, 2002). The main disadvantage of the
latter is, in turn, high complexity of (recursive) implementation, which results
in a large computational ”overhead”. In effect, the FMM can significantly ac-
celerate computations (say of one order ofmagnitude), but only if the number
of the vortex particles increases well over fifty thousands (Gibbon, 2002).

In the currentwork, thenumerical simulations havebeen carriedoutwitha
fewhundreds of the vortex particles injected into the flowdomain in each time
step.Thenumberof theparticles increases almost linearly andachieves several
tens of thousands at the end of the simulation. With the currently available
hardware, we have found it impractical to carry on with a larger number of
particles. In viewof the efficiency limitations of the fast summation algorithms



230 A.Styczek et al.

mentioned above, it seemed reasonable to use a directmethod for the induced
velocity. In order to reduce computational times, partial parallelization of the
codes using MPI libraries have been attempted. The parallelization has been
applied to the most time consuming part of the computations, namely:

• evaluation of the mutual induction of the vortex particles,
• calculations of the stretching effect due to the induced velocity (the
part of the stretching due to potential components has been calculated
sequentially),

• calculations of the right-hand sides of the linear system for newly created
particles.

The parallel computations have included three stages. In the first stage, all
necessary data concerning the particles have been sent over to all processors.
Next, each processor performed calculation for its individual portion of the
particles. Finally, all results have been sent back to the main processor.

6. Numerical test case: flow past a spherical body

The vortex method described above has been tested extensively on the
example of a viscous flowpast a spherical body immersed in a uniform stream
with u∞ = 1. The obtained results include trajectories of the particles and
instantaneous patterns of the velocity and vorticity. In order to visualize the
flow complexity, selected streamlines of an instantaneous velocity have also
been calculated.

In all presented cases, the flow has been accelerated from rest. During the
initial short time interval, the free streamvelocity u∞ was increased fromzero
to unity and then kept fixed. The time step ∆t=0.01was found a reasonable
compromise between the accuracy and computational efficiency.With a smal-
ler time step, the number of the vortex particles grew much faster, making
long lasting simulation not feasible because of the hardware limitations.

As it has been already described, the vortex particles are shed from the
region near the body surface and then moved downstream forming the wake
behind the body. A large portion of the particles drops into the area right
behind the body and remains there for a rather long time. Such behavior
is natural from a dynamical point of view – the area of the near wake is
occupied by a large, nearly toroidal, vortex structures which take a mass of
the flow from the separated boundary layer and push it towards into thewake



Random vortex method for 3D flows. Part II 231

Fig. 1.

interior and into the aft part of the body. From the computational point of
view, this phenomenon is an example of the characteristic feature of all vortex
methods: they have natural ability to self-adapt the spatial resolution in the
dynamically active regions of the flow field. Figure 1 shows instantaneous
locations of the vortex particles obtained for the Reynolds number Re= 200
and the simulation time t=16.Thenumber of the particles generated at each
time step in the vicinity of the body surface was set to 400. The number of
the particles in the flowfieldwas about 40 thousands.Darkly shaded area just
behind the body shows how large portion of the particles is ”trapped” in the
near-wake structures.Another important feature is the asymmetric formof the
wake – its wavy shape, observed also in other simulations and in experiments,
is apparent.

The instantaneous patterns of the velocity field are presented for three dif-
ferent values of the Reynolds number. In Fig.2, the projection of the velocity
on the XZ plane is presented. The Reynolds number is equal: 100 (Fig.2a),
200 (Fig.2b) and400 (Fig. 2c). In all cases, the number of theparticles genera-
ted in each time step was equal to 400, and all velocity patterns corresponded
to the same time instant t=9.2.

The flow field obtained for Re=100 is remarkably symmetric. The vortex
structure in the form of a flattened toroid, extending along the wake at a
distance of about one diameter of the sphere is visible. The obtained wake
symmetry is in good agreement with the experimental observation for the
same Reynolds number.

The comparisonwith the velocity fields obtained for larger Reynolds num-
bers reveals gradualdevelopment of theasymmetry in thewakeflow.Moreover,
the velocity irregular profiles of the reversed flow in the central part of thewa-
ke observed for Re = 400 appear due to a more complicated and containing
more that one large vortex structure flow pattern. The development of the
wavy shape of the wake in the streamwise direction can be inferred from the
most right velocity vectors in Fig.2c (for X coordinate between 6 and 7). The



232 A.Styczek et al.

Fig. 2.



Random vortex method for 3D flows. Part II 233

shape of the flow at Re = 200 exhibits an intermediate character: it is still
very regular but the first symptoms of the symmetry breaking are already vi-
sible. The overall picture of thewake flowdynamics in the considered range of
the Reynolds number seems to be qualitatively similar to the results of recent
vortex simulation by Ploumhans et al. (2002).

Fig. 3.

For better visualization of the flow structures, the contour maps of the
magnitude of the XZ-projection of the velocity and Y -component of the
instantaneous vorticity field are presented in Fig.3 and Fig.4, respectively.
The corresponding flowhas been calculated for Re=100 and the time instant
t = 17. The large-scale symmetry seen before on the vector maps, especially
in the near wake, is now less apparent, whichmay be partly explained by the
choice of a much later time instant. The vorticity map presents some small-
scale details of theflow. It turnsout that regions of largest vorticity are located
near the body surface (as expected) andon the ”edges” of thewake.The latter
observation is rather not trivial as onemight expect that the area of the large
vorticity would be located closer to the center line of the wake, namely in the
vicinity of points where the direction of the wake flow is altered.

In spite of the apparent global symmetry, the small-scale motion obtained
in the vortex simulation can be remarkably complicated. InFig.5we show two



234 A.Styczek et al.

Fig. 4.

Fig. 5.



Random vortex method for 3D flows. Part II 235

streamlines of the instantaneous velocity field calculated for Re= 200 at the
time instant t=17. The streamlines originate from two points located in the
upstreamdirection and close to each other. As long as the streamlines go near
the body surface they remain close and have a similar shape. The situation
changes drastically when they enter the wake region. The dynamics in the
”frozen” velocity fieldbehind the body (here the streamlines are interpreted as
trajectories of themarker particlesmoving in pseudo-time in the fixed velocity
field) is highly sensitive to initial conditions meaning that two, initially close,
trajectories diverge rapidly wiggling around the whole area in a complicated
manner.

While the presented results show the capability of the proposed vortex
method tomodel awake flowbehinda bluffbody realistically, they also reveal
itsmajorweakness: poor quality of representation of thin shear layers near the
material surface. It is evident from the vector maps in Fig.2 than the velocity
gradients at the surface are much smaller that in a real flow (also the general
tendency, i.e. ”the larger Re the larger gradients” isweaklymarked). In effect,
the friction drag is vastly underestimated. Also, it is difficult to say to what
extent the shift of the large vorticity off the wake center is merely the artifact
of the thick boundary layer and its premature separation.

The source of the above difficulty is in the particle’s geometry and, partly,
in hardware limitations. In a sense, each vortex particle occupies the whole
space (see Part I) but it is mostly concentrated in a rather small spherically
symmetric vicinity of its central point characterized by the magnitude of the
”core” radius σ. Such vorticity carriers are very convenient to use because
of simple analytical induction formulae. Besides, in contrast to the ”vortex
particles” used by some researchers, the vorticity of the particles constructed
in this work is a well-defined divergent-free vector field. On the other hand,
the spherical symmetry of these objects does not make them well suited for
approximations of a flow field with highly anisotropic geometry in a kind of
a boundary layer, unless their core radii are very small. The latter, however,
brings huge computational efficiency problems. If the particles are small then
the number of the vortex particles necessary to enforce the no-slip condition
with a reasonable accuracy becomes very large. The reason for this is rather
obvious. If the body surface was loosely covered with a small number of tiny
particles then theaccuracyof cancellation of the tangential velocity component
(set to zero only in a small number of collocation points) would be very poor.

A possible way of overcoming the efficiency limitations and achieving a
satisfactory resolution in the transversal direction of the boundary layer is to
introduce some ”flattened” vortex structures. The crucial question is how to



236 A.Styczek et al.

construct such particles, so that the corresponding induction formulae (pro-
bably not explicit anymore) can be computed quickly. Additionally, such an
idea is a source of further arbitrariness in the method, namely, the problem
ariseswhenandhowa ”flat” vortex particle is to be replaced by an ”ordinary”
(spherical) one.

A radical approach is to resign from the particle approximation of the flow
in the closest vicinity of the body and use the concept of a singular vortici-
ty layer. The no-slip boundary condition is realized by a surface distribution
of discontinuity of the tangent velocity component determined as a solution
to the appropriate boundary integral equation. The magnitude of a local di-
scontinuity determines the equivalent local flux of the vorticity through the
body surface. The vorticity production (or destruction) is realized by injection
(or removal) of the adjacent vortex particles and/or the modification of the
”vorticity charges” of the particles from the neighborhood. There is also a
possibility to generate singular structures (in the form of small planar vortex
sheets), which are replaced by ”ordinary” particles if they move farther away
from the boundary. Some of these ideas have been recently implemented and
tested by Ploumhans et al. (2002).

7. Summary and final remarks

Themethod of vortex particles proposed in this study has been tested on
the standard example of a flowaround a spherical body atmoderateReynolds
numbers. In this part of the paper, we have focused on some technical and
implementation details. The presented results are of preliminary character.
Due to computing hardware limitations we were able to carry out calcula-
tions with a rather small number of the vortex particles. These circumstances
affected the spatial resolution, but to a different extent in different flow re-
gions. The dynamics of structures of the velocity in the wake area seems to
bewell reproduced.The large-scale motion remains in general agreement with
the experimentally observed patterns of the flow. The major problem is the
insufficient resolution of thin shear layers and the boundary layer at the body
surface. It results in a serious underestimation of the surface-transversal velo-
city gradient and, consequently, friction forces. This problem appears because
the vortex particles used in the numerical simulations were too large (in the
sense of the spatial concentrations of their cores). Usingmuch smaller particles
would improve the spatial accuracy, but at the expense of currently intolerable
increase of computational costs. Othermethods of the enforcement of the no-



Random vortex method for 3D flows. Part II 237

slip condition, based on some surface-fitted or even singular vortex element,
can be incorporated in the future development of the proposedmethod.

Acknowledgements

This work has been supported by the State Committee for Scientific Research

(KBN), grant No. 8TO7A02520.

References

1. Blelloch G., Narlikar G., 1997, A practical comparison of N-body algo-
rithms – parallel algorithms, Series in Discrete Mathematics and Theoretical
Computer Science, 30

2. Chorlton F., 1969,Boundary Value Problems in Physics and Engineering,
Van Nostrand Reinhold Co., London

3. Gardiner C.W., 1990, Handbook of Stochastic Methods, 3rd Ed. Springer
Verlag

4. Gibbon P., Sutmann G., 2002, Long-range interactions in many-particle si-
mulations,QuantumSimulations of Complex Many-Body Systems, 10, 467-506

5. Ploumhans P., Winckelmans G.S., Salmon J.K., Leonard A.,
WarrenM.S., 2002,Vortexmethods for direct numerical simulations of three-
dimensional bluff-body flows: application to the sphere at Re=300, 500, and
1000, Journal of Computational Physics, 178, 427-463

6. Pringle G.J., 1994, Numerical study of three-dimensional flow using fast pa-
rallel particle algorithms, Ph.D. Dissertation, Departament of Mathematics,
Napier University, Edinburgh

7. Styczek A., Duszyński P., Poćwierz M., Szumbarski J., 2004, Random
vortex method for three dimensional flows. Part I: Mathematical background,
Journal of Theoretical and Applied Mechanics, 42, 1, 3-20

Stochastyczna metoda wirowa dla przepływów trójwymiarowych.
Część II: Numeryczna implementacja i przykładowe wyniki

Streszczenie

W drugiej częśći pracy omawiane są zagadnienia związane z implementacją me-
tody wirowej w trzech wymiarach oraz prezentowane są wybrane wyniki obliczeń.



238 A.Styczek et al.

Testowymprzykładem jest opływ kuli strumieniem jednorodnym.Opisano konstruk-
cję solwera zewnętrznego zagadnienia brzegowego typu Neumanna, proces generacji
nowych cząstek wirowych w pobliżu powierzchni ciała oraz sposób obliczania ”stret-
chingu”. Omówiono problem efektywnego obliczania prędkości indukowanej. Przed-
stawione wyniki, obejmujące chwilowe pola prędkości i wirowości oraz wybrane linie
prądu, demonstrują złożoność pola przepływu w śladzie za opływanym ciałem.

Manuscript received January 23, 2004; accepted for print February 27, 2004