Thesis

"Vortex particle-mesh with combined immersed
boundary and mesh refinement techniques :

application to bluff-body and wake-vortex flows"
Lonfils, Timothée
Abstract
The present work joins a global effort in the development of efficient and accurate
numerical tools for the simulation of complex fluid flows. More specifically, we
focus on unsteady external flows, a fluid dynamics discipline which actually
pervades applied sciences and engineering. These are the flow encountered
in aircraft or car aerodynamics, wind energy, biological locomotion, etc. In the
problems addressed here, we assume that the flow is incompressible and that the
inertial forces dominate the viscous stresses: we are dealing with both moderate
and high Reynolds number flows. Furthermore, the present work treats the flow
equations in a very distinct approach. We use a vortex particle-mesh (VPM)
method, which belongs to the broader class of “vortex methods”. Such methods
use the vorticity-velocity formulation of the Navier-Stokes equations, rather than
the velocity-pressure formulation: the vorticity field is thus the primary variable.
The major advantage comes from the com...
Document type : Thèse (Dissertation)
Référence bibliographique
Lonfils, Timothée. Vortex particle-mesh with combined immersed boundary and mesh refinement
techniques : application to bluff-body and wake-vortex flows. Prom. : Winckelmans, Grégoire
Available at:
http://hdl.handle.net/2078.1/107701
[Downloaded 2017/08/31 at 21:05:00 ]
Université catholique de Louvain
Ecole polytechnique de Louvain
Institut de Mécanique, Matériaux et Génie Civil
Vortex particle-mesh method with combined immersed

boundary and mesh refinement techniques. Application
to bluff-body and wake-vortex flows
Timothée Lonfils
Mai 2011
Composition du jury: Thèse présentée en vue de

Pr. G. Winckelmans (UCL, promoteur) l’obtention du grade de docteur
Pr. P. Chatelain (UCL) en sciences de l’ingénieur
Dr. P. Geuzaine (Cenaero)
Pr. P. Koumoutsakos (ETH)
Pr. J.-F. Remacle (UCL)
Pr. T. Pardoen (UCL, président)
©2011
Timothée Lonfils
All Rights Reserved
Acknowledgments
Je voudrais avant tout remercier mon promoteur, Grégoire Winckelmans, qui

m’a fait confiance et m’a permis de réaliser cette expérience inoubliable. Sa
rigueur scientifique et son regard critique m’ont poussé jour après jour à me
dépasser. Je remercie également chacun des membres de mon jury, Philippe
Chatelain, Petros Koumoutsakos, Philippe Geuzaine, Jean-François Remacle,
Thomas Pardoen et Grégoire Winckelmans pour les remarques avisées lors de
ma défense privée.
Je remercie le F.R.S.-FNRS (Fonds de la recherche scientifique) qui m’a fi-
nancé durant ces quatre années en m’accordant une bourse de doctorat F.R.I.A.
ainsi qu’en co-financant un modeste cluster parallèle de visualisation.
De plus, je remercie chaleureusement le CENAERO (Center for Excellence
in Aeronautics) et dans une moindre mesure le C.I.S.M (Centre de calcul in-
tensif et de stockage de masse) qui m’ont permis l’utilisation de leur clus-
ter de calcul parallèle afin de réaliser les simulations numériques présentées
dans ce document. Leur support et leur connaissance technique m’ont été fort
bénéfiques.
Ces quatre années n’aurait pas été aussi agréables sans l’ambiance qui
régnait au sein du pôle de recherche. Merci à mes collègues : les Wincky-boys
(Matthieu, Ivan, Laurent, Roger, Goéric, Yves, Nicolas) bien entendu. . . mais
aussi tous les autres: Amel, Xavier, Julien, Francesco, Nicolas, Frank, . . .
Je remercie tout particulièrement Matthieu avec qui j’ai partagé mes débuts
en tant que “chercheur” et qui, au fil des années, est devenu, par son regard
avisé et ses remarques toujours pertinentes, un deuxième promoteur.
Cette thèse aurait été fort différente sans la deuxième vie apportée par
Philippe. La collaboration a été malheureusement trop courte mais néanmoins
vi
fructueuse. Il a été, de facto, mon troisième promoteur!

Merci à Dominique pour la relecture de ce document. Merci à celle qui
partage ma vie, Audrey, et qui m’a soutenu tout au long de ces années. Merci
à toute ma famille et mes amis.
Enfin, je souhaiterais remercier la communauté open source. Aujourd’hui,
un ensemble performant et gratuit de logiciels professionnels sont développés,
et ce toujours dans l’esprit originel! Merci à tous ces anonymes. . .
Contents
1 Introduction 1
1.1 Background . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1
1.2 Some preliminaries . . . . . . . . . . . . . . . . . . . . . . . . . 3
1.3 Objectives . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5
1.4 Outline . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6
2 An IB technique within the combined VPM-PFM method 9

2.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10
2.2 An Immersed Boundary Method Within a Vortex Particle-Mesh
Method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12
2.3 Validation: flow past a sphere at Re = 300 . . . . . . . . . . . . 20
2.4 Flow past a very low aspect-ratio wing . . . . . . . . . . . . . . 27
2.5 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 36
3 A mesh refinement technique for the VPM method 37

3.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37
3.2 A basic multi-block generator . . . . . . . . . . . . . . . . . . . 39
3.3 A redistribution scheme based on average-interpolation . . . . . 41
3.3.1 Multi-level redistribution scheme . . . . . . . . . . . . . 41
3.3.2 Treatment of varying resolutions . . . . . . . . . . . . . 44
3.4 A 1-D multiple resolution example . . . . . . . . . . . . . . . . 47
3.5 Vortex ring test case . . . . . . . . . . . . . . . . . . . . . . . . 49
3.6 Multiple resolution simulation of the flow past a sphere . . . . 62
3.7 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 67
viii Contents
4 Investigation of velocity deficit effects on aircraft wake rollup 69

4.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 69
4.2 Initial conditions . . . . . . . . . . . . . . . . . . . . . . . . . . 70
4.3 Numerical setup . . . . . . . . . . . . . . . . . . . . . . . . . . 78
4.4 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 81
4.4.1 3D dynamics: a qualitative analysis . . . . . . . . . . . 81
4.4.2 Quantitative analysis . . . . . . . . . . . . . . . . . . . . 82
4.4.3 Characterization of the rolled-up vortex sheet . . . . . . 103
4.5 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 110
5 Numerical investigations of end-effect 113

5.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 114
5.2 Numerical method . . . . . . . . . . . . . . . . . . . . . . . . . 117
5.3 Analysis and comparison . . . . . . . . . . . . . . . . . . . . . . 119
5.4 T-D and S-D simulations of an accelerated wing . . . . . . . . 132
5.5 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 159
6 Conclusions & perspectives 161

6.1 Achievements . . . . . . . . . . . . . . . . . . . . . . . . . . . . 161
6.2 Perspectives . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 163
6.2.1 A VPM method using a staggered grid . . . . . . . . . . 163
6.2.2 Further improvements of the Poisson solver for hierarchi-
cally refined grids . . . . . . . . . . . . . . . . . . . . . . 164
6.2.3 Improvement of the handling of no-slip . . . . . . . . . . 165
6.2.4 A coupling approach using an Eulerian “near-wall” solver 165
Bibliography 170
A S-D simulation of a near-wake rollup using a modeled inflow 181
B 3D filament simulations of Crow-like instabilities IGE 187

B.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 187
B.2 Vortex filament method . . . . . . . . . . . . . . . . . . . . . . 188
B.3 Definition . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 190
B.4 Comparison: numerical simulation versus experiment . . . . . . 191
B.5 Analysis of the Crow instability in ground effect . . . . . . . . 194
Chapter 1
Introduction
1.1 Background
The present work joins a global effort in the development of efficient and accu-
rate numerical tools for the simulation of complex fluid flows. More specifically,
we focus on unsteady external flows, a fluid dynamics discipline which actu-
ally pervades applied sciences and engineering. These are the flow encountered
in aircraft or car aerodynamics, wind energy, biological locomotion, etc. In
the problems addressed here, we assume that the flow is incompressible and
that the inertial forces dominate the viscous stresses: we are dealing with both
moderate and high Reynolds number flows.
The governing equations of fluid mechanics were derived between the 15th
and 19th century, culminating with the works of Navier and Stokes some 150
years ago. The resulting Navier-Stokes equations have since then been one of
the pillars of fluid mechanics and have been used universally in all fields of
science.
The development of computers, since the seventies, has enabled new per-
spectives for the study of complex fluid flow problems. This discipline, dedi-
cated to the numerical resolution of the fluid mechanics equations, is called com-
putational fluid dynamics (CFD). The straightforward solution approach con-
sists in capturing all scales of the flow through a very fine space-discretization;
a direct numerical simulation (DNS) is then performed. This approach however
quickly leads in practice to an unaffordable computational time.
2 Chapter 1. Introduction
Two main approaches are commonly used to tackle this issue: the Large
Eddy Simulations (LES) and the Reynolds Averaged Navier-Stokes (RANS)
simulations. The former models the effect of the unresolved scales on the re-
solved scales. The latter directly obtains the time-averaged solution through
the modeling of the Reynolds stress tensor. Because we are interested in ac-
curately capturing transient or unsteady phenomena, we here rely on the first
approach, LES.
Furthermore, the present work treats the flow equations in a very distinct
approach. We use a vortex particle-mesh (VPM) method, which belongs to
the broader class of âvortex methodsâ. Such methods use the vorticity-velocity
formulation of the Navier-Stokes equations, rather than the velocity-pressure
formulation: the vorticity field is thus the primary variable. The major advan-
tage comes from the compactness of the vorticity field for external flows and
wakes. A limited number of particles is thus required to discretize the entire
flow.
In the vortex method used here, the elements of discretization are “vortex
particles” that are convected with the local velocity vector. The major advan-
tage comes from the compactness of the vorticity field for some flows of interest
(e.g., external flows and wakes). A limited number of particles is then required
to discretized the entire flow. Moreover, vortex methods intrinsically take into
account unbounded-flow boundary conditions.
Vortex particle methods, however, have drawbacks, some of them funda-
mental. First, the element of discretization being isotropic, they are not in-
trinsically adapted to capture physical phenomena such as boundary layers.
Second, they require, in general, more computational time. Finally, with time,
the particles can cluster or deplete in some regions. Redistribution schemes
with high order conservation properties are then required; they allow perform-
ing accurate long time simulations avoiding the clustering/depletion of particles
in inherent to any Lagrangian method.
Vortex methods have been developed for quite some time and are now ap-
plied to aerodynamics, fundamental fluid mechanics or biological locomotion
(Hieber and Koumoutsakos [36], Zhang and Eldredge [93]). They are used for
incompressible and compressible flows (Eldredge [28]), reacting flows (Thiri-
fay [82]), viscous and inviscid flows, free vortical flows or wall-bounded flows
1.2. Some preliminaries 3
(Koumoutsakos [47], Ploumhans and Winckelmans [67] and Cottet and Poncet
[17]).
DNS and high quality LES can also be performed. Cocle et al. [12] proposed
an accurate multi-scale subgrid model for such LES.
1.2 Some preliminaries

Before going further, let us quickly go through the fundamentals of vortex
methods.
Formulation equations The Navier-Stokes equations for an incompressible

flow, thus ∇ · u = 0, in vorticity-velocity formulation read
Dω
= (∇u) · ω + ν∇2 ω, (1.1)
Dt
where u is the velocity field and ω , ∇ × u is the vorticity field, thus ∇ · ω = 0.

D ∂
Dt , ∂t + u · ∇ is the material derivative. Using the Helmholtz decomposition,
one can write the velocity field as
u = ∇ × Ψ + ∇φ, (1.2)
where Ψ is the streamfunction and φ is the scalar potential. We take φ such

that ∇φ = U0 , where U0 is the upstream velocity. If we set the gauge of Ψ as
∇ · Ψ = 0, the vorticity field and the stream function are related by the Poisson
equation
∇2 Ψ = −ω. (1.3)
Discretization The vorticity field is discretized using vortex particles which

are responsible for a volume Vp and carry a vorticity vector quantity (αp ,
R
α(xp ) = Vp ωdV ≃ ω(xp )Vp ). The evolution equations for the position and
the strength of those particles read
dωp
= (∇u(xp )) · ω(xp ) + ν∇2 ω(xp ) (1.4)
dt
dxp
= u(xp ) ∀p = [1 . . . N ] (1.5)
dt
Computation Vortex particle methods require to solve Eq. (1.2) and (1.3)
to obtain the velocity field from the particle strengths. For unbounded do-
mains, Eq. (1.2) and (1.3) can be solved using fast summation techniques for
the Green’s function, such as the Parallel Fast Multipole (PFM) method to
reduce cost complexity from O(N 2 ) to O(N log N ) (Barnes and Hut [1], Green-
gard [33] and Greengard and Rokhlin [34]). The parallel implementation used
here has been developed by Salmon and Warren [74] and further extended for
vortex methods by Winckelmans [90]. However these methods are still expen-
sive for the simulation of external flows, even at moderate Reynolds number
(e.g. O(103 )).
The presents approach tackles this problem through the efficient combina-
tion of the “Vortex Particle-Mesh” (VPM) method, also called “Vortex-In-Cell”
(VIC), originally developed by Christiansen [10], see also Cottet and Poncet [17]
and Cocle et al. [13]) method and the PFM method. The VPM methods solve
the Poisson equation (Eq. (1.3)) using fast grid solvers that are much faster
than the fast summation for the Green’s function based Poisson solvers. The
required boundary conditions are here provided by the PFM method. Those
boundary conditions being accurate, a compact grid can be used, that is tight
to the vorticity field. This is the basis of the hybrid VPM-PFM method. From
a parallel computational point of view, we note that there is no need to solve
a global Poisson equation. Rather, multiple Poisson equations can be solved
locally, on each subgrid, using the proper boundary conditions provided by the
PFM method (which has a global view of the whole vorticity field).
These hybrid methods require the source term of Eq. (1.3) to be available
on the grid. This entails the translation of particles strength into a grid vor-
ticity through a redistribution. The evaluation of the right hand side (RHS)
of Eq. (1.4) and Eq. (1.5) are required to integrate in time the evolving equa-
tions. Those are here computed using finite differences. That requires, on the
one hand, to interpolate particle strengths onto the grid in order to compute
the time-variation of each particle intensity and position. On the other hand,
this “grid information” has to be interpolated from the grid to the particles.
This method has been successfully developed by Cocle et al. [13, 11] for (half-
)unbounded vortical flows.
We also note that the Poisson equation can also be solved with unbounded
1.3. Objectives 5
boundary conditions using the Hockney algorithm [37] as done for 2-D wall-
bounded flows by Morgenthal and Walter [61] and for 3-D free vortical flows
by Chatelain and Koumoutsakos [9]. This approach is very efficient because it
uses fast Fourier transforms (see [8] for a massively parallel application using a
billion grid points), but it requires twice more grid points for each unbounded
direction. The introduction of a multi-resolution capability in the latter ap-
proach is likely quite complex.
1.3 Objectives
It is now well established that vortex method definitely overcome classical meth-
ods (i.e. CFD solvers based on the discretization of the Navier-Stokes equations
using the velocity-pressure formulation) when simulating free vortical flows (see
Winckelmans et al. [88], Daeninck et al. [22, 21], Cocle [11], Cottin et al. [19]).
In addition, this is still true whatever the value of the Reynolds number.
In spite of the major developments in vortex method in the last decades, those
methods rapidly become unaffordable when carrying out the simulation of vor-
tical flows interacting with a viscous ground or the simulation of external flows
past bodies. It is still very challenging to perform external flow simulations
using vortex methods, especially when the Reynolds number is high. They
are intrinsically well adapted to capture varying flow resolutions, yet they still
require further improvements in order to benefit from the recent advances:
• Combining a fast grid-based Poisson solver together with the handling of

a varying resolutions. External flows are generally characterized by the
need for multiple resolutions.
• Carrying out the flow past moving and deforming bodies. Vortex methods
require to deal with discontinuous vorticity field, as we here consider
the Navier-stokes equations in vorticity-velocity formulation: the velocity
tends to zero at the wall whereas the vorticity is maximum at the wall
and zero inside (for non rotating bodies), and thus discontinuous.
• Efficiently dealing with boundary layers. Vortex particle methods are

intrinsically not adapted to capture boundary layers as the element of
discretization (i.e., a particle) is isotropic.
The present research work aims at pursuing the development of the vortex
particle-mesh (VPM) method with two outstanding contributions in order to
address those problems.
First an immersed boundary technique has been developed and adapted
to capture the flow past arbitrary shape bodies. The method is based on
the approach first developed by Koumoutsakos [45] and further improved by
Ploumhans [67, 68]. The motivation is to benefit from the enhanced efficiency
provided by the VPM method in terms of computational cost. The required
Poisson equation is solved using an efficient grid-based solver combined with a
parallel fast multipole method (which provides the required boundary condi-
tions on each subdomain) as shown in previous works, also within our group,
see Cocle et al. [11]. The method is here further developed in the spirit of the
immersed boundary method as we model an extra term of the evolving equa-
tions, expressed as grid quantities, in order to take into account the presence
of a no-slip boundary.
Second, an accurate approach handling hierarchically refined meshes has
been developed and validated. The basis of this technique is in the same
spirit as the methodology developed by Bergdorf et al. [2] and Hejazialhosseini
et al. [35]. We use both grid patches and particles of varying resolution.
The originality of this work is the combination of the handling of multiple
flow resolutions together with a full 3-D VPM Navier-Stokes solver to com-
pute incompressible flows. This method also benefits from the versatility of
the parallel fast multipole method which enables the efficient solution of the
Poisson equation and straightforward domain decomposition. We have also
demonstrated the potential of the mesh refinement technique on a Large-Eddy
Simulation of a turbulent vortex wake rollup at very high Reynolds number.
1.4 Outline
The first part of the manuscript provides the description of two advances in
vortex methods and then their application to the numerical investigations of
bluff-body and wake-vortex flows. The two major improvements done in this
work are: i) using the VPM-PFM method with the Immersed Boundary (IB)
method to compute wall-bounded external flows at a moderate Reynolds num-
1.4. Outline 7
ber and ii) including a multi-resolution approach in the VPM-PFM method.
• Chapter 2 presents the numerical algorithm developed to compute the

flow past arbitrarily-shaped geometries (lifting- and bluff-bodies) at a
moderate Reynolds number (i.e. O(102 − 103 )). The methods is based on
an immersed boundary approach to solve the Poisson equation required
to obtain the streamfunction. The vorticity sources, corresponding to
the cancellation of the residual slip velocity and needed by the Poisson
equation, are here computed by means of the vortex panel method. This
approach is validated against the benchmark of the flow past a sphere at
U0 c
Reynolds 300. The method is also used for the DNS (at Re = ν = 3000)
of an impulsively started flow past a profiled body: i.e. a very low aspect
ratio wing.
• Chapter 3 deals with the integration of a multiple space-resolution tech-

nique within the VPM-PFM method. The approach uses the classical
redistribution schemes and the average-interpolation wavelets. It pro-
vides smooth functions across grids of different resolutions. In practice,
it turns out that the linear impulse is quite well conserved. To make
the code more efficient in terms of load balancing (computational cost
and memory cost), both the parallelization and the mesh management
have been redesigned. The method is successfully tested on a vortex ring
simulation using various mesh size transitions. The chapter also focuses
on the gain resulting from the combination of the immersed boundary
approach and refined grids. The simulation of the flow past a sphere at
Reynolds 300 is carried out. Multiple resolutions are used to capture
the wake up to 38 diameters. The computational gain is highlighted in
comparison to the results obtained in Chapter 1.
The second part of the manuscript is devoted to the numerical investigation of

various vortex-systems using vortex methods as efficient numerical investigation
tools, i.e., the combined VPM-PFM (called VIC-PFM in those chapters as they
consist of previous published reports where the term VIC was used instead of
VPM) and the regularized vortex filament method. The focus is on the physics
of the investigated problems.
• Chapter 4 investigates the effect of the velocity deficit due to the fric-
tion drag on the near-wake rollup of a vortex-sheet at Re = 106 up to
10 wingspan downstream of a lifting body. Space-developing LES, with
and without velocity deficit, have been carried out. First, the space-
developing rolled-up vortex pair with friction drag is analyzed and char-
acterized (circulation and tangential velocity profiles, vortex-core radius,
residual inner-core velocity, turbulence, etc.). Second, the two simula-
tions are compared in terms of net-result after rollup. The effect of the
velocity deficit effect is also analyzed in the case of a time-developing
simulation.
• Chapter 5 investigates the “end-effect” phenomena. The end-effect oc-

curs when the total circulation of the vortex system varies longitudinally.
The associated pressure wave induces a longitudinal velocity in the vortex
core. The system possibly becomes unstable, leading to the creation of a
helical instability. A comparison to the experiment of a counter-rotating
vortex-pair (carried out by the CNRS-IRPHE) is done. Moreover the
space-developing end-effect is investigated in the case of an accelerated
vortex sheet.
Finally, the effect of an inviscid ground on the Crow(-like) instability growth is

considered in Appendix B. The regularized vortex filament method was used.
Two cases were investigated. First, a comparison to an experiment (carried out
by CNR-IRPHE) is done. The simulation replicates as much as possible the
experimental setup, that is, the evolution of the long-wave instability of a pair
of fat vortices put in ground effect, and as long as viscous effect are negligible.
Second, the influence of various initial perturbation levels of a thin vortex pair
is analyzed.
Chapter 2
Development of an
immersed boundary
technique within the
combined Vortex
Particle-Mesh/Parallel Fast
Multipole Method
This chapter was written for the ECCOMAS 2010 conference proceeding as a
paper entitled “Development of an immersed boundary method using boundary
elements within a vortex-in-cell/parallel fast multipole method”. This article
is essentially reproduced as it was, yet with some modifications and with some
additional results.
Abstract This paper presents a method for the simulation of incompress-

ible flows past arbitrary shape bodies. The method couples a combined vortex
particle-mesh (VPM) and parallel fast multipole (PFM) method to an Immersed
10 Chapter 2. An IB technique within the combined VPM-PFM method
Boundary (IB) approach. The Navier-Stokes equations in the vorticity-velocity

formulation are solved. The Poisson equation for the streamfunction is solved
efficiently in an unbounded domain by means of a direct solver; in each subdo-
main, the boundary conditions are provided by the PFM method which has a
global view of the whole vorticity field.
Because the vorticity field is a compact field, a tight grid can be used. This
makes an efficient method for unbounded vortical flows.
An IB method is used in order to deal with complex body geometries. In
vorticity-velocity formulation, imposing the no-slip condition requires to solve
a vector integral equation. This is done here using a boundary element method.
The obtained panel strengths provide the equivalent vorticity flux at the wall.
The method has been successfully validated using a well-referenced test case,
which consists in the flow past a sphere at Re = 300. It is then also tested on
the impulsively started flow past a very low aspect ratio wing at Re = 3000.
2.1 Introduction
This paper describes the extension of the VPM-PFM method in order to com-
pute external flows past arbitrary geometries at moderate Reynolds number.
The methodology is based on an Immersed Boundary (IB) approach. This
was first developed by R. Peskin [65] for velocity-pressure based Navier-Stokes
solvers for the simulation of inner-heart flows. See also the review of Mittal
et al. [58] for more details on IB methods.
Several previous works have dealt with combining a vortex method to an
immersed boundary technique. Morgenthal et al. [61] have also developed a
2-D IB method for VPM algorithms. An influence matrix (Particle-Particle
Particle-Mesh algorithm, the P 3 M ) is used in order to capture the under-
resolved field in the near-wall region. Cottet and Poncet [17], and Poncet [70]
developed an IB method in a VPM code using potential sources instead of
vorticity sources. In their approach, the body surface and the potential source
are discretized over a few grid points. These sources are solved at every step.
Brinkman penalization methods have also been used to capture such flows.
(See Coquerelle and Cottet [15] and Rossinelli et al. [72])
We here capture the discontinuities of the velocity field, due to slip, by
2.1. Introduction 11
means of a panel method, thus in a step reminiscent of a sharp interface method

(see [85, 54, 52, 92]). In a second step however, this jump, corresponding to
a vorticity flux has to be passed to the grid. Our approach is thus indeed an
immersed boundary method as the wall discretization is finally mollified.
Vortex methods solve the Navier-Stokes equations in the vorticity-velocity

formulation. The vorticity field is discretized using N vortex particles; each
carrying a strength vector and moving with the local velocity vector. The
convection is thus done in a Lagrangian manner; leading to negligible dispersion
error. The evolving equations for the strength and the position of a vortex
particle p read
dωp
= (∇u(xp )) · ω(xp ) + ν∇2 ω(xp ) (2.1)
dt
dxp
= u(xp ) ∀p = [1 . . . N ] . (2.2)
dt
The velocity field u is computed from the vorticity field ω through the stream-
function Ψ by solving the following Poisson equation
∇2 Ψ = −ω. (2.3)
The velocity is then computed using the relation
u = ∇ × Ψ + U0 . (2.4)
This method entails to “translate” the particle strength onto a grid strength.
This is here done using the redistribution procedure. The redistribution is also
used to handle Lagrangian distortion in particle methods. With time, particles
tend to cluster or deplete some region of the flow. Redistribution proceeds by
periodically generating a new set of particles (from the old set) on a regular
lattice.
This chapter is structured as followed. First, the algorithm is described in

details. Secondly, the proposed methodology is validated against the flow past
a sphere at Re = 300; a well-referenced test case in the literature. Finally, the
method is applied to the impulsively started flow past a profiled-body consisting
in a low aspect ratio NACA0012 wing.
2.2 An Immersed Boundary Method Within a

Vortex Particle-Mesh Method
The coupling of an IB method is well-suited for codes using the velocity-pressure
formulation because velocity components tend to zero at the wall. An IB
method in a vorticity-based Navier-Stokes solver however has to handle a dis-
continuous field. Indeed, the vorticity is non zero on the fluid side and maxi-
mum at the wall and it is zero inside the body up to the wall. In the present
approach, the wall treatment relies on a symmetrized vorticity field with re-
spect to the wall. The presented algorithm is mainly based on the approach
developed by Ploumhans [66, 68] and adapted to the VPM-PFM method [13].
Figure 2.1: General setup of the IB method.
The algorithm reads: At time tn , the vorticity field is discretized using a

set P1 of particles characterized by their position and their strength: (xnp , αnp )
for p ∈ [1 · · · N ] where N is the number of particles. Moreover, the arbitrary
body shape is discretized using a set S of flat triangular panels forming a closed
surface and described by the position of the three vertices: (xn1,j , xn2,j , xn3,j ) for
j ∈ [1 · · · M ] where M is the number of panels. The following 10 steps are
performed to obtain the solution at time tn+1 .
1. Mesh overlay & Random grid shifting. A Cartesian grid with a

2.2. An Immersed Boundary Method Within a Vortex Particle-Mesh Method13
uniform resolution h, is laid over the support of vorticity, i.e. it covers

the location of all particles and panels.The particle/panel position are
randomly shifted between −0.5h and 0.5h, in each direction, in order to
distribute the IB discretization error over the body surface. Indeed, the
body discretization crosses arbitrarily the Cartesian grid.
2. Wall distance. The grid points are tagged relative to the wall distance.
For each grid point xijk = (xi , yj , zk ) = (x0 + ih, y0 + jh, z0 + kh), a point
xwall on the wall surface exists such that d = min(|xijk −xwall |), d is then
the distance from xijk to the wall. In practice, the distance to the wall d
is here computed as done in Daeninck et al. [21]: first the nearest panel
to the point xijk is computed by means of the nearest neighbor search
functionality of the fast multipole panel method. Second, the distance
is defined as the normal component in the local reference frame of the
nearest panel.
3. Evaluation of the vorticity field on the grid. We construct a set P2

of image vortex particles inside the body such that the spurious vorticity
flux at the wall, corresponding to the case without those particles, is
cancelled. Of course, there is always a small spurious effect due to the
non-zero local curvature, R, of the wall, as compared to h (need R/h ≫
1). The methodology is sketched in Fig. 2.2. For each near-wall particle,
an equivalent image particle is created inner the body which has the
same tangential vorticity components and an opposite normal vorticity
component. As we use a “halo” region (see Ploumhans [66, 68]), the
union of the sets P1 and P2 does not cover the entire flow and a small
gap exists around the wall (i.e. d/h < 0.25 where d is the distance to the
wall). For these near-wall grid points, using the particle redistribution
scheme provides an “incomplete” field. The near-wall particle weight thus
has to be extrapolated from the flow along the normal to the wall. This
is here done using one-side second-degree polynomial functions P:





 Pt,n (xwall + hn) = αt,n (xwall + hn)




 Pt,n (xwall + 2hn) = αt,n (xwall + 2hn)
(2.5)
 dPt
dn (xwall ) = 0 for tangential components







 Pn (xwall ) = 0 for normal component
where Pn = ∇P · n and Pt = P − Pn n.
Figure 2.2: General description of the methodology for computing the near-wall
vorticity field in the local reference frame. The “bullets” (•) and the “stars” (∗)
respectively represent the set P1 of outer-body particles and the set P2 of inner-body
image particles. To compute the vorticity field at the near-wall grid point in the “halo”
region (here shaded region), one-sided extrapolating functions are used (differentiating
the normal an the tangential components).
Afterward, the weights of the particle sets (P1 and P2 ) are redistributed
′
onto the complete grid (no halo region) using the M4 scheme:
′
M4
αnp −−−−−→ αnijk . (2.6)
More precisely,
N
X ′ ′ ′
αnijk = M4 (xi − xp ) M4 (yj − yp ) M4 (zk − zp ) αnp . (2.7)
p=0
where, 



 1 − 25 |x|2 + 32 |x|3 , |x| ≤ 1,
′

M4 (x) = 1 2 (2.8)
 2 (2 − |x|) (1 − |x|), 1 < |x| < 2,



 0, |x| ≥ 0.
′
The M4 scheme has good conservation properties, see Monaghan [60],
Cotted and Koumoutsakos [16] and Winckelmans [87]. The vorticity, the
linear impulse and the angular impulse are conserved.
4. No-slip condition. The vortex panel vector intensity is computed such

that the tangential vector velocity underneath each panel is canceled (see
Fig. 2.1).
Imposing the no-slip flow condition at the wall requires to solve a vector
integral equation. This is done here using a Boundary Element Method
(BEM). For vortex panels of vector strength ∆γ:
Z
1 1
½∆γ(x) × n + · (x − x′ ) × ∆γ dS(x′ ) = U0 + uω (x)
4π S |x − x′ |3
= ∆Uslip (2.9)
where uω is the solution of Eq. (1.2) and (1.3), the vortex-field induced ve-
locity. This equation is solved using an iterative method (a under-relaxed
Jacobi type algorithm). The surface body is discretized using triangles.
The matrix-vector multiplications (i.e. particle-panel and panel-panel in-
teractions) in the linear solver are evaluated efficiently by means of the
PFM method.
The RHS, i.e. the residual velocity ∆Uslip , is also evaluated using the
PFM method. The particles are considered for this step as “point” par-
ticles so as to have an entire view of the vorticity field. A “halo” re-
gion, where no particle can exist, is defined to avoid any spurious vor-
ticity flux at the wall see Ploumhans et al. [67, 68]. Typically, one fixes
dhalo /h = 0.20 · · · 0.25. Particles entering this zone are rebounded with
respect to the wall.
5. Vorticity flux. A Lighthill’s model of vorticity flux at the wall is

used [16, 67, 68]. For each grid point xijk , the net circulation increase
∆αnijk,cons , due to the diffusion of the vortex panels within the flow be-
tween tn and tn+1 , is computed using a conservative scheme. Let’s con-
sider a panel with an area Ω and an intensity ∆γ. Defining an infinites-
imal surface dS = dξ dη, the net circulation increase over a time ∆t is
given by the following formula
Z 1
dαijk
∆αijk = d(t/∆t),
0 d(t/∆t)
where
√
dαijk 1 (z −hl /2)/ 4νt
= ∆γ [erfc(u)](zk +h/2)/ √
d(t/∆t) 2 k 4νt
Z √ √
(x −h/2−x(ξ,η))/ 4νt (y −h/2−y(ξ,η))/ 4νt
[erfc(u)](xi +h/2−x(ξ,η))/√4νt [erfc(u)](yj +h/2−y(ξ,η))/√4νt dξ dη .
i j | {z }
S
dS
(2.10)
ξ and η are the local integration coordinates. hl /2 is equal to zk if

zk < h/2 + dhalo and equal to h/2 otherwise. This formula is here nu-
merically integrated i) in time using four-point Gauss quadrature and ii)
over the triangle using the three-point Hammer rule. This formula can be
analytically integrated if the panel is a rectangle (S = f × b). Eq. (2.10)
can then be rewritten as followed (see Ploumhans [66, 68])
√
dαijk (z −hl /2)/ 4νt
= νt∆γ [erfc(u)](zk +h/2)/ √
4νt
d(t/∆t) k
√ √
(x +h/2−b/2)/ 4νt (x +h/2+b/2)/ 4νt
[ierfc(u)](xi −h/2−b/2)/√4νt − [ierfc(u)](xi −h/2+b/2)/√4νt
i i
√ √
(yj +h/2−f /2)/ 4νt (yj +h/2+f /2)/ 4νt
[ierfc(u)](y −h/2−f /2)/√4νt − [ierfc(u)](y −h/2+f /2)/√4νt .
j j
(2.11)
In practice, if the triangular panels are equilateral enough, the panels

√
can be considered as rectangular (i.e. f , b , S△ ). Otherwise, the
integral is done numerically. In the present work, both are available and
needed (as triangles are not always “almost equilateral”) such as in the
discretization of an airfoil, see later. Moreover, the scheme can be made
conservative through a weighted correction as developed by Ploumhans
et al. [68],
!
|∆αijk |2 X
∆αwall
ijk,cons = ∆αijk +P 2
S∆γ − ∆α(r,s,t) .
r,s,t |∆α(r,s,t) | r,s,t
6. Poisson equation. Solving the Poisson equation ∇2 Ψnijk = −(ωijk

n
+
n
∆ωijk ) using the mathematical library FISHPACK (for more details con-
cerning this library, see [80, 78, 79]). ∆ωi,j,k is the “representation” of the
panel vorticity onto the grid. It is computed either using the redistribu-
tion scheme or using the diffused-panel field. The former provides robust
behavior for impulsively started flows and the latter provides better re-
sults in regime mode. The required boundary conditions are obtained
using the PFM method.
7. Evolution equation. The RHS of Eq.(2.1) and (2.2) are computed

using Finite Differences (FD)
unijk = ∇ × Ψnijk
and
n
dωijk
= (∇unijk ).ωijk
n
+ ν∇2 ωijk
n
.
dt
Notice that the velocity and the stretching terms are computed using
centered second-order FD. The diffusion term is computed using a 27
points “cube” scheme instead of the “cross” scheme (see [13]), which
provides less sensitivity to the grid orientation.
8. Vorticity divergence. With time, the vorticity field does not remain
an image of the velocity curl for several reasons. i) The interpolation
′
procedure (particle-mesh, mesh-particle) using the M4 scheme does not
conserve the divergence of the vorticity field. ii) The stream function Ψ
obtained by solving the Poisson equation, using FISHPACK, is not diver-
gence free. Moreover, a “collocated” lattice is used instead of a “stag-
gered” lattice. The numerical scheme is thus not fully consistent: the
numerical Laplacian operator does not match with the numerical double
rotational operator (∇h × (∇h × Ψ) 6= −∇2h Ψ = ω).
The vorticity field thus has to be regularly corrected in order to project it

into a divergence free space, following a methodology first developed by
Cotted and Poncet [17], and Poncet [69]. Such a methodology was also
used by Cocle et al. [13]
The vorticity field ω can be corrected to a divergence free field ω̃:
ω̃ = ω − ∇D, (2.12)
by taking D such that

∇2 D = ∇ · ω. (2.13)
Moreover, the presence of a wall entails boundary conditions for the above
equations:
ω̃ · n = 0 ⇔ ∇D · n = ω · n. (2.14)
The Poisson equation is solved using the same combination of the FISHPACK
library and the PFM method. For the same reason that the vorticity
field is not exactly divergence free (consistency between Laplacian op-
erator versus double rotational operator), the above correction scheme
(Eq. (2.13) and (2.12)) does not exactly provide a solenoidal field. How-
ever, it turns out that this step provides satisfactory results, if the vor-
ticity field is smooth enough.
However the near-wall vorticity field is often locally weakly-resolved, due
to the “halo” region. The evaluation of the divergence sources then leads
to spurious corrections of the vorticity field. For this reason, the near-
wall divergence sources (i.e. d < 2h) and the boundary condition are not
taken into account during this correction step.
9. Interpolation. Grid quantities (i.e. velocity, time-derivative of strength

and net increase of intensity due to panel diffusion) are interpolated onto
′
the particles using the M4 scheme
′
M4
unijk −−−−−→ unp
n ′ n
dα M4 dα
−−−−−→
dt ijk dt p
′
M4
∆αn,wall
ijk −−−−−→ ∆αn,wall
p
10. Time integration. The particle positions and strengths are integrated
in time using respectively the second-order Leap-Frog (LF) scheme and
the second-order Adams-Bashford (AB) scheme
xn+1
p = xpn−1 + 2∆tunp .
n−1 n !
1 dα dα
αn+1
p = n
αp + ∆t 3 − + ∆αn,wall
p
2 dt p dt p
tn+1 = tn + ∆t.
After particle redistribution (typically every 5 time steps), the AB in-

tegration scheme cannot be used. It is replaced by the second-order
Runge-Kutta (RK) scheme. then the first step (predictor) is
x∗p = xnp + ∆tunp

" n #
dα ∆αn,wall
p
α∗p n
= αp + ∆t +
dt p ∆t
t∗ = tn + ∆t,
and the second step (corrector) is
1
xn+1
p = xnp + ∆t unp + u∗p
2 " #
n ∗
1 dα ∆αn,wall
p dα ∆α∗,wall
p
αn+1
p = n
αp + ∆t + + +
2 dt p ∆t dt p ∆t
tn+1 = t∗ .
2.3 Validation: flow past a sphere at Re = 300

The present VPM-PFM-IB methodology has been validated on a well-known
test case for which experimental and numerical data are available: the flow past
U0 D
a sphere at Re = ν = 300. At such a Reynolds number, the flow is unsteady,
thus making it a relevant test case, even though at moderate Reynolds number.
Figure 2.3: General view of the flow past an arbitrarily shape.
The force acting on a body can be computed from an exact formula devel-
oped by Noca et al. [62, 63] and applied to PIV (Particle Image Velocitmetry)
measurements. This formula was also successfully used in the context of vortex
particle method by Ploumhans et al. [68] and Daeninck [21]. Let us consider a
control volume VC with a surface SC which includes the body described by the
simply-connected closed surface Sb , see Fig. 2.3. The force acting on S reads
Z
F 1 d
= − ((u · x)n − u(x · n) + (N − 1)x(u · n)) dS
ρ N − 1 dt SC +Sb
Z
1
+ (u · u)n − (n · u)u dS
SC 2
Z Z
1 1
− (n · u)(x × ω)dS + (n · ω)(x × u)dS
N − 1 SC N − 1 SC
Z Z
1
+ x × (n × ∇ · T)dS + n · TdS, (2.15)
N − 1 SC SC

where T = ν ∇u + (∇u)T and N is the number of dimension (here N = 3).
2.3. Validation: flow past a sphere at Re = 300 21
This formula is particularly well adapted for vortex methods, because it only
requires the knowledge of the velocity and its derivatives; the pressure field is
not required. Note that the surface integral on Sb vanishes when u = 0 at the
wall (as for the flows considered in this thesis). The evaluation of the various
fields are thus only needed on the control surface SC . Here, the control volume
is a cube centered at the sphere center with the edge size fixed to 1.6D. The
size of the control volume has no effect, as shown in Daeninck [21].
The mesh resolution is here fixed to h/D = 1/75. An adaptive time step
is used such that both the maximal CFL criteria are |u|max ∆t/h ≤ 0.90 and
|ω|max ∆t ≤ 0.75. The resulting time step is ∆tU0 /D ≃ 8.4 10−3. The wake
is simulated up to 20 diameters downstream. The sphere is discretized using
20 480 quasi-equilateral triangles obtained by recursive subdivision of an icosa-
p
hedron, see Fig. 2.4. The average surface of the triangles is S△ /h = 0.93.
The flow was discretized using ∼ 100 106 grid points and ∼ 30 106 particles.
To perform the complete simulation, 128 processors were required for 3 weeks
in order to accumulate the time-history of the force.
Fig. 2.5 shows various views of the wake behind the sphere. The vortical
structure, the so-called “hairpin” vortices, are clearly present (here visualized
using the λ2 criterion [41]. Fig. 2.6 shows a 2-D cut of the vorticity field for
the entire flow (a) and a zoom in the near-wall region (b).
Figure 2.4: Sphere discretization using recursive subdivision of an icosahedron.
One defines the Drag coefficient (CD )
F · ey
CD = 1 2 2
(2.16)
2 ρU 0 (πD /4)
where ey is the streamwise unit vector. The Lift coefficient (CL ) is defined as
||F − (F · ey )ey ||
CL = 1 2 2
(2.17)
2 ρU0 (πD /4)
The time-evolution of the diagnostics (drag/lift coefficients) are presented in

Fig. 2.7. Table 2.1 shows a comparison with reference results. The time-
averaged value of CD and CL , as well as the Strouhal number(Str = f D/U0 ),
compare well with those of the literature. As for the oscillation amplitude
1 1
(∆CD , 2 (max CD − min CD ) and ∆CL , 2 (max CL − min CL )), they also
compare quite well with those of the literature. It is worth noting that those
latter values are more difficult to obtain. Moreover, results with higher quality
are obtained when using the highest allowed time step. As a matter of fact,
the larger the time step, the larger the number of grid points on which the
diffusion process from the panels is captured. For a given a resolution, it thus
exists a minimum time step. The ratio of effective diffusion length to the grid
size can be estimated as
√
ldif f 4ν∆t
≃ ndif f .
h h
The effective diffusion should be such that it covers at least ndif f grid points for
the diffusion to be well discretized. For a given Reynolds number (Re = U0 D/ν)
and a set value for ndif f (here we use ndif f = 4), the previous equation leads
to 2
U0 1 h
∆t ≃ Re .
D 4 D
It turns out that using a time step too small leads to oscillation of the force
diagnostics due to an inconsistency between the solution of the Eq. (2.11) (that
should cancel the residual slip velocity) and what the grid really receives after
the diffusion process.
During the simulation, the mesh Reynolds numbers, based on the vorticity
|ω|max h2 |u|max h
field (Re ω
mesh = ν ) and based on the velocity field (Re u
mesh = ν ),
were both maintained at a low level, of O(3 · · · 7) (this ensuring a good resolu-
tion of the flow physics), see Fig. 2.8.(a). The divergence of the vorticity field
was also maintained at an acceptable low level, of O(5 10−4 ), see Fig. 2.8.(b).
(a)
(b)
(c)
Figure 2.5: Isocontour of the λ2 criterion: (a) top view, (b) side view and (c) offset
view.
Table 2.1: Summary table for quantities of interest compared to the literature for the
flow past a sphere at Reynolds 300.
CD ∆CD CL ∆CL Str

−3 −2
Present work 0.672 3.4 10 0.065 1.5 10 0.134
−3 −2
Johnson & Patel [42] 0.656 3.5 10 0.069 1.6 10 0.137
Tomboulides & Orszag [83] 0.671 2.8 10−3 - - 0.137
Constantinescu et al. [14] 0.655 - 0.065 - 0.136
Ploumhans et al. [68] 0.683 2.5 10−3 0.061 1.4 10−2 0.135
Georges [31] 0.661 2.5 10−3 0.066 1.3 10−2 0.134
(a)
(b)
Figure 2.6: 2D slice of the z-component of the vorticity field UD0 ω: (a) global view
and (b) zoom in the near-wall region. Saturation was used in order to see both the
near wall vorticity and the wake vorticity.
0.7
0.6
Force Coefficient [-]
0.5
0.4
0.3
0.2
0.1
0
110 120 130 140 150 160 170 180 190
tU0 /D
Figure 2.7: Time-evolution of the force acting on the sphere: drag coefficient (solid
line) and lift coefficient (dashed line).
(a)
10
6
Re mesh
0
110 120 130 140 150 160 170 180 190
tU0 /D
(b)
−3
x 10
1
0.9
0.8
|(∇ · ω)|dV
0.7
0.6
0.5
Vω
R
0.4
U0 Vω
D2 1
0.3
0.2
0.1
0
110 120 130 140 150 160 170 180 190
tU0 /D
Figure 2.8: (a) Time-evolution of the mesh Reynolds number Re ω mesh (solid line)
and Re u
mesh (dashed line). (b) Time-evolution of the divergence of the vorticity field.
Vω is the volume where ω is non zero.
2.4. Flow past a very low aspect-ratio wing 27
2.4 Application: flow past a very low aspect-

ratio wing at Re = 3000
This section describes the results obtained when computing the flow past a
lifting body such as a rectangular wing at moderate Reynolds number.
Initial condition & numerical setup The IB technique within the VPM-
PFM method has been also applied to the flow past a profiled body (i.e. a
“wing”) at a moderate Reynolds number Re = cU0 /ν = 3000, where c is
the wing chord. More precisely, we here consider a body consisting in a very
low aspect-ratio rectangular wing (AR , b/c = 1) using a NACA0012 profile
section. The angle of attack is fixed to 5°. The wing surface is discretized using
∼ 152 000 triangular panels. For that purpose, the Gmsh mesh generator was
used (for more details concerning this software, see Geuzaine and Remacle [32]).
√
The panel sizes are twice finer ( S△ /c ≃ 3.5 10−3) in the wingtip region than
in the wing center. Fig. 2.9 shows the top view and the side view of the wing
discretization. The flow is captured by a resolution fixed to h/c = 3.0 10−3 .
Notice that, in this section the field

U0
sign(ωk ) log 1 + |ωk |
c
is used in color plots instead of the vorticity component ωk itself. It allows to

use the same colorbar limits when the vorticity component varies by several
orders of magnitude.
Results The wing is impulsively started and the simulation was carried out
up to tU0 /D = 3.0. At this time, the total number of grid points is roughly
120 106. 256 processors were used to computed the flow and ∼ 72 hours of
computational time were required. Notice that the origin of the coordinate
system is centered with respect to the wingspan but shifted at the trailing
edge. The coordinate system is such that x, y and z point to the spanwise,
streamwise and normal directions respectively.
If we follow the Prandlt lifting line theory which is in fact only valid for
high aspect-ratio wings (AR ≥ 10) and high Reynolds number flows, the lift
(a) (b)
Figure 2.9: Discretization of the NACA0012 very low aspect ratio wing using ∼
152 000 triangles. The resolution is twice coarser at the center. (a) top view and (b)
side view.
0.06
0.05
0.04
0.03
ΓC
0.02
2 U0 c
1
0.01
1
−0.01
−0.02
−0.03
−0.5 −0.4 −0.3 −0.2 −0.1 0 0.1 0.2 0.3 0.4 0.5
x/b
Figure 2.10: Flow past a very low aspect ratio wing: circulation around each profile
at t Ub0 = 3.
Figure 2.11: Flow past a very low aspect ratio wing: illustration of the contour C
used for the circulation evaluation and of the separation point ys /c.
would then be given by:

Z b/2
Ltot
= U0 ΓC (x)dx
ρ −b/2
where ΓC (x) is the circulation around a contour C enclosing the wing section
at x: Z Z
ΓC (x) = − u · dl = − ωx dydz.
C SC
Fig. 2.10.(b) and Fig. 2.11 show the obtained circulation and the contour chosen
for the calculation respectively. One sees the separation point of the boundary
layer (found at the position ys /c = 0.42). Surprisingly, we observe that the
circulation around the wing for the section at the center is negative. If one
estimates the lift coefficient (reminding the severe assumptions associated to
the Prandlt theory) through the integration of the circulation, one obtains
CL ≃ 0.016. This value has to be compared to the 2D lift coefficient: Cl ≃
5π

2π 180 ≃ 0.55. The performance of the present wing is thus, as expected, very
low (even if it is likely better than CL ≃ 0.016 in reality). The preponderant
factor is indeed the aspect ratio. The wingtip effect appears to be of high
importance as the streamlines have to go around the leading edge extremities.
However, the efficiency is even worse because the Reynolds number is too low;
the boundary layer is laminar and separates sooner.
The slice of the spanwise vorticity component, at the center line, is plotted in
Fig. 2.12.(a). The boundary layer detachment is observed as well as the starting
vortex. Fig. 2.12.(b) shows a 3D perspective view of 6 streamlines and 6 cuts

of the streamwise vorticity component at tU0 /D = 3.0. The global formation
of the trailing vortices is visualized. The counter-rotating vortices, due to
the wing lift, are surrounded by opposite vorticity. This is generated by the
flow circumventing the leading edge extremities. Fig. 2.13 and 2.14 respectively
present various slices of the streamwise vorticity and the streamwise velocity at
various downstream distances. Close to the wing, different recirculation zones
are created. The vorticity field appears to be well captured and is smooth
up to the wall in spite of the immersed boundary. In the wake region, the
complex formation of the counter-rotating vortex pair is observed. The opposite
vorticity remains quite strong. A residual trace of velocity deficit, created by
the boundary layers, is still observed in the vortex-core. The inner-core velocity
deficit is roughly uy /U0 − 1 ≃ −0.1.
(a)
(b)
Figure 2.12: Flow past a very low aspect ratio wing. (a) Slice of the spanwise
vorticity field at the center: sign(ωx ) log(1 + Uc0 |ωx |). (b) 3D perspective view of
6 streamlines and 6 cuts of the streamwise vorticity component, located at x/b =
[−1.0, −0.5, 0.0, 0.5, 1.5, 2.5].
y/c = −1.0
y/c = −0.8
y/c = −0.5
y/c = −0.2
Figure 2.13: Flow past a very low aspect ratio wing. Downstream slices of the
streamwise vorticity field sign(ωy ) log(1 + Uc0 |ωy |).
y/c = 0.0
y/c = 1.0
y/c = 2.0
y/c = 3.0
Figure 2.13: Cont’d.

y/c = −1.0
y/c = −0.8
y/c = −0.5
y/c = −0.2
Figure 2.14: Flow past a very low aspect ratio wing. Downstream slices of the
streamwise velocity field uy /U0 − 1.
y/c = 0.0
y/c = 1.0
y/c = 2.0
y/c = 3.0
Figure 2.14: Cont’d.

2.5 Conclusion
An IB method within the VIC-PFM method has been exhaustively described.
This approach provides satisftactory results for the simulation of the flows past
arbitrary body shape, at moderate Reynolds number, in terms of flow dynamics
and forces acting on the body.
The IB approach conserves the main advantages of the VIC-PFM method[13],
i.e., i) the easy and efficient parallelization of any problem by locally solving
the Poisson equation on each subdomain and ii) the unbounded boundary con-
ditions are implicitly taken into account.
The method has been successfully validated against the test case of the
flow past a sphere at Re = 300. Satisfactory results were obtained in terms of
averaged values, amplitude- and frequency-oscillations.
The method has also been successfully applied to the simulation of the
impulsively started flow past a very low aspect-ratio wing (NACA0012 profile,
AR = 1) at moderate Reynolds number. The complete formation of lift-induced
vortices is captured.
So far, the method uses uniform isotropic grids. The extension to a multi-
resolution framework is described in the next chapter.
Chapter 3
A mesh refinement
technique for the vortex
particle-mesh method
3.1 Introduction
Numerical methods for external wall-bounded flows, that use a uniform resolu-
tion, rapidly lead to an unaffordable number of grid points due to the constraint
of capturing the boundary layers (for which the thickness scales with Re−1/2 ) as
well as the wake. “Classical” vortex particle (VP) methods (i.e. using the Biot-
Savart fast evaluation to compute the velocity field from the particle strength,
as opposed to the hybrid vortex particle-mesh method, that uses a grid-based
Poisson solver) are well-suited to automatically adapt to the local resolution.
However, an appropriate redistribution is performed since a reset of the particle
on a regular lattice has to be done in order to tackle the clustering/depletion
issue. Three main approaches exist to handle varying resolution. i) A global
mapping is used in order to apply the classical redistribution. ii) The mesh
is composed of multiple uniform blocks which have different resolutions. The
redistribution procedure is then adapted to this particular mesh. iii) Wavelet-
particles (see Bergdorf and Koumoutsakos [3]) are used. That allows very fine
local refinement.
38 Chapter 3. A mesh refinement technique for the VPM method
Ould Salihi [64] (see also Cottet and Koumoutsakos [16]) used a mapping of
the redistribution lattice to perform non-uniform simulations. Ploumhans [68]
extended this approach in 3-D to achieve the simulation of the flow past a
sphere for various moderate Reynolds numbers.
Daeninck [21] developed a redistribution technique on hierarchically refined

grids. Classical redistribution scheme are used in order to redistribute a particle
of size hl /h = 1/2, 1 and 2 on a grid h. The length scale, hl , and the strength
of this particle are modified to remain conservative on the grid h. Moreover,
a local and conservative smoothing of the solution is performed. The method
was successfully applied to the flow past a hemisphere at Re = 3000. This
method provides good results for “classical” vortex particle methods. The two
first derivatives, however, oscillate too much when they are evaluated using
finite differences. Another idea is based on average-interpolating wavelets (see
Bergdorf et al. [2, 3] and Rossinelli et al. [73]). This approach was successfully
applied to 2-D and 3-D convection-diffusion equations. In this chapter, a mesh
refinement technique, based on the latter, is presented for the resolution of
the 3-D Navier-Stokes equations using either Direct Numerical Simulations
(DNS) or Large Eddy Simulations (LES). It is also worth mentioning that, in
the context of bluff-body flows, Rossinelli et al. [72] used the GPU-enhanced
performance to significantly reduce the computational cost.
We here use hierarchically refined grids to handle the multiple resolution.

The methodology for solving the Poisson equation remains the same; i.e. the
Poisson equation is solved using a fast grid solver (here using the FISHPACK
library) and the required boundary conditions are provided using the Parallel
Fast Multipole (PFM) method. However, the mesh management (generation &
parallelization) has to be redesigned in order to have a better load balancing.
The current VPM implementation only allows one block per cpu.
The present chapter is structured as follows. A mesh-block generator is first

described. Second, a redistribution scheme with multiple resolution handling
capability is depicted. And, finally, the overall methodology is illustrated and
validated.
3.2. A basic multi-block generator 39
3.2 A basic multi-block generator

The parallel efficiency of a refined mesh technique lies in load balancing. In the
first version of the VPM-PFM code, when using uniform resolution in the entire
computational domain, the mesh was simply split over each processor by equal
subdivisions. This approach leads to important drawbacks and limitations. A
new multi-block mesh generator has thus been developed.
The mesh generator described hereinafter is based on hierarchically refined
grids in the spirit of Adaptive Mesh Refinement (AMR, first developed by
Berger and Oliger [4]). However, we here only present results using a prescribed
mesh size function; it is thus not “adaptive”. The algorithm reads as follow:
1. choice of a mesh size function: the user provides a function δ(x) that
defines the maximal mesh size to be respected over the whole grid.
2. multi-block initialization: a mother-box (level 0) is first created using

the user-provided bounding box.
This mother-box is then equally split in [n0x , n0y , n0z ] in each direction.
The n0x × n0y × n0z boxes are set with a (user-provided) initial resolution
h0 .
3. refinement strategy: if the mesh resolution of a box (at a level l )

is too coarse with respect to the prescribed function δ(x), this box is
successively split into 8 child boxes with a resolution hl+1 = ½hl . This
process is continued until the resolution constraint is satisfied for every
point in the box.
4. coarsening strategy: if the mesh resolution is twice too fine for every
point in a box (typically when the initial mesh size h0 is too constraining),
the mesh size of the box is doubled.
5. smoothing: the ensemble of boxes are smoothed to ensure that the mesh
size of any box is never twice higher/lower than that of its neighbors.
6. block agglomeration: in order to reduce the computational cost as-

sociated with the evaluation of the boundary conditions (required for
solving the Poisson equations and evaluated using the parallel fast multi-
pole PFM method), adjacent equal-resolution boxes are merged such as

to form a bigger parallelipedic block with the same resolution.
7. parallel implementation: the boxes are finally distributed over each

processor using the METIS Library[44, 43]. It optimizes the load bal-
ancing, by both equally sharing the computation and/or memory costs
and by minimizing the communication cost. The former, associated to a
block i, is computed as follows:
i

Ti = C1 C2 NSurf i
log(NPtotart ) + (1 − C2 )NPi art + (1 − C1 )NGrid ,
i
NSurf , NPi art , Ngrid
i
, NPtotart are the number of grid points located on the
surfaces of the block i, the number of particles located in the block i,
the number of grid points located on the block i and the total number
of particles respectively. C1 and C2 lie between [0 1]. C1 determines the
relative importance of the computational cost with respect to the memory
cost. C2 is intrinsic to the numerical code implementation and provides
the relative computational time in the various steps of the algorithm.
The latter, associated to the pair of blocks (i, j), is simply defined to be
proportional to the number of grid points on the common block surfaces:
(i,j)
T(i,j) , T(j,i) = C3 Nsurf ,
C3 fixes the relative importance of the communication cost and depends

on the hardware architecture. The values of C1 , C2 and C3 are empirically
found and also depend on the considered applications.
This new mesh generator library leads to major improvements:
• the number of subdivisions is independent of the number of processors,
• the blocks, when they become useless (no vortex particle in the block),
are disabled “on-the-fly”,
• the blocks can be non-contiguous,
• the parallelization of the multi-block mesh over each processor is versatile

and performed using the METIS [44, 43] library that allows to take into
account both the computational cost and the communication cost.
3.3. A redistribution scheme based on average-interpolation 41
3.3 A redistribution scheme based on average-

interpolation
This section aims at describing the methodology to redistribute the strength

of a particle, which has a size hl , on a uniform grid which has a different
resolution h. The basis of the following approach is an extension of the well-
′
known redistribution schemes (Λ3 , M4 , etc.). Those schemes redistribute on
a uniform grid h a particle p of strength αp and of size hl = h located at the
position xp using a compact kernel K. In 1-D,

X xi − xp
W (xi ) = αp K
p
h
is the fraction of the particle strength αp set to the grid location xi . Depend-
ing on the choice of the kernel, the representation of the particle strength on
the grid h conserves a certain number of moments; e.g., moments up to the
third order are conserved when using the M4′ scheme. The 3-D extension is
straightforward, using tensor products:

X xi − xp yj − yp zk − zp
W (xijk ) = αp K K K .
p
h h h
We here extend this methodology to handle multiple sets of particles which have
various sizes hl /h = 1/2, 1, 2. In particular, patches of particles of different
resolutions can lie over a grid of parameter h. One thus has to i) “detect”
what is the size of the particle that lives above every grid points in order to ii)
redistribute the sets of particles with respect to which particle size is seen by
the grid at those locations.
3.3.1 Multi-level redistribution scheme
One defines (αlp , xp ), the strength and the position of a particle p at a level
l and living on a grid for which the mesh size is h. At a level l, the particle
size is, by definition, 2l h. The value of the level l takes either the value 0 or
±1. Three cases can occur: the particle is redistributed over a grid which has
either the same resolution, or a resolution twice as fine, or a resolution twice
as coarse.
The basic idea is to combine two methodologies: 1) the classical redistribu-
tion procedure and 2) the average-interpolating wavelet (for more information
concerning the wavelets, see Sweldens and Schröder [81]). First generation
average-interpolating wavelets constitute a very efficient mathematical toolbox
if one wants to interpolate uniformly distributed data with an accurate control
over the resolution. Second generation wavelets deal with unevenly distributed
data. It basically constructs the primitive function of the variable of interest
(taken as the averaged value over its area) in order to compute the average on
another area by evaluating this primitive function. The scheme is presented in
1-D; it can be readily extended to three dimensions through the use of tensor
products, just as for uniform redistribution.
1. classical redistribution: Each particle of size hl is first redistributed

onto a grid of the same resolution hl . For that purpose, classical redistribution
′
schemes (such as M4 , Λ3 schemes) are used, see Fig. 3.1.(a). These grid points
are virtual, meaning that they are not created, unless they are at the same
level as the particle.
2. average-interpolating wavelet: Consider a particle not aligned with

the grid. Its strength is redistributed on four points (xi−1 , xi , xi+1 , xi+2 )
l
with four associated strengths (Wi−1 , Wil , Wi+1
l 0
, Wi+2 ). We define Wil as the
redistributed strength at the position xi at the resolution level l,
Z xi +hl /2
Wil = ω(x′ )dx′ .
xi −hl /2
Rx
The primitive function P (x) , −∞
ω(x′ )dx′ is thus known between the redis-
tribution points xi , i.e. at xi + hl /2, since
Z xi +hl /2
P (xi + hl /2) , Pi = ω(x′ )dx′
−∞
i Z
X xi +hl /2
= ω(x′ )dx′
n=0 xi −hl /2
i
X
= Wnl .
n=0
(a)
(b)
(c)
(d)
Figure 3.1: Sketch of the redistribution procedure on refined grids. (a) Redistribution
of the strength of the particle of size hl onto a mesh of size h = hl , (b) Piece-wise
interpolation of the primitive P (x) to provide 2 new strengths on the mesh of size
h = 21 hl , (c) Entire redistribution of the particle strengths of size hl onto a mesh of
size h = 21 hl and (d) redistribution of a particle strength of size hl onto a mesh of
size h = 2hl .
The function P (x) can then be used to compute the strengths at various loca-
tions for different resolution levels, l′ , using
Z xi′ +hl′ /2
Wil′ = ω(x′ )dx′
xi′ −hl′ /2
= P (xi′ + hl′ /2) − P (xi′ − hl′ /2).
• Lifting scheme: We here consider the case with coarse-to-fine average-

interpolation between level l and l − 1. Two weights (Wjl−1 and Wj+1
l−1
)
are created between the points xi and xi+1
Wjl−1 = Pi+½ − Pi
l−1
Wj+1 = Pi+1 − Pi+½

xi +xi+1 1
where, Pi+½ = P 2 = 16 (−Pi−1 + 9Pi + 9Pi+1 − Pi+2 ) is ob-
tained using piecewise third-order polynomial interpolation. This step is
illustrated in Fig. 3.1.(b). The procedure can be done for each interval
[xk−1 , xk ] in order to construct the complete piecewise redistribution as
shown in Fig. 3.1.(c). This lifting procedure provides a smooth redistri-
bution for which the two first moments are conserved. Concerning the
numerical implementation, no intermediate grid hl is created; the parti-
cles of size hl are directly redistributed on the grid hl−1 without using a
grid hl .
• Restriction scheme: the fine-to-coarse average-interpolation simply

consists in concatenating multiple blocks in order to create a coarser
one
Wjl+1 = Pi+1 − Pi−1 .
This step is illustrated in Fig. 3.1.(d). Only the zero order moment is
conserved.
3.3.2 Treatment of varying resolutions

We have just introduced a scheme which is able to redistribute a particle of
size hl on a grid of size hl , hl /2 or 2hl . The direct application of this scheme
to a set of particles with different resolutions would lead to oscillations. In-
(a)
(b)
Figure 3.2: Sketch of the ghost particle (empty circles) creation from refined grids
using “real” particles (filled circles). (a) Creation of coarse ghost particles (b) creation
of fine ghost particles.
deed, the smooth interpolation of particle quantities requires a continuous and

overlapping particle arrangement.
To tackle this problem, Daeninck [21] applied a global smoothing step which
conserves the first three moments of the solution. This approach unfortunately
proved unadapted to the use of grid-based methods (as it is the case here). We
indeed need to preserve the smoothness of the solution. The present approach
will rather rely on ghost particles. This buffer region maintains the overlap
of particles up to the resolution interface and preserves the smoothness of the
solution. We illustrate the procedure on a 1-D example (see Fig.3.3.(a) for the
h
illustration). The physical or “real” particles of size 2, h and 2h lie over a grid
of size h.
First, we generate our ghost particles in a buffer region to extend the “real”
ones (as described in Fig. 3.2). Knowing the complete solution on the “real”
domain, we initialize the “real” and the “ghost” particles. Notice that both
“real” and “ghost” particles are initialized using the same characteristic length
as that of the generating grid. The procedure in order to create those ghost
particles depends on the resolution of the neighbor grid. The procedure can be
summarized as follows:
1. When the neighbor grid uses the same resolution, “ghost” particles are
created as the “real” particles; i.e. by identifying the particle strength
with the grid strength located at the same position.
2. When the neighbor grid is twice finer (see Fig.3.2.(a)), neighbor “real”
particles (e.g. 2 in 1-D and 8 in 3-D) are concatenated into a “ghost”
particle.
3. When the neighbor grid is twice coarser (see Fig.3.2.(b)), neighbor “real”
particles together with previously created neighbor “ghost” particles are
used for the interpolation to create fine “ghost” particles. It is worth
noting that this step should be necessarily done at the end as it uses
“ghost” particles created during the step 2.
Each of these ghost particles therefore carries either the same vorticity as
several finer real particles, or a fraction of the vorticity carried by a coarser
particle. This can be visualized in Fig.3.3.(a) where we represent the handling
h
of the interface between resolutions h and 2: ghost particles of size h are
h
created above the real particles of size 2. In the same way, ghost particles of
h
size 2 are created above the real particles of size h.
Notice that the coarse-to-fine redistribution (i.e. redistribution of particles
with size 2h onto grid h) require a larger support. The number of created
ghost particles thus also depends on the support of the redistribution scheme.
In the previous section, we have seen that a consistent M4′ interpolation will re-
quire two neighbor particles at the same level, assuming a moderately distorted
particle set. When carrying out a coarse-to-fine interpolation, this number of
required neighbors is increased by the support of the average-interpolation as
the support of M4′ is broader than that of the third order cubic spline. For the
fine-to-coarse interpolation, it is however not the case.
In a second step, the size of the real particles are redistributed onto the
grid; one basically tags each grid points with the size of the real particles living
above. The third step then consists in a selective particle-to-mesh interpolation.
A mesh node with a tag level l can only receive the redistributed information
3.4. A 1-D multiple resolution example 47
of particles at a level l. This step plays the same role as the indicator function
introduced by Bergdorf.
The extension in two- and three-dimension is readily done by a straightfor-
ward application of the procedure described here-above. Fig.3.3.(b) illustrates
the easy 2D extension of the procedure. One determines region of the grid in
order to only use a unique particle size for the redistribution at a given grid
point.
Real and ghost particles are convected, ensuring a proper overlap of particles
with the same size up to the interface between two resolutions. Because the sets
of real and ghost particles are a self-consistent representation of the vorticiy
field, no limitation on CFL is required, e.g. a particle can venture arbitrarily far
into a patch of a different resolution. This is actually still a condition for this
self-consistency. The thickness of the ghost particle buffer must be sufficient
for the interpolation wavelet procedure. Those ghost particles can indeed get
deformed by the flow, requiring instead a Lagrangian CFL condition, just as
for the real particles.
3.4 A 1-D multiple resolution example
As an illustration, let us consider a 1-D Gaussian function f (x) with σ 2 /h2 = 30

and discretized using three sets of particles over the interval x/h = [−30, 30].
Particles of size 2h, h, h/2 are discretizing the sub-intervals [−30, −10], [−10, 10]
and [10, 30], respectively. These particles are then redistributed on a uniform
grid h using the methodology described before. Fig. 3.4 shows the three sets
of particles (and their respective strength αp = f (xp )h), the Gaussian function
and the discrete function f˜(xi ) = αi /h obtained after redistribution: the func-
tions f and f˜ are superposed. The absolute errors made on the function f and
on its two first derivatives,

Ef (x) = f − f˜,
df
df˜
Ef ′ (x) = − ,
dx dx
d2 f
d2 f˜
Ef ′′ (x) = 2 − 2 ,
dx dx
(a)
* * * * *
* * * * * * * * * * *
* * *
(b)
Figure 3.3: Sketch of the grid point tagging and the use of the ghost particles: (a) in
1D and (b) in 2D. Filled in circles “real” particles and circles “overlapping” particles.
are plotted in Fig. 3.5. The errors made on f and its first derivative (see
Fig. 3.5.(a)) follow the same evolution and are seen to be of the same order
of magnitude. However, the errors made by the redistribution of a particle on
a grid of a different size are at least 10 times higher than those made by the
redistribution on a grid of the same size. Concerning the second derivative (see
Fig. 3.5.(b)), the error made for each redistribution case is roughly continuous
and remains relatively low, but yet ten times higher than those made on the
function and its first derivative. Concerning the interpolation of coarse particles
on a finer grid, the errors made are indeed higher, yet are better than what
we could expect as we started with twice less information. Regarding the
3.5. Vortex ring test case 49
1.5
[−]
0.5
0
−30 −20 −10 0 10 20 30
x/h
Figure 3.4: 1-D redistribution of a Gaussian function. Set of particle with various
size and their strength αp : the filled circles correspond to the “real” particles
2 and the
empty circles to the ghost particles. The Gaussian function f (x) = exp − xV , (solid
line) and the redistributed function f˜(x) (superposed dashed line).
interpolation of fine particles on a coarser grid, we can definitely conclude that

the added average interpolation error does not constitute an issue.
The errors, made on the function f and its two first derivatives, are low
enough. This is essential as, in our case, they correspond to the discretization
of velocity, the vortex-stretching term and the diffusion term.
3.5 Vortex ring test case

In order to assess the quality of the refined mesh technique, the simulation of
a vortex ring has been carried out. The vortex ring is first discretized using a
uniform resolution, then it crosses various resolution regions.
Initial condition The initial azimuthal component of the vorticity field is

taken as Gaussian, 2
Γ0 s
ωθ = 2 exp − 2 ,
πσ0 σ0
(a) (b)
−3
−4
x 10 x 10
2.5
2
Ef (x), Ef ′ (x)h
Ef ′′ (x)h2
2
1.5
1
1
0.5
0 0
−30 −20 −10 0 10 20 30 −30 −20 −10 0 10 20 30
x/h x/h
Figure 3.5: 1-D redistribution of a Gaussian function. (a) error on the function f
(solid), on its first derivative (dash) and (b) on its second derivative.
where s2 = z 2 +(r −R)2 , Γ0 is the initial circulation, σ0 the initial core size and
R the ring radius. One defines the time-scale as t0 = R2 /Γ0 . The Reynolds
number based on the circulation is Re = Γ0 /ν = 5500. The initial core size
is fixed to σ0 /R = 0.4131 (same case as Shariff et al. [75, 76] and Cocle [11]).
The time step used for the simulation is ∆t/t0 = 0.05. Three simulations
have been done: i) one reference simulation using uniform resolution (fixed to
h/R = 0.04) in the entire domain, ii) a simulation where the vortex ring starts
from a resolution h then goes through a region of resolution 2h and then finally
back to a region with the initial resolution h, and iii) a simulation where the
vortex ring starts from a resolution h then goes through a region of resolution
h/2 and then finally back to a region with the initial resolution h. Fig. 3.6
shows a restricted view of a cut through the mesh for each investigated case.
Results Fig. 3.7 presents a cut of the vorticity field for each case at various
time: t/t0 = [0.0, 4.0, 6.0, 9.0, 13.0, 20.0]. The vortical structures smoothly
cross the different resolution areas, from a finer to a coarser (FTC) resolution as
well as from a coarser to a finer resolution (CTF). No reflection or oscillations
are observed at the transitions. Notice that the FTC transition is the easiest
case because the redistribution acts like a smoother. The structures still remain
accurate, even if the vortex-core is captured using only σ0 /(2h) ≃ 5 grid points.
This can be especially observed at time t/t0 = 9.0; at that time the vortex-core
is clearly under-resolved. Concerning the CTF transition, the lifting scheme
works properly: the interpolation does not create spurious small scales. At time
t/t0 = 20.0, the vortices reach the initial resolution region, their core is finally
convected at the same position for each of the three cases. That indicates that
the linear impulse is qualitatively well conserved (as seen below). No major
difference is seen in terms of the physics; even the thin wake is still captured.
Fig. 3.8, 3.9 and 3.10 respectively show the time-history of the three compo-
R R
nents of the total vorticity Ω , V ωdV , of the linear impulse I , V x × ωdV
R
and of the angular impulse A , 13 V x × (x × ω)dV . Since a free vortical ring
is considered, those quantities have to be conserved.
However, because the grid-to-particle interpolation procedure is not con-
′
servative (as we use the M4 scheme to interpolate information from the mesh
to deformed scattered particles) the first three moments are then not exactly
conserved. This is particularly due to the particle-set distortion (the so-called
“clustering/depletion” effect); the particles do not remain on a regular lattice.
Those moments are nevertheless qualitatively well conserved. More precisely,
concerning the time-evolution of the total vorticity (Ω) and of the linear im-
pulse (I), the global trend is quite similar to each case. Concerning the angular
impulse (A), the FTC case provides better conservation. It is worth noting
that the present code works in simple-precision: those quantities are at least
10−8 . Moreover, it is worth noting that no re-projection of ω into a divergence
free field was done in the case with uniform grid, because the divergence mea-
sure did not reach the user-prescribed divergence threshold. That also explains
why the total vorticity is not mathematically conserved to zero.
52
Chapter 3. A mesh refinement technique for the VPM method
Figure 3.6: Vortex ring at Re = 5500. Slice of the mesh at y/R = 0: case using uniform resolution (center), case using locally a coarser
resolution (left) and case using locally a finer resolution (right).
t/t0 = 0.0
3.5. Vortex ring test case

Figure 3.7: Vortex ring at Re = 5500. Slice of ωy t0 at y = 0 for the simulation using a uniform resolution (center), using locally a coarser
resolution (left), and using locally a finer resolution (right).
53
Figure 3.7: cont’d

t/t0 = 4.0

t/t0 = 6.0

t/t0 = 9.0
t/t0 = 13.0
3.5. Vortex ring test case

57

t/t0 = 20.0
−5
x 10
(a) 3
0
Γ0 R
1
−1
Ωx
−2
−3
−4
−5
0 2 4 6 8 10 12 14 16 18 20
t/t0
−5
x 10
(b)
0
−2
−4
Γ0 R
1
Ωy
−6
−8
−10
0 2 4 6 8 10 12 14 16 18 20
t/t0
−5
x 10
(c) 2
−1
Γ0 R
1
−2
Ωz
−3
−4
−5
−6
0 2 4 6 8 10 12 14 16 18 20
t/t0
Figure 3.8: Time-evolution of the three components of the total vorticity: i) using a
uniform resolution (solid line) ii) using locally a coarser resolution (dashed line) and
iii) using locally a finer resolution (dash-dotted line).
−4
x 10
(a) 2
1
Γ0 R2
1
0
Ix
−1
−2
0 2 4 6 8 10 12 14 16 18 20
t/t0
−4
x 10
(b)
1
Γ0 R2
1
Iy
−1
0 2 4 6 8 10 12 14 16 18 20
t/t0
(c) −3.39
−3.395
−3.4
Γ0 R2
1
Iz
−3.405
−3.41
−3.415
0 2 4 6 8 10 12 14 16 18 20
t/t0
Figure 3.9: Time-evolution of the three components of the linear impulse: i) using a
uniform resolution (solid line) ii) using locally a coarser resolution (dashed line) and
iii) using locally a finer resolution (dash-dotted line).
−4
x 10
(a)
2
Γ0 R3
1
Ax
0
0 2 4 6 8 10 12 14 16 18 20
t/t0
−4
x 10
(b)
2
Γ0 R3
1
1
Ay
−1
0 2 4 6 8 10 12 14 16 18 20
t/t0
−6
x 10
(c)
2
Γ0 R3
1
1
Az
−1
−2
0 2 4 6 8 10 12 14 16 18 20
t/t0
Figure 3.10: Time-evolution of the three components of the angular impulse: i)

using a uniform resolution (solid line) ii) using locally a coarser resolution (dashed
line) and iii) using locally a finer resolution (dash-dotted line).
3.6 Multiple resolution simulation of the flow

past the sphere at Re = 300
U∞ D
The DNS of the flow past a sphere at Re = ν = 300 has also been carried out
using the refined mesh technique described before. The immersed boundary
technique is also used in order to perform this simulation. This is thus the
VPM-PFM-IB code with mesh refinement added. The present validation case
is basically the same as that using the immersed boundary technique (described
in Chap. 2, Sect. 3) with uniform resolution. The latter simulation is definitely
too resolved in the wake region, and led to too many grid points (∼ 120 106).
Setup The multi-resolution simulation uses the same resolution in the near-
wall region h0 /D = 1/75. The wake is captured up to 38 diameters down-
stream, thus comparable to the case with uniform resolution. The time-step is
taken as before ∆tU∞ /D = 8.5 10−3. The maximal CFL criterion were seen

to be CF L = |u| h ∆t ≤ 0.89. The resolution function δ(x) is taken as
max
follows: 
 h0 r


 D if D ≤ 0.8


h0
r
r
δ(x) = D 25 D − 19 if 0.8 < D ≤ 1.0

 r −1


 0 D
6 hD 10
 19

3 otherwise
where r = |x − xc | is the distance of a point x to the sphere center xc . Vortex

particles, that go out of the mesh, are simply removed. In opposition to the
case using uniform resolution (see Chap.2, Sect. 3, or in [21, 11]), no specific
outflow condition is used. In fact, the varying resolution method can itself
be used as a way to properly do wake flow with outflow boundary condition:
simply running a very long domain with resolution getting coarser and coarser
provides a natural way of having a good outflow condition.
Results Fig. 3.11.(a), 3.11.(b) and 3.11.(c) show a slice of the vorticity field
for each component. Typical results are obtained. They are also compara-
ble to those obtained by Ploumhans [68]. A cut of the mesh is presented in
Fig. 3.11.(d). The blank squares located in the up and down corners correspond
to disabled block because there are particle free. The ratio between the coarsest
3.6. Multiple resolution simulation of the flow past a sphere 63
and the finest mesh size is 32. The total number of grid points needed for this
simulation is ∼ 11 106, to be compared to the case using uniform resolution
(i.e. ∼ 120 106). The gain is thus substantial. A zoom cut of the vorticity field
close to the surface of the sphere together with the mesh is shown in Fig. 3.12.
It is seen that the vorticity field is smooth up to the surface, indicating that
the simulation is well-resolved. Fig. 3.13 presents 3 views of a λ2 surface. The
typical hairpin vortices shed behind the sphere are visualized; four of them
are observed. The time-history of the lift and drag coefficients are shown in
Fig. 3.14. Table. 3.14 summarizes the global values of interest. Those values
are computed as for the case using uniform resolution; i.e. using the Noca’s
formula over a control volume. The control volume is also taken as done for
the case with uniform resolution (i.e. a cube centered in the sphere center and
with a edge length fixed to 1.6D). Those are compared with our simulation
using uniform resolution and with those of literature. The values are seen to be
in fair agreement. The averaged lift is a bit too high; but one must recall that
it is quite difficult to evaluate. The values could also be affected by a maybe
too coarse level of the resolution after 17.5D
(a)
64
(b)
(c)
(d)
Figure 3.11: Multi-resolution

simulation of
the flow past
a sphere at Re = 300:
global view of a slice of the three vorticity components: (a)
D D D
sign(ωz ) log 1 + U0 |ωz | , (b) sign(ωy ) log 1 + U0 |ωy | and (c) sign(ωx ) log 1 + U0 |ωx | . (d) cut of the refined mesh.
3.6. Multiple resolution simulation of the flow past a sphere

D
Figure 3.12: Multi-resolution simulation of the flow past a sphere at Re = 300: zoom of the vorticity component sign(ωz ) log 1 + |ωz | .
65
U0
A cut of the mesh used for the simulation is also shown.
(a)
66
(b)
(c)
Figure 3.13: Multi-resolution simulation of the flow past a sphere at Re = 300: surface of the λ2 criterion: (a) top view, (b) side view and
(c) offset view.
3.7. Conclusion 67
0.7
Force Coefficient [-] 0.6
0.5
0.4
0.3
0.2
0.1
0
150 155 160 165 170 175 180 185 190 195 200
tU0 /D
Figure 3.14: Time-evolution of the drag coefficient (solid line) and lift coefficient
(dashed line).
Table 3.1: Flow past a sphere at Reynolds 300: summary table of global quantities
compared to the literature
CD ∆CD CL ∆CL Str

−3 −2
Multi-resolution 0.677 3.3 10 0.070 1.4 10 0.134
Uniform resolution 0.672 3.4 10−3 0.065 1.5 10−2 0.134
Johnson & Patel [42] 0.656 3.5 10−3 0.069 1.6 10−2 0.137
Tomboulides & Orszag [83] 0.671 2.8 10−3 - - 0.137
Constantinescu et al. [14] 0.655 - 0.065 - 0.136
−3 −2
Ploumhans et al. [68] 0.683 2.5 10 0.061 1.4 10 0.135
Georges [31] 0.661 2.5 10−3 0.066 1.3 10−2 0.134
3.7 Conclusion
We have developed a mesh refinement technique for particle-mesh methods.
It combines the approach of Bergdorf et al. [2, 3] and Rossinelli et al. [73] to
achieve, at the same time, block-wise grid refinement, high order accuracy, and
CFL relaxation.
We have extended this approach for the simulation of 3D incompressible
flows, possibly with solid bodies. A parallel Poisson solver has thus been in-
tegrated. We here successfully applied the approach of Cocle et al. [11] in the
mesh refinement context. That is, the required Poisson equation is here solved
using an efficient combination of a fast grid-based Poisson solver with a parallel
fast multipole method. The method has been validated for 3D incompressible
flows with solid bodies. A simulation of a vortex ring crossing different resolu-
tion regions has been done. It was seen to qualitatively conserve well the three
first moments for each transition.
In order to handle refined meshes, a mesh generator based on hierarchically
refined grids was also developed. It enables to have a good load balancing
taking into account both the computational cost and the memory cost, using
the METIS library.
The multi-resolution simulation of the flow past a sphere at Reynolds 300
has also been carried out. The results have been compared to the equivalent
simulation using uniform resolution and to those of literature; a good adequacy
has been obtained in terms of the time-evolution of the forces acting on the
body.
Chapter 4
Investigation of velocity
deficit effects on aircraft
wake rollup
This chapter is the follow-up of a technical report of the FAR-Wake project

entitled “Rollup of a temporally-evolving wing wake with velocity deficit“ [19].
European project AST4-CT-2005-012238, Report TR2.2.2-3, January, 2007.
The authors are C. Cottin, G. Winckelmans, T. Lonfils, R. Cocle.
4.1 Introduction
Vortex wake (WV) dynamics have been extensively investigated for decades. It
is now well-established that the wake behind a lifting-body essentially consists,
after rollup, in a pair of counter-rotating vortices. When created by an aircraft,
those vortices can be potentially hazardous for a following aircraft. The WV
intensity scales with the lift and thus with the weight. The WV persistence
is at the origin of the separation standards between two aircraft for landing
and takeoff phases. This limits the capacity of the busiest airports. Much re-
search focuses on the way to diminish their intensity by means of active/passive
control systems. Another avenue of research concerns the reduction of those
standards with respect to the meteorological data, measured and/or forecast
70 Chapter 4. Investigation of velocity deficit effects on aircraft wake rollup
(wind, turbulence, temperature). Indeed, the standard separation from ICAO

(International Civil Aircraft Organization) are often too conservative. Exper-
imental and numerical investigations are thus required to further understand
the mechanism of formation and of destruction of wake vortices. See the gen-
eral review of Spalart [77], for more information concerning aircraft trailing
vortices, including the formation, motion and persistence of wake vortices and
relevance to air travel applications.
This work is dedicated to the investigation of the effect of a velocity deficit,
due to the friction drag, on the near-wake rollup dynamics on the net-result
after rollup. It is the follow-up of the research carried out by Cottin et al. [19]
in the framework of the EC funded European FAR-Wake project. Three-
dimensional simulations were then performed. Although they were interesting,
they lacked relevance since only time-developing (T-D) simulations were con-
sidered. A T-D simulation evolves in-time an initial condition using periodic
condition in the axial direction and thus assumes a streamwise homogeneous
direction. The reality is much different; the flow evolves in-time as well as
spatially. In this chapter, space-developing simulations (S-D) were performed,
using the refined mesh handling capability as described in Chapter 3, in order
to obtain more realistic results. A comparison to the T-D case is also presented
for completeness.
A characterization of the net result of wake rollup, in terms of tangential
velocity profile, axial velocity, and their evolution in the wake, is provided,
using the results of the S-D simulation. The fits of these profiles with the most
usual vortex models is also done. In particular, we show that the two length-
scale models of Proctor-Winckelmans and of Jacquin provide a very good fit
(both in terms of maximum tangential velocity and effective core size).
4.2 Initial conditions

This chapter aims at investigating the effect of the axial velocity deficit due
to boundary layers (friction drag) on the roll-up phase and on the resulting
vortex system. Two S-D simulations have been carried out with and without
velocity deficit. The former uses a modeled external-force corresponding to the
wing presence, and the latter uses an inflow modeled vorticity field. A further
4.2. Initial conditions 71
T-D simulation with velocity deficit has also been performed using the modeled
vorticity field as an initial condition.
The first approach consists in directly adding an extra force-term in the

RHS of the Navier-Stokes equations, corresponding to the force exerted by the
wing on the fluid. That entails to model the force distribution. The second
approach consists in modeling the vorticity field shed by a regularized lifting
line at the trailing edge of a generating wing in order to use it as an inflow
condition.
Concerning the SD simulation with velocity deficit, we noticed that non-

physical results were obtained when using a modeled inflow. Hence, we modified
our approach and developed, instead, a modeled external-force. Indeed, when
using the modeled vorticity field with velocity as an inflow, the vortex-pair
circulation is found to significantly increase with respect to the downstream
distance, by ∼ 20%. This is due to a net spanwise vorticity flux at the inflow.
The imposed wake vorticity is up-and-down symmetric but the streamwise ve-
locity is not (because of the lifting-line). As a consequence, there is a spurious
vorticity into the domain. The injected vorticity is not divergence-free, the suc-
cessive reprojection redirects the spanwise vorticity to the streamwise direction.
Those results (diagnostics and visualizations) are provided in Appendix A. It is
worth noting that the problem does not occur for the S-D simulation without
velocity deficit as the spanwise vorticity is zero; indeed there is only streamwise
vorticity resulting from the regularization of the lifting line.
The two approaches use the Prandtl’s lifting line theory. This theory con-
sists in replacing the lifting wing of lift distribution, l(y), by a vortex line
defined by its span loading, Γ(y), using the relation: l(y) = ρ U0 Γ(y) (ρ being
the fluid density and U0 the flight speed). For an elliptically loaded wing, one
has: s 2
y
Γ(y) = Γ0 1−
b/2
where Γ0 is the half plane total circulation. The vortex spacing after roll-up,
b0 is also obtained from
Z b/2
Γ(y)dy = Γ0 b0 ,
−b/2
π
with b0 = 4 b for elliptical loading. Considering the definition of the total wing
lift, and the lift coefficient, respectively:

Z b/2 Z b/2
L = l(y) dy = ρ U0 Γ(y) dy = ρ U0 Γ0 b0
−b/2 −b/2
L 2 Γ0 b 0
CL = 1 2
=
2 ρ U0 S
U0 S
b2
with S = AR the wing surface (AR being the wing aspect ratio), the charac-
teristics of the wake can be expressed as a function of the wing characteristics
as:
1 CL
Γ0 b 0 = U0 b2
2 AR
One also obtains, for the downwash velocity of the vortex wake of after rollup:
2
Γ0 U0 b CL
V0 = =
2 π b0 4π b0 AR
and thus, for an elliptical loading,
V0 4 CL
= 3
U0 π AR
CL
In the present work, the wing characteristics are fixed such that AR = 0.20 and
V0 −2
thus, U0 = 2.58 10 . The Reynolds number is taken as very high (i.e. nearly
realistic):
Γ0
ReΓ = = 106 .
ν
Hence, subgrid-scale modeling will be used.
Description of thervorticity field inflow model For a wing of the cir-

2
y
culation Γ(y) = Γ0 1− b/2 , the corresponding shed vorticity field is a
vortex sheet with circulation per unit length :
dΓ
γ(y) = − (y)
dy
This singular vortex sheet is here regularized using a Gaussian kernel in order
to produce a regular vorticity field:
Z b/2
1 (y − y ′ )2 + z 2
ωσ (y, z) = exp − γ(y ′ ) dy ′
πσ 2 −b/2 σ 2
where σ is a regularization parameter, set in the present work to σ/b = 1/75,

hence still giving a fairly thin vortex sheet.
The resulting vorticity field is presented in Fig. 4.1. It corresponds to the

vorticity field imposed at the inflow for the S-D simulation without velocity
deficit. The time-developing simulation (with velocity deficit) uses this axial
vorticity field and adds a spanwise vorticity field that corresponds to an added
wake component.
Since the axial velocity deficit wake is a model for the net result of the
boundary layers developed along the chord of the wing surface, the same elliptic
distribution (as for the chord and lift) is here used :
s 2
y
uw (y, 0) = uw0 1−
b/2
where uw = ux −U0 is the velocity deficit, with uy the axial velocity component
and U0 the freestream velocity. The parameter uw0 determines the magnitude
of the velocity deficit. The axial direction, x, is defined positively when going
downstream of the wing: the velocity deficit, uw , is thus negative. It is also
further regularized to obtain a regular axial velocity field uw,σ (y, z). The ve-
locity deficit is thus defined by two parameters: the magnitude, uw0 , and the
thickness (set by the regularization parameter σ).
The approach used to find realistic simulation parameters consists in con-

sidering the different contributions to the total drag. As mentioned in [24], the
induced drag contributes to 40 to 50% of the total drag in cruise. Since the
induced drag is well approximated by the cross-flow kinetic energy in a close
plane behind the wing [25], it is obtained by computing the cross-flow kinetic
energy per unit length of the initial vorticity field described above:
Z Z
1 1
Di ≃ Ec0 = ρ (u2y + u2z ) dS = ρ ψx ωx dS
2 S 2 S
Figure 4.1: Initial axial vorticity field, ωx t0 , for the reference case (wake without
velocity deficit).
where uy and uz are the velocity components in y− and z−direction, S is the

cross plane (orthogonal to the flight direction x), ωx is the x−component of
vorticity and ψx is the stream-function satisfying the Poisson equation, ∇2 ψx =
∂ 2 ψx ∂ 2 ψx
∂y 2 + ∂z 2 = −ωx .
The total drag is the sum of the induced drag and of the “profile” (viscous)
drag. In the configuration investigated here, the velocity deficit is a direct
consequence of the development of the boundary layers on the airfoil. The
velocity deficit must thus be defined in a such way that the corresponding
“profile” drag contributes to 50 to 60% of the total drag. For the spatially-
evolving wake, the “profile” drag is obtained as:
Z
Dp = ρ (U0 − ux ) ux dS
S
with U0 the freestream velocity and ux axial velocity. In the present case,
when a time-developing simulation is considered, the relative velocity deficit is
as uw = −(U0 − ux ), and thus, the profile drag is obtained as :
Z
Dp = −ρ (U0 + uw ) uw dS
S
In order to determine uw0 , the fact that the initial flow field corresponds to
that shed at the trailing edge of the generating wing is considered. Concerning
the axial velocity field, the no-slip condition imposes the velocity to decrease
to zero at the trailing edge (centerline of the initial vortex sheet). Thus, from
this point of view, the natural value for uw0 is −U0 , corresponding to u = 0.
Concerning the regularization parameter: the same value of σ is used in
both models. Thus, the boundary layers have roughly the same thickness as the
vortex sheet shed by the wing (from lifting line theory). This set of parameters
leads to an induced drag, Di , and a profile drag, Dp , contributing respectively
to 42.4% and 57.6% of the total drag, which is realistic for cruise conditions.
The initial velocity profiles, and the corresponding vorticity and velocity
fields obtained are presented in Figs. 4.2 and 4.3, respectively. The initial
condition presented in figure 4.3 corresponds to the velocity deficit configura-
tion. The combination of the initial conditions presented in figures 4.1 and 4.3
corresponds to the wake with velocity deficit.
0.04
0
0.03
0.02 −0.2
ux /U0 − 1
0.01
−0.4
z/b
0
−0.6
−0.01
−0.02 −0.8
−0.03
−1
−0.04
−1 −0.8 −0.6 −0.4 −0.2 0 0.2 −0.4 −0.2 0 0.2 0.4 0.6
ux /U0 − 1 y/b
Figure 4.2: Initial axial velocity profiles: (a) (ux − U0 )/U0 at y/b = 0 (the symbols
represent the numerical field as discretized with the grid spacing h); (b) ux /U0 − 1 at
z/b = 0.
Description of the force distribution model As explained before, the

lift per unit length, l(y), is related to the circulation through the relation
l(y) = ρU0 Γ(y). For an elliptically loaded wing,
s 2
y
l(y) = ρU0 Γ0 1− ,
b/2
s
2
2 2 CL y
= ρ U0 b 1− .
π AR b/2
The lift force is further smoothed using a Gaussian kernel. The lift force per
unit of volume is taken as:
2
1 x z2
fL (x, y, z) = l(y) exp − 2 − 2 ,
πσx σz σx σz
(a)
(b)
(c)
(d)
Figure 4.3: Initial velocity and vorticity fields for the velocity deficit configuration:
(a) Axial velocity, ux /U0 − 1 ; (b) Spanwise vorticity, ωy t0 ; (c) Vertical vorticity,
ωz t0 ; (d) transverse vorticity norm, ωyz t0 = (ωy2 + ωz2 )1/2 t0 .
where σx and σz are regularization lengths. We here set σx /b = 1/15 and

σz = 1/75 (same value of the vertical regularization parameter as that used
for the inflow model). The regularized induced-drag force per unit volume is
computed as a “fraction” of the lift using the downwash ǫ(y) = ǫ; here taken
as uniform since an elliptically loaded wing is considered which, according to
the Prand’l theory, gives a uniform downwash:
fi (x, y, z) = ǫ(y) fL (x, y, z)

CL
= fL (x, y, z)
πAR
It further makes sense to take into account the velocity deficit effect by adding
a model for the friction drag force ff (x, y, z) proportional to the chord, and
thus here also elliptically loaded. The friction drag force per unit volume is
thus taken of the same form as that for the induced drag:
ff (x, y, z) = ξfi (x, y, z)
ξ is here chosen to be unity, consistent with the fact that, for cruise condition,
the induced drag contributes to 40 · · · 50% of the total drag (see de Bruin [24]).
That also leads to a “lift-to-drag ratio” (or “glide-ratio”) equal to:
L 1
= ≃ 7.9
Dtot (1 + ξ) πACL
R
We also connect the 2-D rollup simulations and the 3-D S-D simulations
using:
x = U0 t = U0 t0 τ,
with τ = t/t0 the dimensionless time for the 2-D simulation, with t0 = b0 /V0
(where V0 = Γ0 /(2πb0 ) is the wake descent velocity). A station at x/b in the
3-D S-D simulation thus has an equivalent dimensionless time equal to
x/b
τ= 2
= 3.29 10−2 x/b.
U0 2πb0
b Γ0
4.3 Numerical setup
The numerical setup of the three simulations is herein described. We first

present the two S-D simulations, performed with and without modeling of a
velocity deficit (due to the boundary layers). Second, the parameters used for
the time-developing simulation are depicted.
S-D simulation without velocity deficit The simulation uses an inflow

which is the corresponding axial vorticity field (i.e. vortex sheet) shed by an
elliptically loaded wing without the influence of the friction drag. This vortex
sheet is regularized using the value σ/h = 1/75. An outflow condition is also
used for which one imposes a “no-through flow” at the outgoing plane, located
x
at b = 10; the resulting constraints on the vorticity field read (for a y − z plane
located at x = xo ):
ωn (xo − x, y, z) = ωn (xo + x, y, z)
ωt (xo − x, y, z) = −ωt (xo + x, y, z)
where the subscripts n and t stand for normal and tangential. This leads to
un (xo , y, z) = U0 . A uniform resolution is used and fixed to h/b = 1/160. The
time-step is chosen as ∆t/t0 = 6.45 10−5. The total number of grid point is
roughly ∼ 45 106.
S-D simulation with velocity deficit An external force effect is added

to take into account the force exerted by the wing on the fluid; i.e. the lift
and the total drag (friction and induced). In practice, as the vorticity-velocity
formulation of the Navier-Stokes equations is used, ones adds the curl of this
external force, computed on the mesh.
An adaptive time step is used (see Koumoutsakos [46]). The time constraint

is based on the strain tensor S , 12 ∇u + (∇u)T :
∆t|S|max < 0.1,
where |S|max is the maximum of the sum of the absolute values of the strain
4.3. Numerical setup 79
tensor lines:  
X
|S|max , max  |Sij |
i
j
For this LES, we here use a “complete-small” subgrid-scale model (referred as

RVM model in Cocle et al. [12] ). The turbulent viscosity scales with the com-
plete strain tensor field while it only acts on small scales of the LES field. For
this simulation, one benefits from the multi-resolution capability as described
in Chapter 3. The mesh size function δ(x) is defined as follows

 h0 x


 b if b ≤ 2.0



 2 h0

if 2 ≤ x
≤6
b b
δ(x) =
0
4 hb x



 if 6 ≤ b ≤ 10



 10(x − 10) + 4 h0

if x
≥ 10
b b
where h0 /b = 1/250. Notice that the multi-resolution methodology permits

to use a much better resolution in the near-field; here defined as being the
region up to two wingspans. The entire flow is captured up to 16 wingspans.
From x/b = 10, the resolution is forced to deteriorate, constituting an outflow
condition. No further outflow condition is used: the vortical structures are
simply coarsened, convected and, finally, removed. One defines the total energy,
E, and the enstrophy, ε, of the system
Z
1
E = ω · ΨdV,
2 V
Z
1
ε = ω · ωdV.
2 V
The time-history of the energy and of the enstrophy are presented in Fig. 4.4.(a)
and 4.4.(b) respectively. At the beginning, the energy and the enstrophy
rise. They eventually reach a statistically stationary regime from the time
tU0 /b ≃ 17. Yet the simulation was run up to tU0 /b ≃ 47. The average and
the fluctuation average of the flow are post-processed using 91 instantaneous
fields in the time frame [17 − 47]. The divergence level is maintained at a
low level . 4 10−4 (see Fig. 4.4.(c)). The history of the time step is shown in
Fig. 4.4.(d). For information, the time step, scaled with the maximum vortic-
ity, is also presented and was found to be roughly ∆t max(|ω|) ≃ 0.20 · · · 0.35.
The simulation was carried out using 128 cpu’s during roughly 2 weeks. The
total number of grid point was ∼ 90 106.
(a) (b)
0.18 3.5
0.16
3
0.14
E/(U02 b3 )
0.12 2.5
ε/(U02 b)
0.1
2
0.08
0.06 1.5
0.04
1
0.02
0 0.5
0 5 10 15 20 25 30 35 40 45 50 0 5 10 15 20 25 30 35 40 45 50
tU0 /b tU0 /b
(c) (d)
−4
x 10
10 0.4
9 0.35
|∇ · ω|dV
8
0.3
7
0.25
[−]
6
0.2
V
5
R
b2 1
0.15
U0 V
3 0.1
2 0.05
0 5 10 15 20 25 30 35 40 45 50 0 5 10 15 20 25 30 35 40 45 50
tU0 /b tU0 /b
Figure 4.4: S-D simulation of the wake rollup with velocity deficit: (a) time-history
of the energy, (b) of the enstrophy, (c) the divergence error and (d) the time con-
straints ∆tmax(|ω|) (thin line) and ∆t|S|max (thick line).
T-D simulation with velocity deficit The simulation is initialized using

the vorticity field model consisting in the vorticity sheet (due to the lift) com-
bined with transverse vorticity (due to the friction drag). A periodic condition
is used in the axial (i.e. streamwise) direction. The periodic length is fixed
to 0.5 b; it has been shown that it is sufficient to capture the phenomena of
interest, see Cottin [19]. The time step is set to ∆t/t0 = 3.23 10−3. A uni-
form mesh size was used: h/b = 1/160. The flow was computed up to the
4.4. Results 81
time τ = t/t0 = 0.3. Those results have been extensively investigated in the
FAR-Wake technical report [19]. One should refer to this report for the related
results presented in this chapter.
4.4 Results
First, the vortex dynamics are qualitatively investigated. For each case, instan-
taneous 3-D surfaces and averaged 2-D slices are shown. Second, a complete
characterization of the resulting rolled-up wake is carried out; circulation pro-
files as well as diagnostics are provided.
4.4.1 3D dynamics: a qualitative analysis
In Fig. 4.5, one presents a 3-D isocontour of the vorticity field for the S-D
simulation without velocity deficit. The “classical” simple vortex-sheet rollup is
observed; i.e. without generation of turbulence by instabilities. For a CL /AR =
0.2, the rollup of the vortex system is almost complete at 8 · · · 10 wingspans
downstream; i.e. when the vortex spacing is equal to the vortex centroid b0 =
π
4 b. The flow is fully stationary.
Fig. 4.6 shows a 3D view of the λ2 iso-surface (see Jeong and Hussain [41])
of the S-D simulation with velocity deficit. Fig. 4.8.(a) and 4.8.(b) show the
related top-view and side-view. The λ2 criterion is a better vortical-structure
tracer than the vorticity field, especially for multi-resolution simulations. A
zoom in the near field is also presented in Fig. 4.7. A fast transition to tur-
bulence is observed within the first two wingspans downstream. In the central
part of the flow (i.e., in the low vortex-stretching region), jet-like instabili-
ties develop: very-short longitudinal wavelength instabilities and streamwise
Kelvin-Helmoltz instabilities. Those instabilities quickly grow, saturate and
eventually generate turbulence. In the meanwhile, a supplementary helical in-
stability develops in the primary vortex region; associated with a high vortex-
stretching. In this region, the vortex sheet, containing high transverse vorticity,
becomes unstable and is wrapped around the vortex core. Eventually, the flow
becomes fully turbulent. Moreover, the multi-resolution transition is quali-
tatively well performed: the small vortical structures are properly convected
through the common interface of two blocks with different resolution. Concern-
ing the T-D rollup with velocity deficit, Fig. 4.9 presents the 3-D vortex wake
using a vorticity iso-surface. The same global dynamics as the S-D case are
observed. For an exhaustive investigation of this case, refer to Cottin et al. [19].
Concerning the S-D simulation using a modeled wing force, it is worth

noting that the “wing” is half-formed at the station x/b = 0 because of the
symmetry of the external force with respect to x = 0. The initial time τ = 0 is
defined at the “trailing edge” arbitrarily fixed at x/b = 2 σx .
4.4.2 Quantitative analysis
Here, a quantitative analysis of the vortex-wake rollup is presented for each

case. Various slices at different downstream distances of the averaged stream-
wise vorticity are shown in Fig.4.10, 4.11, 4.12 for the S-D simulation with and
without velocity deficit and for the T-D simulation with velocity deficit respec-
tively. The spirally rolling up of the vortex-sheet is clearly observed especially
for the case without velocity deficit. The physics, associated to both cases with
velocity deficit, are similar. However, the vortex structures in the central part
decay much faster in the S-D simulation. Fig. 4.13 shows the related trans-
verse vorticity. We see the double shear layer emanating from the wing and
corresponding to the effect of two modeled boundary layers. Those are at the
origin of the transition to turbulence in the central region.
The averaged velocity deficit, uw , ux /U0 −1, is presented in Fig. 4.14, 4.15
and 4.16 for the SD simulation with and without velocity deficit effect and for
the T-D simulation with velocity deficit. The velocity deficit, occurring in the
S-D simulation without velocity deficit, is due to the combination of the S-D
effect associated to the spirally rolling up of the wake. This leads to a velocity
deficit in the vortex inner-core, which is relatively low ux (xc )/U0 − 1 ≃ −0.1
(where xc is center of the primary vortices) surrounded by a velocity surplus.
The global physics of the S-D case with velocity deficit are very different. As the
primary vortices form, a strong velocity deficit is entrained inside the vortices.
A velocity surplus (i.e. a “jet”) shortly occurs, between τ = 0.1 · · · 0.2, in the
vortex center. This “jet” component quickly diffuses, eventually vanishes and
4.4. Results 83
definitely gives way to a net significant deficit component ux /U0 − 1 ≃ −1.8.

A helical instability inside the vortices can occur in the condition to have
a low swirl number Swirl , max(uθ )/ max(uw ), where uθ is the tangential
velocity (see therein for more details). For a q-vortex model, the vortex becomes
unstable if Swirl . 1.8, consistent with the significant wake drag. Here, the
value is much higher ∼ 4 · · · 5, the vortices are thus stable.
The deficit in the vortex inner-core for the equivalent T-D simulation is
much different, see Fig. 4.16. The velocity deficit, in the vortex core as well
as in the central part, remains quite high ux /U0 − 1 ≃ −5. The related swirl
number is Swirl ≃ 1.3 · · · 1.8; a helical instability then develops, see Fig. 4.9
and Cottin et al. [19].
We define the turbulent kinetic energy (TKE), k, as follows
1 ′ 2
k, (ux ) + (u′y )2 + (u′z )2 ,
2
where u′α , uα − uα , see Fig. 4.17 and Fig. 4.18 show the TKE-planes for both
cases. The turbulence level is found to be of the same order of magnitude in
both cases.
Figure 4.5: Space-developing simulation of the vortex-sheet rollup without velocity deficit: 3D perspective view of two isocontours of the
vorticity norm kωkt0 = [200 400] colored with the axial vorticity component.
4.4. Results
Figure 4.6: Space-developing simulation of the vortex-sheet rollup with velocity deficit: 3D perspective view of the λ2 surface colored with the
85
streamwise vorticity component.
Figure 4.7: Space-developing simulation of the vortex-sheet rollup with velocity

deficit: 2 zooms in the near-field region of the 3D perspective view of the λ2 sur-
face colored with the streamwise vorticity component.
4.4. Results
(a)
(b)
Figure 4.8: Space-developing simulation of the vortex-sheet rollup with velocity deficit: top view (a) and side view (b) of the λ2 surface
colored with the streamwise vorticity component.
87
τ = 0.0
τ = 0.023
τ = 0.065
Figure 4.9: Time-developing simulation with velocity deficit: 3D view of two isocon-
tours of the vorticity norm kωkt0 = [200, 400] colored by the axial vorticity.
4.4. Results 89
τ = 0.13
τ = 0.20
τ = 0.26
Figure 4.9: Cont’d

Figure 4.10: Space-developing simulation with velocity deficit: slices of the time-
averaged streamwise vorticity ωx t0 ; from top to bottom, xb = [0, 0.1, 2, 4, 8], equiva-
lent time τ = [−0.003, 0.0, 0.07, 0.1, 0.3].
4.4. Results 91
Figure 4.10: cont’d.

Figure 4.11: Space-developing simulation without velocity deficit: slices of the

streamwise vorticity ωx t0 ; from top to bottom, xb = [0, 4, 8], equivalent time τ =
[0.0, 0.1, 0.3].
4.4. Results 93
Figure 4.12: Time-developing simulation with velocity deficit: longitudinally-

averaged streamwise vorticity ωx t0 ; from top to bottom, xb = [0, 4, 8], equivalenttime
τ = [0.0, 0.1, 0.3].
averaged spanwise vorticity ωy t0 ; from top to bottom, xb = [0, 0.1, 2, 4, 8], equivalent
time t = [−0.003, 0.0, 0.07, 0.1, 0.3].
4.4. Results 95

averaged streamwise velocity uw /V0 ; from top to bottom, xb = [0, 0.1, 2, 4, 8], equiva-
lent time τ = [−0.003, 0.0, 0.07, 0.1, 0.3].
4.4. Results 97

Figure 4.15: Space-developing simulation without velocity deficit: slices of the time-
averaged streamwise velocity uw /V0 ; from top to bottom, xb = [0, 4, 8], equivalent time
t = [0.0, 0.1, 0.3].
4.4. Results 99
Figure 4.16: Time-developing simulation with velocity deficit: longitudinally-

averaged streamwise velocity uw /V0 ; from top to bottom, time τ = [0.0, 0.1, 0.3].
Figure 4.17: Space-developing simulation with velocity deficit: slices of the tur-
bulent kinetic energy k/V02 ; from top to bottom, xb = [0, 4, 8], equivalent time
τ = [−0.003, 0.0, 0.07, 0.1, 0.3].
4.4. Results 101

Figure 4.18: Time-developing simulation with velocity deficit: turbulent kinetic en-
ergy k/V02 ; from top to bottom, time τ = [0.0, 0.1, 0.3].
4.4. Results 103
4.4.3 Characterization of the rolled-up vortex sheet in-

cluding the friction- and induced-drag effects
In this section, we characterize the net-result after rollup of the space-developing

vortex wake taking into account the friction drag. The circulation profile Γ(r):
Z 2π Z r
Γ(r) = ωx (r′ , θ)r′ dr′ dθ
0 0
and the azimuthal velocity profile uθ (r) = Γ(r)/(2πr) are defined. The “effec-
tive vortex core size, rc , is the radius where the azimuthal velocity is maximum,
uθ (rc ) = max (uθ (r)) = uθ,max .

r
The total vorticity of the half-plane, Γ+ , and the vortex-centroid spacing, yc ,

at a cross-plane x are also defined:
Z ∞ Z ∞
+
Γ (x) = ωx (x, y, z)dy dz
−∞ 0
and Z Z
∞ ∞
2
yc (x) = +
y ωx (x, y, z)dy dz
Γ (x) −∞ 0
respectively. The circulation profiles Γ(r) and the azimuthal velocity profiles
uθ (r) for the stations x/b = [2, 4, 6, 8] are presented in Fig. 4.19.(a) and
Fig. 4.19.(b). The maximum tangential velocity slightly decreases by diffusion
from 8.7V0 to 8.1V0 . Fig. 4.20 and Fig. 4.21 show the tangential velocity profiles
and the related circulation profiles for various fitted models at x/b = 4.0 and
x/b = 8.0, where the rollup is almost completed. All profiles have been fitted
on the tangential velocity profiles in the least-squared sense, without weight
function. Table 4.1 provides the model fitted parameters. We observe that the
one-length scale models (i.e. the Low Order Algebraic, the High Order Alge-
braic and the Gaussian) are not able to correctly capture either the circulation
profile or the tangential velocity profile. Nevertheless, the best one-length scale
model is the Low Order Algebraic one. The two-length scale models (i.e. the
Jacquin model and the Proctor-Winckelmans model) are able to match prop-
erly both the maximum tangential velocity and the core radius size rc . It has to
be stressed that the “global parameter” of the Jacquin model, which is usually
assumed to be “universal” (for wake at high Reynolds number and with thin
vortex sheet before rollup), was here modified to have a better fit. The param-
eter α has been taken as α = 0.5. The exponent p in the Protor-Winckelmans
model has been taken as p = 3, which is in the usual range p = 3 − 4. Finally,
the Regularized Kaden-Winckelmans model fits better than all one-length scale
models but doesn’t quite reach a high enough maximum tangential velocity. For
a long time, it seems to be close to the Low Order Algebraic model.
Fig. 4.22.(a) presents the evolution of the core size, rc . It grows from
rc /b = 1.5% to 5% at time t ≃ 3. One clearly observes the influence of the
x
resolution that is twice coarsened from b = 2.2, i.e. τ = 0.075. At this location,
rc /b = 2.4%: the vortex core is then captured using only 6 points, which is
still “acceptable” yet at the limit. The captured vortex-core is still fine enough
for the results to be of interest. In particular, it still compares well to those of
real aircraft system. Fig. 4.22.(b) shows the maximum velocity deficit inside
the vortex core. It confirms that the swirl number is ∼ 4.5.
One now recalls some of the assumptions made to end up with the classical
relation (originally developed using the inviscid theory) between the lift and the
wake circulation at a station x. It basically consists in realizing a momentum
balance around a control volume; as done in Chapter 2 using the formula by
Noca et al. [62, 63]. This formula determines the force, F, exerted by a fluid
on a body knowing the velocity field (and its derivatives) on the surface of a
control volume, SC .
Starting from that formula, if one assumes that the flow is stationary and
that the diffusive flux is negligible (as the Reynolds number is very high), this
formula reads
Z
F 1
= (u · u)n − (n · u)u dS
ρ SC 2
Z Z
1 1
− (n · u)(x × ω)dS + (n · ω)(x × u)dS (4.1)
2 SC 2 SC
we neglect the turbulent contributions. We also choose the control volume such
that the surface perpendicular to the flight path is fixed at the downstream
position x and all the other surfaces tend to infinity. Consequently, only one
4.4. Results 105
term, lying in this perpendicular surface, contributes to the lift (neglecting the
pressure contribution to the lift). Simplifying Eq.4.1 leads to
Z
L
≃ ux uz dS
ρ S(x)
Z
= (uw + U0 )uz dS
S(x)
Z Z
= U0 uz dS + uw uz dS .
S(x) S(x)
The second integral of the right hand side (RHS) corresponds to the contribu-
tion to the lift of the velocity deficit due to friction drag. Integrating by part
the first integral, leads to:
Z Z
L ∂uz
≃ −U0 y dS + uw uz dS .
ρ S(x) ∂y S(x)
∂uz ∂uy
Given that ωx = ∂y − ∂z :
Z Z Z
L ∂uy
≃ U0 y ωx dS − U0 y dS + uw uz dS .
ρ S(x) S(x) ∂z S(x)
∂uy
The second term vanishes as ∂z is left-right symmetric. Thus we finally obtain:
Z Z
L
≃ U0 y ωx dS + uw uz dS
ρ S(x) S(x)
Z
= U 0 y c Γ+ + uw uz dS .
S(x)
Γ+ yc
Fig. 4.22.(c) shows the evolution of those two terms. The former, Γ0 b , con-
stitutes most of the total lift since L/(ρΓ0 bU0 ) = π/4 ≃ 0.785. The latter
R
( Γ0 1bU0 uw uz dS) contributes to 0.5%. The turbulent terms of the momentum
R
balance, u′x u′z dS is also verified to be negligible. It is worth mentioning that
the way the forces (i.e. lift and drag) was imposed could drive the results.
Given the imposed lift force was stationary, the averaged field remains a good
image of a time-averaged quantity. Finally, the total circulation of the half-
plane, Γ+ /Γ0 , (defined above) is also presented in Fig. 4.22.(c) and is very well
conserved.
(a)
0.8 x/b=2.0
x/b=4.0
x/b=6.0
Γ(r)/Γ0
0.6
x/b=8.0
0.4
0.2
0
0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4
r/b
(b)
14
x/b=2.0
x/b=4.0
12
x/b=6.0
x/b=8.0
10
uθ (r)/V0
0
0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5
r/b
Figure 4.19: Circulation profiles (a) and tangential velocity profiles (b) at various
distances in the wake.
4.4. Results 107
(a)
12
Low Order Algebraic
High Order Algebraic
10 Lamb−Oseen (Gaussian)
Jacquin
Kaden−Winckelmans
8 Proctor−Winckelmans
Space−Developing LES
uθ (r)/V0
0
0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5
r/b
(b) 9
Low Order Algebraic
8 High Order Algebraic
Lamb−Oseen (Gaussian)
7 Jacquin
Kaden−Winckelmans
6 Proctor−Winckelmans
uθ (r)/V0
0
0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5
r/b
Figure 4.20: Tangential velocity profiles for various vortex models (a) at x/b = 4.0
(τ = 0.13) and (b) at x/b = 8.0 (τ = 0.26).
(a)
0.8
0.6
Γ(r)/Γ0
0.4 Low Order Algebraic

0.2 Jacquin
Kaden−Winckelmans
Proctor−Winckelmans
0
0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5
r/b
(b)
1
0.8
0.6
Γ(r)/Γ0
0.4 Low Order Algebraic

0.2 Jacquin
Kaden−Winckelmans
Proctor−Winckelmans
0
0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5
r/b
Figure 4.21: Circulation profiles for various vortex models (a) at x/b = 4.0 (τ =
0.13) and (b) at x/b = 8.0 (τ = 0.26).
4.4. Results
Table 4.1: List of usual vortex models.The values of the model parameter(s) fitted on the space-developing rollup simulation at x/b = 4.0 and
x/b = 8.0 are provided
Γ(r)/Γ0 a1 /b a2 /b other
(r/a1 )2
Low Order (r/a1 )2 +1 0.042 rc = a1
Algebraic 0.048
(r/a1 )2 ((r/a1 )2 +2γ)

High Order ((r/a1 )2 +γ)2
0.048 rc = a1
Algebraic 0.056 γ = 1.781

Gaussian 1 − exp − β(r/a1 )2 0.056 rc = a1
0.064 β = 1.256
2
Jacquin (r/a 2)
(1+α)/4 (1−α)/4 0.022 0.19 α = 0.5
(a1 /a2 )4 +(r/a2 )4 1+(r/a2 )4
0.075 0.17
!
(r/a1 )2
Proctor- 1 − exp − p 1/p 0.037 0.072 rc ≃ (a81 /a32 )1/5
(r/a1 )2
1+
(r/a2 )3/4
Winckelmans 0.049 0.065 p=3
109
α(r/a2 )1/2
Regularized 1 − exp − (r/a1 )3/2 1+(α−1)(r/a 2 ))
0.020 0.46 α = 1.78
Kaden-Winckelmans 0.030 0.40
4.5 Conclusion
This chapter presents the investigation of a velocity deficit effect (due to the
friction drag) on the rollup of a thin vortex-sheet at very high Reynolds num-
ber (Re = 106 ). Large eddy simulations were carried out and three cases was
investigated: two space-developing (S-D) simulations, with and without this ve-
locity deficit, and one time-developing (T-D) simulation with velocity deficit.
The first S-D simulation uses the multi-resolution ability developed and de-
scribed in Chapter 3 while the other two simulations use a uniform resolution.
A “in cruise condition” is here considered. This chapter is the complementary
follow-up of a Far-WAKE technical report (Cottin et al. [19]).
First, it is worth noting that the refined mesh handling technique has been
applied using large eddy simulation at high Reynolds number and using mul-
tiscale subgrid modeling.
Second, it has been shown that the velocity has a major influence on the
resulting vortex-system. The flow is intrinsically unstable and rapidly becomes
turbulent. A strong velocity deficit remains inside the vortex-core of order
−[1.5 · · · 1.9]V0 (towards the wing).
The wake-rollup with velocity deficit using the T-D and the S-D simulations
have been also compared. The global dynamics are similar. In the T-D case,
the velocity deficit inside the vortex-core is higher, leading to the development
of a helical instability.
Finally, the net result of the S-D wake rollup with velocity deficit has
been characterized. Most of vortex models have been fitted for two station
x/b = [4, 8]. It appears that the two-scales models (Proctor-Winckelmans and
Jacquin) are both good candidates to model wake vortices. Indeed, two length
scales are required to both fit the vortex-core size and the maximal azimuthal
velocity.
4.5. Conclusion 111
(a) 0.05
0.045
0.04
0.035
rc /b
0.03
0.025
0.02
0.015
0.01
0 0.05 0.1 0.15 0.2 0.25
τ
(b) 0
−1
−2
uw /V0
−3
−4
−5
0 0.05 0.1 0.15 0.2 0.25
τ
(c)
1
0.8
0.6
[−]
0.4
0.2
0
0 0.05 0.1 0.15 0.2 0.25
τ
Figure 4.22: (a) Evolution of the vortex effective core size rc . (b) Evolution of the
maximum velocity deficit inside the vortex core and (c) Evolution of the linear-impulse
yc Γ+ /(Γ0 b) (solid line), of +
R the total circulation of the half-plane Γ /Γ0 (dashed line)
1
and of the integral Γ0 bU0 uw uzdS (dotted line).
Chapter 5
Numerical investigations of
end-effects associated with
accelerated/decelerated
wings : time-developing
and space-developing
simulations
This chapter was written within the FAR-Wake project as a technical report.
European project AST4-CT-2005-012238, Report TR1.1.2-6, February, 2007.
The authors are T. Lonfils, R. Cocle, G. Daeninck, C. Cottin and G. Winck-
elmans.
Abstract In this report, we investigate the end-effect phenomenon occurring

due to the variation of the circulation in a vortex wake. We first focus on com-
paring a numerical simulation to an experiment made by CNRS/IRPHE[57].
Two types of instabilities occur : an axisymmetric one produced by the circu-
114 Chapter 5. Numerical investigations of end-effect
lation variation, and generating an axial velocity which can lead to a helical
one, as it is presently the case, if this velocity is strong enough. We also
compare some results to the theory derived for the Lamb-Oseen vortex. We
then investigate a more realistic case of end-effects associated with an acceler-
ated/decelerated wing. Analyzing the corresponding time-developing and space-
developing (S-D) simulations, an explanation of the different behavior in the
accelerated and in the decelerated cases is proposed.
5.1 Introduction
Due to the aircraft lift, wake vortices are inevitably generated. They can be
hazardous for a follower aircraft and are then an important issue to ensure a
high security level in air traffic. The purpose of this work is to characterize
the typical vortex instabilities occurring when aircraft accelerate/decelerate,
which could lead to the disappearance of the vortex coherence : this is why
such instabilities are of interest. Theoretical and numerical investigations have
already been carried out, but they remain academic as they are applied to
the Lamb-Oseen vortex (see [40, 30, 59, 29, 39]). This phenomenon has been
extensively investigated by Leibovich [48, 49]. See also the review done by
Lucca-Negro and O’Doherty [53].
In the first part of this report, we perform a numerical simulation to be
compared to an experiment also carried in the framework of FAR-Wake [57] in
which the axial propagation of vortex perturbations was investigated.
Two other simulations are also carried out. The first one is time-developing,
and thus assumes periodicity. The second one simulates the complete spatial
evolution of the wake, which evolves using an inflow condition. The latter
constitutes a better representation of the reality, but it requires much more
computational resources. Both cases are analyzed and compared.
Experimental investigations of end-effect generated by a accelerated/dece-
lerated wing have also been carried out using the UCL/TERM facilities ([27],
[26]). The influence of the lift coefficient, the flight speed and the acceler-
ation/deceleration instensity have also been studied. They confirm that the
propagation speed is of the order the maximal tangential velocity. Fig. 5.1
shows a helical instability propagation for a deceleration case.
5.1. Introduction 115
Because vortex lines have to be closed, the increase of the wake circulation
is associated with the creation of spanwise vorticity component. The general
description of the corresponding wake configuration is depicted in Fig. 5.2. The
global dynamics are essentially dependent on the ratio of the net increase of
circulation (i.e. ∆Γ = (Γ1 − Γ0 )/Γ1 ) over the acceleration distance E/b0 ([27]).
The so-called “end-effect” phenomenon (relative to a sudden acceleration/de-
celeration of the wing) corresponds to the particular case with Γ0 = 0; but the
term will here be also used to qualify the general case (i.e. Γ0 6= 0).
Figure 5.1: Dye visualization of instabilities generated by a decelerated wing. The

fluorescein is injected at the wingtips (for more details, see [27]).
5.2. Numerical method 117
5.2 Numerical method
For the simulations, we use the Vortex-in-Cell method combined with the Par-
allel Fast Multipole Method (called VIC-PFM, for short). It is based on a
combination of Lagrangian and finite difference methodologies, and it solves
the vorticity form of the incompressible Navier-Stokes equations
Dω
= ∇ · ( uω ) + ∇ · W , (5.1)
Dt
in which D/Dt = ∂/∂t + (u · ∇) . The divergence of the diffusive vorticity flux

tensor W is written here

T
∇ · W = ν ∇2 ω + ∇ · νsgs ∇ω s + (∇ω s ) (5.2)
where νsgs is the effective subgrid-scale viscosity and ωs is the “small-scale”

part of the LES field (since we here use the multiscale subgrid-scale approach).
The numerical solution of Equation (5.1) is sought in a two-fold approach.

First, the discrete vorticity field is constructed from a sum of regularized vortex
particles of strength α such that

X 1 1 |x − xp (t)|2
ωσ (x, t) = αp (t) √ exp − , (5.3)
p ( 2π σp )3 2 σp2
R
in which αp = ω dx = ωp h3 , xp is the particle position and σp its regu-
larization parameter [16, 87]. The convective part, or lhs of Equation (5.1),
is evaluated using a Lagrangian approach with d xp /dt = u(xp ) . For the sec-
ond part, the time variation of the vortex particle strength, i.e., the rhs of
Equation (5.1) that includes both the vortex stretching and dissipation terms,
it is solved on a regular grid using 2nd order finite differences. Interpolation
from the Lagrangian particles to the Eulerian grid is done using the M4′ scheme
[16, 87]. Once on the Eulerian grid, the stream function vector ψ is evaluated
by solving the Poisson equation
∇2 ψ = −ω . (5.4)
The velocity field, needed for convection and stretching, is then obtained by
evaluating u = ∇ × ψ using again 2nd order finite differences. The global

time marching procedure is carried out using a 2nd order Leap-frog scheme for
the convection and the 2nd order Adams-Bashford scheme for the stretching
and dissipation terms. Finally, the divergence-free character of the vorticity
vector field is insured by a proper reprojection of the vorticity field (which also
requires solving a Poisson equation).
Particular to the present implementation of this procedure is the treatment

of the boundary conditions. The idea is here to work on a domain that contains
tightly the vorticity field and that can be decomposed in several subdomains on
which the exact boundary conditions are obtained using a Parallel Fast Mul-
tipole (PFM) method. This amounts solving a 3-D Poisson equation without
any iteration between the subdomains (e.g., no Schwarz iteration required):
this is so because the PFM method has a global view of the entire vorticity
field. So, at any given time step, the boundary conditions for ψ are obtained
using the Green’s function approach via the PFM method. This allows for a
significant reduction of the computational cost for solving Equation (5.4) —as
it retains proper open domain conditions and the capability to make parallel
subdomain decomposition simulations— when comparing with more classical
VIC methodologies [89, 13].
The LES modeling is here done using a multiscale subgrid-scale model: the
model acts on the small scale part of LES field, ω s , only. It is obtained from
the complete LES field, ω, using a compact discrete filter (applied iteratively).
We here use the Regularized Variational Multiscale (RVM) model [12]. The
advantage of such model is that it preserves the inertial range while providing
dissipation at the high wave numbers: it is only active during the complex
phases of the flow, while remaining inactive during the phases when the vortices
are coherent and well-resolved.
5.3. Analysis and comparison 119
5.3 Numerical simulation of end-effects within

vortex pairs : analysis and comparison with
the IRPHE experiments.
Experimental setup
In this section, we present the numerical investigation carried out in order to
compare with one of the CNRS/IRPHE experiments. In their experiments, the
vortex is generated in a water tank using a plate which is rotated (for more
details, see [57]). We here focused on the case at moderate Reynolds number,
i.e. ReΓ = 3960, in order to be able to afford an equivalent Direct Numerical
Simulation (DNS).
Numerical setup & initial condition

We investigated the case of a vortex pair spaced by b0 and linked at the ex-
tremities as depicted in Fig. 5.2 and Fig. 5.3.
Figure 5.2: Schematic of the wake vortex lines corresponding to an accelerated and
decelerated wing. Flow periodicity is here assumed over the length L0 + L1 .
This case is also quite similar to the wake generated by a wing suddenly
accelerated and decelerated (i.e. Γ0 = 0), assuming that the vortex pair is
immediately formed. The vortices are initialized using the Low Order Algebraic
(LOA) vortex model :
Γ(r) r2 Γ(r)
= 2 uθ (r) = , (5.5)
Γ1 r + rc2 2πr
where rc is the effective vortex core size (defined as the radius where the tan-
gential velocity uθ (r) is maximum). Here we use rc = 0.05 b0. We thus initially
(a) êy
6
-êx
(b)
êz
êy êx
Figure 5.3: (a) 2-D cross section of the initial vorticity field (vortex lines) at z/b0 =
0.0. (b) 3-D view of the vorticity norm isocontour |ω|b20 /Γ1 = 1.0 and |ω|b20 /Γ1 = 5.0.
have
Γ1 1 Γ1 rc 1 rc
uθ,max = = / = V0 / = 10V0 , (5.6)
4πrc 2 2πb0 b0 2 b0
where V0 is the wake descent velocity. The Reynolds number, based on the cir-
culation, is here chosen as Re = Γ1 /ν = 5000. The distance of the reconnection
E and the vortex length L are respectively equal to 1.0 b0 and 20 b0 . The time
step and the mesh size are fixed to ∆t Vb00 = 3.98 10−4 and h/b0 = 1/64. The
simulation was carried out using 60 million grid points. It ran on 16 processors
during 240 hours.
Analysis of the results
Since the problem is left/right symmetric and begin/end symmetric, we focus

on only one vortex, without loss of generality.
Evolution of the vortex core size The instantaneous circulation profile

Γ(r) is defined by
Z 2π Z r
Γ(r) = ωx (r′ , θ) r′ dr′ dθ. (5.7)
0 0
It is well known that the end-effect phenomenon much increases the effective
vortex core. Moreover, in the current simulation, the Reynolds number is
moderate and thus the increase of the vortex core by diffusion is also not
negligible. Using the analogy of the Gaussian vortex diffusion, one characterizes
the evolution of the LOA vortex core by diffusion as
rc2 (t) ≃ rc2 (0) + ανt, (5.8)
where α is found to be equal to roughly 7.9 (using the current simulation, we

analyzed the evolution of the vortex core at the cross plane x/b0 = 10 when
it is not yet affected by the end-effect) as shown in Fig. 5.4.(a), instead of 4.0
for the Gaussian vortex. We can also estimate a priori the evolution of the
maximum tangential velocity, assuming that the vortex remains close to the
LOA model :
Γ1
uLOA
θ,max (t) = . (5.9)
4πrc (t)
Fig. 5.4.(b) shows the evolution of the measured maximum tangential velocity
as based on the circulation profile assuming a LOA model. It turns out that
the maximum tangential velocity decreases significantly during the simulation,
and that both cases are quite different. The case based on the LOA model
underestimates the velocity during the whole simulation : this means that the
vortex structure doesn’t remain like the LOA one.
Fig. 5.5 presents the evolution of the vortex core size rc for three cross
planes : x/b0 = 5.0, x/b0 = 7.5 and x/b0 = 10.0. Each shows the same
global evolution. At first, the vortex core increases purely by diffusion and
then suddenly magnifies when the perturbations cross the plane. We can then
calculate the propagation speed of the instability as the ratio of the distance
between two cross planes, to the time difference when the vortex core increases
:
uc = ∆x/∆t.
Here we find uc /uθ,max(0) ≃ 2.5/2.8 = 0.89.
(a) (b)
−3
x 10
14 1
0.9
12 0.8
uθ,max /uθ,max(0)
10 0.7
0.6
rc2 /b20
8 0.5
0.4
6 0.3
4 0.2
rc2 (t)/b20 = rc2 (0)/b20 + 7.9νt/b20 0.1
2 0
0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1
2πb2 2πb2
t Γ1 0 t Γ1 0
Figure 5.4: (a) Time evolution of the effective core size rc (thick line) and its
least-square fitted line (thin line). (b) Time evolution of the maximum tangential
velocity uθ,max (t) (solid line), and also that assuming that the LOA model remains
valid (dashed line).
We also analyze the potential hazard of the vortex. The total circulation
of the vortex Γ1 is not a good estimation since it is conserved. We use Γ5−15 :
the mean of the vortex circulation between 5 and 15 meter for an aircraft with
0.25
∆tuθ,max(0)/b0 = 2.8
0.2
rc /b0
0.15
0.1
0.05
0 2 4 6 8 10 12 14
Γ1
tuθ,max(0)/b0 = 10t 2πb 2
0
Figure 5.5: Evolution of the vortex core size rc for various cross plane : x/b0 = 2.5
(dash-dotted line), x/b0 = 5.0 (dashed), x/b0 = 7.5 (dotted line) and x/b0 = 10.0
(solid line).
1
0.9
0.8
0.7
Γ5−15 /Γ1
0.6
0.5
0.4
0.3
0.2
0.1
0
0 2 4 6 8 10 12 14
Γ1
0
Figure 5.6: Circulation Γ5−15 as a function of time for the cross plane x/b0 = 5.0.
a 60 meter wingspan. It reads

Z b/4
1
Γ5−15 = Γ(r) dr, (5.10)
b/6 b/12
where b is here taken as b = (4/π)b0 (assuming a wing with elliptical load-

ing). Fig. 5.6 shows its evolution : we clearly observe a net diminution at
tuθ,max/b0 ≃ 4.8, due to the passage of the perturbations.
Propagation of the perturbations There are several manners to measure

the speed of the perturbation. The conclusion of the observations of Meunier in
its experiment[57] is that there would be two main instabilities. The first one
is supposed to be a “front wave” with a strong axial velocity within the vortex
core (as shown in Fig. 5.7.a); the second would be a helical instability which
strongly deforms the vortex (as shown in Fig. 5.7.b) and would propagate at a
lower speed .
(a)
(b)
Figure 5.7: (a) Longitudinal section of the vortex at y = b0 /2. Axial component of
the velocity field and front velocity the non-zero axial velocity (uw ). (b) Side view of
the vorticity norm isocontours ||ω|| and description of the isophase point speed (uip )
and its front speed(uif ).
Fig. 5.8 shows the time evolution of the front position of the non-zero axial
velocity. For the nondimensionalization, we used the characteristic length scale
b0 and the maximum tangential velocity as the velocity scale. We compared
the propagation speed using three different tangential velocity : (a) the time
evolving maximum tangential velocity uθ,max(t) and (b) the initial maximum
tangential velocity uθ,max(0). As foreseen, the propagation speed is of the or-
der of the maximal tangential velocity in each case : 1.5 and 0.86 respectively.
The viscous diffusion of the vortex core here plays a significant role in slowing
down the propagation, as the tangential velocity also diminishes. We can also
compare the propagation speed to the theoretical one predicting the propaga-
tion of the group velocity of an axisymmetric mode for a Gaussian vortex (see
Fabre et al.[29] and Jacquin et al.[40]) :
Γ1
uw ≃ 0.63 (5.11)
2πσ
where σ is the Gaussian vortex core size. To compare to that, we here consider
the Gaussian vortex which has the same long wavelength dynamics as the LOA
Γ1
vortex : this give σeq = 1.34 rc . Here we find uw /( 2πσeq
) ≃ 0.57 indeed close
to the theoretical value. In the experiment, the wave perturbation speed was
above : roughly 0.9.
Another important information is the swirl parameter S which is the ratio

between the maximum tangential velocity and the maximum axial velocity in
the vortex. Here we find S being nearly equal to 0.9 . . . 1.2. A value included
in the interval [0, 1.5] means that the helical Kelvin mode of the vortex is
unstable.
(a) (b)
10 10
8 8
6 6
x/b0
x/b0
4 4
slope : 0.86
2 slope : 1.5 2
0 0
0 2 4 6 0 5 10
tuθ,max(t)/b0 tuθ,max (0)/b0
Figure 5.8: Position of the non-zero axial velocity front (solid) and its least-squared
fitted line (dashed) for various time scaling based on : (a) the time evolving maximum
tangential velocity uθ,max (t), (b) the initial maximum tangential velocity uθ,max (0).
We then analyze the propagation of this helical perturbation. The time

evolution of the isophase points (as described in Fig. 5.7.b with the red lines)
and the position of their front are presented in Fig. 5.9. We find a group ve-
locity equal to uif /uθ,max(0) = 0.85. This velocity is the same as the wave
speed uw (which creates the axial velocity). The reason could be the fact that
the axial velocity generated by the wave is strong enough to permit the helical
instability to immediately occur. This is the main difference with the exper-
imental results for which the wave speed was twice the speed of the isophase
front. Meunier suggests that the speed of the isophase point compares with
the phase velocity of the helical Kelvin mode (m = 1). We find in our case
a perturbation wavelength of λ/rc = 16 (taken where the instability is inside
the vortex core) and thus a wavenumber k rc = 0.39. The theoretical phase
velocity of the helical mode of a Lamb-Oseen vortex with this wavenumber is
Γ1
0.33 . (5.12)
2πσ
In our case we find uip ≃ 0.23 Γ1/(2πσeq ).
uif /uθ,max
5 (0) ∼ 0.85 uip /uθ,max(0) ∼ 0.21
x/b0
4 uip /uθ,max(0) ∼ 0.15

3 uip /uθ,max(0) ∼ 0.12
2 uip /uθ,max(0) ∼ 0.10
1 uip /uθ,max(0) ∼ 0.041
0
0 2 4 6 8
Γ1 10 12
0
Figure 5.9: Evolution of the deformation front (dash) and of the isophase points
(dash-dot). The solid lines correspond to the least-squared fitted lines.
Visualization of the perturbation We define the vorticity perturbation

field as
ω = ω (0) + ω ′ , (5.13)
where ω(0) is the “base” field (i.e. the vorticity field that would exist without
the perturbation). Here, ω (0) is taken as the field at the center slice, i.e. at
x/b0 = 10. ω ′ is thus the vorticity perturbation field. Isocontours of ωx′ for
two different times are presented in Fig. 5.10.
Fig. 5.11 shows the comparison between the numerical simulation and the
experiment for the axial vorticity component at a fixed slice for various times.
The global dynamics is roughly the same as in the experiment. The first phase
consists in the simple laminar diffusion of the vortex. The effective vortex core
then grows suddenly when the wave crosses the slice, with loss of coherence
after that.
3-D view of vorticity norm isocontours are also presented in Fig. 5.12, clearly
slowing the helical instabilities.
−0.3
−0.9
−0.4
−1
−0.5
−1.1
z/b0
z/b0
−0.6
−1.2
−0.7
−1.3
−0.8
−1.4
−0.9
−1.5
0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.2 0.3 0.4 0.5 0.6 0.7 0.8
y/b0 y/b0
Figure 5.10: Isocontours of the vorticity perturbation axial component ωx′ b20 /Γ1 at
x/b0 = 5.0 at tuθ,max /b0 = 6.0 (left) and 12 (right) : ten equally spaced levels from
−4.5 to 4.5; positive levels are in red and negative levels are in blue. The vortex center
and the effective core radius rc are also shown.
Figure 5.11: Visualization of the vorticity field axial component for the cross plane
x/b0 = 5.0. Numerical simulation (left) and IRPHE experiment (right). The size
of the visualization boxes are 15.5 rc and 11.1 rc respectively. The time between two
plots is ∆tuθ,max /rc ≃ 24. For the numerical simulation, a circle using a radius rc (t)
centered on the vorticity centroid is also plotted.

Γ1
tuθ,max (0)/b0 = 0.0 t 2πb 2 = 0
0

Γ1
tuθ,max (0)/b0 = 1.6 t 2πb 2 = 0.16
0

Γ1
tuθ,max (0)/b0 = 3.2 t 2πb 2 = 0.32
0
Figure 5.12: Visualization of the vorticity norm isosurfaces for ||ω||b20 /Γ1 = 1.0 (low
opacity), 5.0 (medium opacity) and 10 (high opacity), colored with the axial vorticity
component.

Γ1
tuθ,max /b0 = 4.8 t 2πb 2 = 0.48
0

Γ1
tuθ,max /b0 = 6.4 t 2πb 2 = 0.64
0

Γ1
tuθ,max /b0 = 8.0 t 2πb 2 = 0.080
0
Figure 5.12: [Cont’d]


Γ1
tuθ,max /b0 = 9.5 t 2πb 2 = 0.95
0

Γ1
tuθ,max /b0 = 11 t 2πb 2 = 0.11
0

Γ1
tuθ,max /b0 = 13 t 2πb 2 = 0.13
0

5.4 Time-developing and space-developing sim-

ulations of an accelerated/decelerated wing.
We now investigate the phenomenon of the acceleration/deceleration of a wing.

The behavior of the flow over the wing will not be captured by the simulations
but modeled using the lifting line approach. This theory (i) neglects the effects
of the boundary layers and (ii) “replaces” the wing by a vortex line which has
a certain circulation profile Γ(y) along the spanwise direction. The lift profile
and the circulation profile are related by l(y) = ρU Γ(y) (where U is the flight
speed and ρ the fluid density.). Since the vorticity field is divergence free, a
vortex sheet of circulation per unit length γ(y) has to be shed such that
dΓ
γ(y) = − (y) (5.14)
dy
We can then express the total lift of the wing as

Z b/2 Z b/2
L= l(y) dy = ρ U Γ(y)d y, (5.15)
−b/2 −b/2
We define Γ0 and b0 , respectively the total circulation of the half plane and the
spacing between the vorticity centroids (and thus also the vortex spacing after
rollup). Then
Z b/2
Γ(y)dy = Γ0 b0 (5.16)
−b/2
and
L = ρ U Γ0 b 0 . (5.17)
Let’s now consider an elliptic wing which has the following characteristics :
• wingspan b
• Surface S
• Aspect ratio Ar , b2 /S
The lift can also be related to the aerodynamic characteristic of the wing CL ,
L/( 21 ρU 2 S). We can thus express the characteristics of the wake (Γ0 and b0 )
5.4. T-D and S-D simulations of an accelerated wing 133
as a function of the wing characteristics :
1 1 CL
Γ0 b 0 = U SCL = U b2 (5.18)
2 2 Ar
Here we use
CL
CL = 1.5, AR = 7.5, thus = 0.20, Re = Γ0 /ν = 104 .
AR
The generated vortex sheet rolls up and eventually forms two counter-rotating
wake vortices. In the previous section the initial condition was two vortex tubes
reconnected at the extremities. We don’t do it anymore. Indeed, we have to
compare the relative characteristic time of each phenomenon to be sure that the
time/space-evolution of the wake is properly investigated. The characteristic
time of the vortex wake t0 = b0 /V0 has to be compared to the characteristic
time of the perturbation propagation tp = b0 /Vp . Therefore,
Vp
t0 /tp =
V0
uθ,max
≃
V0
where Vp is the propagation speed of the perturbation which is of the order

of the maximum tangential velocity uθ,max of the vortex; a typical value is
uθ,max/V0 ≃ 10. It means that we have to take into account the phenomenon
of rollup.
Moreover, the third phenomenon is the flight itself which has a characteristic
b0
time tf = U. The ratio of the latter characteristic tf and tp is
Vp Vp V0
tf /tp = =
U V0 U
uθ,max 4 CL
≃
V0 π 3 Ar
∽ 0.2 . . . 0.3
It means that, a priori, the end-effect of the vortex will clearly act in the wake
region and won’t be affected by the space-developing evolution of the wake.
However, this “global” analysis neglected the influence of an axial velocity in
the vortex core region due to the 3-D space-developing rollup (see [59] and
[30]). Hence, we also investigate the space-developing case.

So, we propose to investigate and to compare the two cases : (i) the
space-developing (S-D) simulation and (ii) the time-developing(T-D) simula-
tion. Fig. 5.13 presents the comparison between those two concepts.
(i) The S-D simulation consists in three phases. We first simulate the com-
plete space evolving rollup of the vortex sheet shed by a wing flying at
a constant speed U0 . This is the initial condition. During the second
phase, the wing accelerates to the speed U1 . And, finally at the time of
the third phase, the wing flies at the constant speed U1 . The instabil-
ity evolves timely and is also convected in the wake moving towards the
higher circulation region. Since there are an inflow and an outflow, the
period of validity, and thus of analysis, is limited by the time needed to
convect the instability out of the simulation domain.
(ii) The T-D simulation assumes that the whole space evolution of the wake,
including the three phases, is initially generated. We then use periodic
boundary conditions in the flight direction. The initial condition is thus a
vortex sheet over a distance L0 (corresponding to a wing flying at a speed
U0 ), an acceleration to the speed U1 over a distance E, and another piece
of vortex sheet over the distance L1 corresponding to the one of the flight
speed U1 . To match with the periodic boundary conditions, we also add
a deceleration phase (also over distance a E).
The vortex sheet generated by the wing has to be regularized in order to

produce a regular vorticity field to be provided to the VIC-PFM code. For
that, we use convolution with a Gaussian :
Z b/2
1 (y − y ′ )2 + z 2
ωx (y, z) = exp − γ(y ′ )dy ′ , (5.19)
πσ 2 −b/2 σ2
where σ is the regularization parameter of the vortex sheet. For an elliptic

wing, the span loading (i.e. circulation profile) also has an elliptic profile
!1/2
2 U CL b 2
Γ(y) = 1 − y/(b/2) (5.20)
π A
| {z r }
Γ0
Γ1
t 2πb 2 = 0.0 Lx
0
S-D :
T-D :
L0 L1
E
Γ1
xref
t 2πb 2 = 0.052
0
S-D :
T-D :
êx êz
?
êy
Γ1
t 2πb 2 = 0.34
0
S-D :
T-D :
Γ1
t 2πb 2 = 0.68
0
S-D :
T-D :
Figure 5.13: Comparison of time-developing (T-D) and space-developing (S-D) sim-

ulations: top view of the vorticity field at various times. The zone in gray corresponds
to an example of zoom boxes, moving in the S-D case and fixed in the T-D case, used
for the 3-D visualizations in Fig. 5.22, 5.23 and Fig. 5.24.
Notice that γ(y) becomes infinite at y = ±b/2, but not Γ(y). To avoid numer-
ical problems, one rewrites Eq.5.19 using integration by part:
Z b/2
2 (y − y ′ )2 + z 2
ωx (y, z) = exp − (y − y ′ )Γ(y ′ )dy ′ . (5.21)
πσ 4 −b/2 σ 2
Numerical setup & initial condition for the space-develo-

ping simulation
We investigate the S-D behaviour of the wing acceleration/deceleration with

a constant lift coefficient CL . It means that the lift increases/decreases: for
instance, it is usually the case for an experiment in a towing tank. The velocity
of the wing will follow the expression:

t − ta t − ta
U (t) = U0 + a ∆ta g = U0 + ∆U g , (5.22)
∆ta ∆ta
where U0 , a, ta , ∆ta and ∆U = U1 − U0 are respectively the wing speed before

the acceleration, the acceleration itself, the time when the acceleration starts,
the duration of the acceleration and the total velocity variation. The function
g(ζ) is defined as followed :




 0 if ζ ≤ 0

g(ζ) = ζ if 0 < ζ ≤ 1 (5.23)




 1 if ζ > 1
The wing motion is composed of three phases. The first phase consists in
the flight at a constant velocity U0 , then followed by an acceleration or an
deceleration phase achieved in a time ∆ta and the third one is the flight at a
higher/lower velocity U1 .
The wake circulation increases following U (t):

∆U t − ta
Γ(y, t) = Γ(y) 1 + g . (5.24)
U0 ∆ta
Here we use an acceleration since the instability propagates from the low cir-
culation Γ0 region to the higher circulation Γ1 region: we can thus observe the
end-effect longer due to the inflow/outflow.
The speed of the wing varies from U0 to U1 over a distance E. We thus have:
2E
∆ta = (5.25)
2U0 + ∆U
The S-D simulation considered here is, at first, a complete rollup of the vortex
sheet shed by a wing having an elliptic span loading and flying at velocity U0 .
The wing then accelerates to the velocity U1 over a distance E = b/2. The
component of the vorticity field in the flight direction (êx , the axial direction)
also increases :

∆U t − Ta
ωx (y, z, t) = ωx (y, z) 1 + g .
U0 ∆Ta
A component of the vorticity field in the spanwise direction (êy ) is then also
produced corresponding to the temporal variation of the circulation. Starting
using the total derivative of ωy (y, z, t), we can write
Z Z Z
∂ωy ∂ωy ∂ωy
ωy (y, z, t) = dy + dz + dt. (5.26)
∂y ∂z ∂t
Since the vorticity field is divergence free and since the vertical vorticity com-
ponent is zero, this gives :
Z Z Z
∂ωx ∂ωy ∂ωy
ωy (y, z, t) = − dy + dz + dt. (5.27)
∂x ∂z ∂t
One further assumes that the spanwise component of the vorticity due to the
∂ωy
acceleration is uniform in the vertical direction ( ∂z = 0) and over the time
∂ω
step dt ( ∂ty = 0). One then obtains:
Z
∂ωx
ωy (y, z, t) ≃ − dy
∂x
Z
dt ∂ωx
= − dy
dx ∂t
Z y
1 ∆U ′ t − ta
= − g ωx (y ′ , z)dy ′ .
U (t) ∆ta U0 ∆ta 0
We define the dimensionless time with the characteristic time of the vortex pair
2πb20
during the high velocity regime: t0 = Γ1 (region towards which the instability
propagates).
A computational region of interest for the T-D simulation corresponds to

a region of the S-D simulation moving at U1 as sketched in Fig.5.13. To make
their comparison, we follow a reference cross plane xref : moving in the S-D
case and fixed in the T-D case.
The time step of the simulation is fixed to ∆t/t0 = 6.45 10−5. The mesh size
is h/b = 1/100 and the regularization parameter is σ = 2h. The length of the
computational box has been chosen long enough to observe the phenomenon :
Lx /b = 15. The simulation required 20 processors for 13, 500 time steps (run
for ∼ 125 hours). We carried out the simulation using ∼ 40M grid points.
Numerical setup & initial condition for the time-develo-

ping simulation
The initial condition of the T-D simulation is in the same spirit of the S-D
simulation. Over the distance E (also fixed to b/2), the circulation of the
half vortex sheet increases from Γ0 to Γ1 in space (and not in time as in the
space-developing simulation) :

∆U x
ωx (x, y, z) = ωx (y, z) 1 + g (5.28)
U0 E
To make the problem periodic, we finally use :

∆U x L1 − x
ωx (x, y, z) = ωσ (y, z) 1 + g g . (5.29)
U0 E E
Since the vorticity field is divergence free and since the vertical component of
the vorticity is zero, the spanwise vorticity component is obtained from :
Z
∂ωx
ωy (x, y, z) = − dy., (5.30)
∂x
leading to
′ Z y
∆U x L1 − x
ωy (x, y, z) = − g g ωx (y ′ , z)dy ′ . (5.31)
U0 E E 0
Fig. 5.14 presents the different time-evolution of the wake circulation for the
S-D and for the T-D simulations which corresponds to the wake circulation of
Γ1
Γ(t)
Γ0
0
0
ta ta + ∆ta
t
Figure 5.14: Schematic of the wake circulation evolution for the S-D simulation
(solid) and for the T-D simulation (dash).
a plane moving in the initial condition at the wing speed. The inital condition
for the T-D case is thus a good candidate to approximate the S-D case. The
distances corresponding to the two regimes have been fixed to L0 = 10b and
L1 = 20b. The other numerical parameters are the same as in the S-D simu-
lation, i.e. the mesh size h/b = 1/100, the time step ∆t/t0 = 6.45 10−5 and
the regularization parameter σ = 2h. The simulation required ∼ 65M of grid
points and was carried out using 40 processors for 4, 800 time steps (it ran 170
hours).
Comparison and analysis of the space-developing and of

the time-developing simulations
Fig. 5.15 shows a general 3D view of the S-D vortex sheet rollup at Re =
Γ0 /ν = 104 . Information about the circulation profile is an important issue in
order to model trailing vortices. Many analytical circulation profile formula,
also called “vortex models”, and their related velocity profile, have been devel-
oped. The most usual models are itemized in details in [25] (see also Table 4.1
for a summary). We also connect the 2-D rollup simulations and the 3-D S-D
simulations using:
x = U0 t = U0 t0 τ,
with τ = t/t0 the dimensionless time for the 2-D simulation, with t0 = b0 /V0
(where V0 = Γ0 /(2πb0 ) is the wake descent velocity. A station at x/b int the
3-D S-D simulation thus has an equivalent dimensionless time equal to
x/b
τ= 2

U0 2πb0
b Γ0
Figure 5.15: 3-D visualization of a S-D vortex sheet rollup. Isocontour of the vor-
ticity field norm ||ω||b2 /Γ0 = 8.0 colored by the axial vorticity component. A set of
streamlines are also plotted in the wingtip region, going above/under the wing (or-
ange/cyan lines).
Characterization of the wake Fig. 5.17 presents the evolution of the vortex
core radius rc for the S-D case and the T-D case. During the formation of the
vortex pair, the vortex core quickly increases until x/b ≃ 5 (τ = 0.16). As
foreseen, the S-D simulation and the T-D simulation are very similar in terms
of rc and uθ,max.
Fig. 5.16 shows the circulation profiles Γ(r) and the tangential velocity
profiles uθ (r) = Γ(r)/(2πr) that we have obtained for various cross planes.
The maximum tangential velocity slightly decreases by diffusion from (8.2V0 to
7.0V0 ). Fig. 5.18 and Fig. 5.19 show the axial velocity. We analyze the value
of the maximum axial velocity (positive and negative) in Fig. 5.18 and provide
cross-views in Fig. 5.19. It is interesting to note that both are of the same order
of magnitude. We find a peak of axial velocity deficit (i.e. towards the wing)
in the vortex core roughly valued to 0.5V0 (when the rollup is completed and
increasing shortly after). Moreover, an axial velocity in the opposite direction
(away from the wing) surrounds the vortex core at the peak value of roughly
0.4V0 .
(a)
1
0.9
0.8
0.7
0.6
Γ(r/b)/Γ0
x/b = 4.0
x/b = 5.0
0.5
x/b = 6.0
x/b = 7.0
0.4
x/b = 8.0
0.3 x/b = 9.0
x/b = 10
0.2 x/b = 11
x/b = 12
0.1 x/b = 13
x/b = 14
0
0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5
r/b
(b)
9
x/b = 4.0
8 x/b = 5.0
x/b = 6.0
7 x/b = 7.0
x/b = 8.0
6 x/b = 9.0
x/b = 10
uθ (r/b)/V0
5 x/b = 11
x/b = 12
x/b = 13
4
x/b = 14
3
0
0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5
r/b
Figure 5.16: Circulation profiles (a) and tangential velocity profiles (b) at various
distances in the wake.
(a) (b)
8.5
0.06 8
7.5
uθ,max /V0
0.055
7
rc /b
0.05 6.5
6
0.045
5.5
0.04 5
0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5
τ τ
Figure 5.17: Evolution of the vortex effective core size rc (a) and of the maximum
tangential velocity uθ,max (b) for the T-D simulation (tick line) and for the S-D
simulation (thin line).
1.5
0.5
ux /V0
−0.5
−1
−1.5
0 2 4 6 8 10 12 14
x/b
Figure 5.18: Evolution of the maximum axial velocity (solid) and the minimum axial
velocity (dash) with respect to the distance in the wake.
(a) x/b = 0.0 (τ = 0.0)
(b)
Figure 5.19: (a) Slice of the axial component of the vorticity field ωx t0 . (b) Slice of
the axial component of the velocity field ux /V0 .
(a) x/b = 4.0 (τ = 0.13)
(b)

(a) x/b = 8.0 (τ = 0.26)
(b)

(a) x/b = 12.0 (τ = 0.39)
(b)

Global dynamics Fig. 5.20 presents the 3-D view of the S-D evolution of the
accelerated wing and the resulting end-effect. While the wing accelerates, the
vorticity generated by the acceleration rolls up arround the two wake vortices.
This creates an external deformation of the vortices. Afterwards, a helical
instability develops that moves towards the stronger vortex and irremediably
deforms the wake vortices. Fig. 5.21 shows the 3-D view of the end-effect for
the T-D simulation. The same physics occur.
Visualization of the end-effect Fig. 5.22 shows a zoom, in the propagation

region, of the vorticity isocontour at a low level of vorticity norm. In the
S-D case, the end-effect appears to be more confined. The reason could be
the fact that, in this case, an axial velocity (in the êx direction) surrounds
the vortex and then convects, in the wake, the external deformation which
propagates in the opposite direction. However, the pressure wave is not affected
by the slowing down of the external helical instability. Fig. 5.23 shows two
isocontours of vorticity norm for higher levels in order to visualize the inner
vortices. A pressure wave propagates along each vortex which creates a strong
axial velocity, destabilizes each vortex and allows the development of a helical
instability.
Propagation of the instability Fig. 5.24 compares the evolution of the

axial velocity front (uw ). This velocity is slightly higher in the S-D case :
uw /uθ,max = 1.30 instead of 1.20 for the T-D case. The difference might
partially come from the additional axial velocity in the vortex core, valued at
≃ 0.5V0 . However, the global structure of the axial velocity is quite similar
between the two cases, as shown in Fig. 5.25. The simulations were carried
out at relatively high Reynolds number; the vortex core diffusion doesn’t play
much of a role and the front location evolves quasi linearly in time.
As done in the first section, we can analyze the axial component of the
vorticity perturbation field ωx′ . Fig. 5.26 compares the isocontours of this field
between the T-D simulation and the S-D simulation. The vorticity base fields
(0)
ω x are here taken as the vorticity field at the center slice of the high speed
region for the T-D case, and as the vorticity field at the same time-evolving slice
but for a rollup simulation at high speed regime for the S-D case. It confirms
the conclusion of Moet[59] : an axisymmetric instability is at first created

which magnify the vortex core. This instability is followed by a helical one (the
positive and the negative levels of vorticity perturbations are interlaced) that
leads to the loss of the vortex core coherence.

Γ1
tuθ,max/b = 0.0 t 2πb 2 = 0.0
0

Γ1
tuθ,max/b = 0.42 t 2πb 2 = 0.065
0
Figure 5.20: 3-D view of the wake for the space-developing simulation of an ac-
celerated wing. Isocontour of the vorticity norm ||ω|| b2 /Γ1 = 8 colored by the axial
vorticity component.

Γ1
tuθ,max/b = 1.2 t 2πb 2 = 0.19
0

Γ1
tuθ,max/b = 2.1 t 2πb 2 = 0.32
0


Γ1
tuθ,max/b = 2.9 t 2πb 2 = 0.45
0

Γ1
tuθ,max/b = 3.3 t 2πb 2 = 0.52
0


Γ1
tuθ,max/b = 0.0 t 2πb 2 = 0.0
0

Γ1
tuθ,max/b = 0.33 t 2πb 2 = 0.052
0
Figure 5.21: 3-D view of the wake for the time-developing simulation of an acceler-
ated wing. Zoom of the vorticity norm isocontour ||ω|| b2 /Γ1 = 8 colored by the axial
vorticity component.

Γ1
tuθ,max/b = 1.0 t 2πb 2 = 0.15
0

Γ1
tuθ,max/b = 1.7 t 2πb 2 = 0.26
0


Γ1
tuθ,max/b = 2.7 t 2πb 2 = 0.41
0

Γ1
tuθ,max/b = 4.0 t 2πb 2 = 0.62
0

tuθ,max /b = 0.33
S-D :
T-D :
êx êz
tuθ,max /b = 1.2
?
S-D : êy
T-D :
tuθ,max /b = 2.0
S-D :
T-D :
tuθ,max /b = 2.8
S-D :
T-D :
Figure 5.22: Top view of the vorticity field isocontour ||ω||b2 /Γ1 = 8 colored by the
axial vorticity component. Comparison of the instability evolution between the space-
developing (S-D) and the time-developing (T-D) simulations for various times. The
length of the visualization box is 4b centered relatively to xref .
Space-developing Time-developing
xref /b − 3.0
xref /b − 2.5
êz -
êy
?
êx
xref /b − 1.0
xref /b
Figure 5.23: Top view of the vorticity isocontours colored by the axial vorticity com-
ponent ||ω||b2 /Γ1 = 40 (low opacity) and ||ω||b2 /Γ1 = 60 (high opacity). Compari-
son between the space-developing simulation (left) and the time-developing simulation
(right) at the time tuθ,max /b = 2.8. Also shown are the slices used Fig. 5.26.
3.5
3
(xref − x)/b
2.5
1.5 uif /uθ,max(0) ∼ 1.20

uif /uθ,max(0) ∼ 1.30
1
0.5
0.8 1 1.2 1.4 1.6 1.8 2 2.2 2.4 2.6
Γ1
tuθ,max/b = 6.4 2πb 2
0
Figure 5.24: Evolution of the non-zero axial velocity front for the space-developing
simulation (dash) and the time-developing simulation (solid). The thin lines corre-
sponds to the least-squared fitted lines.
(a)
z/b
(xref − x)/b
(b)
z/b
(xref − x)/b
Figure 5.25: Cross plane view of the axial velocity field at y/b = 0.5π/4 for (a)
the space-developing simulation and (b) the time-developing simulation at the time
tuθ,max /b = 2.8.
Space-developing Time-developing
(a) 0 −0.1
−0.1 −0.2
−0.2 −0.3
z/b
z/b
−0.3 −0.4
−0.4 −0.5
−0.5 −0.6
0.1 0.2 0.3 0.4 0.5 0.6 0.1 0.2 0.3 0.4 0.5 0.6
y/b y/b
(b) 0 −0.1
−0.1 −0.2
−0.2 −0.3
z/b
z/b
−0.3 −0.4
−0.4 −0.5
−0.5 −0.6
0.2 0.3 0.4 0.5 0.6 0.7 0.1 0.2 0.3 0.4 0.5 0.6
y/b y/b
0
(c)
−0.1
−0.1
−0.2
−0.2
−0.3
z/b
z/b
−0.3
−0.4
−0.4
−0.5
−0.5
−0.6
−0.6
0.1 0.2 0.3 0.4 0.5 0.6 0.2 0.3 0.4 0.5 0.6 0.7
y/b y/b
Figure 5.26: Ten equally spaced levels of the axial vorticity isocontour ωx′ b2 /Γ1 from
−5.5 to 5.5, at tuθ,max /b = 2.8, for various cross planes : (a) x/b = xref /b − 3, (b)
x/b = xref /b − 2.5 and (c) x/b = xref /b − 1. Positive levels are plotted in red and
the negative levels are plotted in blue. The circle with r = rc is also plotted in black.
5.5. Conclusion 159
5.5 Conclusion
The VIC-PFM method was shown to be suitable to simulate the dynamics
of end-effects phenomena, both for the time-developing (T-D) and the space-
developing (S-D) simulations. We observed comparable results to those of
CNRS/IRPHE. The propagation speeds are found to be in good agreement
with the theoretically predicted values. The same dynamics were observed for
each simulated case, confirming the observations of Moet et al.[59] : at first
an axisymmetric perturbation, corresponding to a pressure wave, propagates
along the vortex core; if the resulting axial velocity is strong enough, as in
the configurations investigated here, this wave can destabilize the vortex and
a helical instability then develops.
The S-D simulation with wing acceleration points out a difference compared
to the time-developing simulation in terms of propagation speed. The differ-
ence however remains minor and does not modify the dynamics.
Note however that the present S-D simulation, and the resulting 3-D rollup,
underestimates the intensity of the axial velocity within the vortex core. It
would likely be much higher considering a 3-D rollup with momentum deficit
just behind the wing corresponding to the boundary layers and profile drag.
This was not done in the present work; yet it constitutes an interesting follow
up. It also can explain the different behaviour of the “end-effect” associated
with an acceleration or a deceleration of the wing, as observed experimentally.
Chapter 6
Conclusions & perspectives
6.1 Achievements
This thesis has been devoted to the simulation of vortical flows using the vor-
tex particle-mesh (VPM) method, also sometimes called “vortex-in-cell” (VIC)
method. The contributions of this work are twofold: they consist of method-
ological improvements and of investigations of complex vortical flows.
Our work in the vortex particle-mesh method itself has brought two out-
standing advances:
• In the spirit of the immersed boundary (IB) method, an approach has

been successfully developed for the simulation of external flows past
arbitrary-shape bodies (i.e. bluff- and profiled-bodies). The approach
is based on the methodology proposed by Ploumhans et al. [68] that
was then developed for the“classical” vortex particle method (VP, with-
out grid). It was adapted and extended to the VPM framework where
the Poisson equation is solved on a grid; here using the FISHPACK li-
brary. The high fidelity simulation of external flows benefits from the
reduction of the computational cost associated to this combined ap-
proach. The developments were done in the framework developed by Co-
cle et al. [11, 13] where the VPM method was combined with the parallel
fast multipole method (PFM) for efficient domain decomposition simula-
tions. The present VPM-PFM-IB was successfully validated against the
162 Chapter 6. Conclusions & perspectives
well-referenced test case of a flow past the sphere at Re = 300. Further-

more, the method was used to capture the flow past a low aspect-ratio
wing at Re = 3000.
• A mesh refinement handling within this VPM method has also been de-
veloped; including the IB feature. The algorithm implements a wavelet-
based conservative interpolation scheme of vortex particles patches with
various resolutions. It also permits a consistent reinitialization of the vor-
tex particle strengths. Good conservation properties were observed. The
method has also been validated on the flow past a sphere at Re = 300.
Furthermore, this multi-resolution capability has also been used, together
with multiscale subgrid modeling, to perform a large eddy simulation of a
wake-vortex flow at high Reynolds number (Re = 106 ). It is worth noting
that this simulation would not be affordable using a uniform resolution.
This contribution thus definitely increases the possibilities of turbulent
vortical flow simulations, for which we typically observe multiple scales
within the flow.
In the context of air traffic management, a preponderant factor is the aircraft

induced wake-vortex hazard. It is at the origin of the separation standard for
landing and take-off phases. Our work on the investigation of such vortical
flows has also brought outstanding advances:
• We have investigated the influence of a velocity deficit (due to the friction

drag) on the wake rollup. For that purpose, two space-developing (S-D)
simulations, with and without velocity deficit effect, have been carried
out. They have been compared. The former becomes quickly turbulent
while the latter remains stationary. In addition to those simulations,
a complementary time-developing (T-D) simulation with velocity deficit
has been performed. We found some differences in the instability develop-
ment. In opposition to the S-D simulation, a helical instability develops
around the primary vortex-cores. Finally, the net-result of the rolled-
up vortex-wake was characterized and the most common vortex-models
were fitted: it was confirmed that two-length scale models are required to
properly fit both the vortex-core size and the maximal tangential velocity.
6.2. Perspectives 163
• We have studied the end-effect phenomenon, associated to an axial varia-

tion of the vortex/wake circulation. Fundamental differences were found
between the results obtained by the simulations and those obtained by
simplified theory. The application of the VPM method to these flows has
shown that the VPM method does constitute the method of choice for
such problems.
6.2 Perspectives
Finally, the research work reported in this thesis raises a few questions and
hints at potential future developments. We close this thesis with a discussion
of these open issues.
6.2.1 A VPM method using a staggered grid
VPM methods are based on the resolution of the Poisson equation using an
efficient grid solver such as FISHPACK: i) to obtain the streamfunction from the
vorticity field and ii) to reproject the vorticity field in a solenoidal space. This
provides non-fully consistent streamfunction: i.e., we do not exactly recover the
vorticity field by taking twice the grid-evaluated curl of the obtained stream-
function. ∇h × (∇h × Ψ) is not exactly equal to ω. However, it turns out that
for usual applications, the results are accurate enough if the vorticity diver-
gence is maintained at an acceptable, and low, level. The same conclusion is
drawn concerning the correction procedure of the vorticity field. Nevertheless,
the initial level of vorticity divergence can be quite high for wall-bounded flows
using immersed-boundary techniques. Moreover, this high level of divergence
is convected in the wake. This can lead to local and/or global inaccuracies. An
elegant way to address this issue would be to solve the Poisson equation using a
staggered grid, in the same spirit as what is done in classical velocity-pressure
based Navier-Stokes solvers, for which a Poisson equation is solved to make the
velocity field divergence-free.
This approach was recently also proposed and tested by Uchiyama and
Yoshii [84]. They definitely obtain good results in terms of self-consistency and
accuracy.
6.2.2 Further improvements of the Poisson solver for hi-

erarchically refined grids
The presented approach to solve the Poisson equations is not necessarily opti-
mal in terms of computational cost and parallel efficiency. We here combined a
fast direct solver (e.g. FISHPACK) with the parallel fast multipole (PFM) solver;
in order to obtain the required boundary conditions. This approach is very ver-
satile and robust but it requires that each processor be highly loaded. It results
in a poor strong-scaling and an acceptable weak-scaling (Cocle et al. [13]). This
behavior is explained by the following reason: the computational cost associ-
ated to the PFM evaluations is relatively high; given a fixed problem size, when
increasing the number of processors, the PFM code is used more and more and
that finally ends up with a computational time comparable to that of the classi-
cal vortex method using PFM only (i.e. no grid). This is even more true when
using multi-resolution because many grid blocks that have few grid points are
created leading to a less efficient load-balancing.
A good strategy for improvement would be the use of a multigrid solver for
hierarchically refined grids; as first suggested by Brandt [5] and then improved
by Brown and Lowe [7] and Meglicki et al. [56]. The Poisson equation would
be solved on a bigger grid, requiring then fewer PFM evaluations.
An interesting project is the open development platform called “CHOMBO”
which handles multi-physics and AMR abilities as well as a multigrid solver
for Poisson equation (see Martin et al. [55]). An equivalent platform is the
Overture framework (see Brown et al. [6]). We suggest to develop a VPM
method using these platforms. This could be fairly easily implemented, lead-
ing to an even more efficient numerical tool for vortical flows. It is worth
noting that another interesting method exists for solving a Poisson equation
on hierarchically refined grids, and which is still a direct solver (see Huang and
Greengard [38] and Ricker [71]).
The overall algorithm can also be further improved in terms of accuracy.
So far, on a common surface between two blocks with different resolution, the
Laplacian operator is computed independently from each side. Meaning that
the diffusion flux is likely not identical across the surface. This potentially leads
to some vorticity source in the flow. The Laplacian scheme can be modified to
avoid this issue (as done by Martin et al. [55]). We note that this feature is
already implemented in the AMR frameworks cited above.
6.2.3 Improvement of the handling of no-slip

The approach to capture the flow past a no-slip boundary is summarized as
follows: i) the panel sources that cancel the residual slip velocity are computed
using the Boundary Element Method (BEM), ii) those panel sources are re-
distributed on the grid (and thus smoothed), in order to obtain the velocity
by solving the Poisson equation ∇2 Ψ = −ω and iii) adding the panel sources
into the flow through an equivalent diffusion flux. The major drawback of this
approach is a certain inconsistency between, on the one hand, the velocity field
used to compute the panel sources and, on the other hand, the velocity field
obtained after solving the Poisson equation. This leads to a time oscillation of
the vorticity flux at the wall and thus of the force acting on the body.
The consistency of the method can be improved by means of the evaluation
of the boundary sources that are already redistributed on the grid (as done
by Cottet and Poncet [17] and Poncet [70]). This however leads to a limited
accuracy since the discontinuous vorticity solution is still not captured. A
better approach is likely the use of a sharp-interface method (see Udaykumar
et al. [85] and Marella et al. [54] and Yang and Udaykumar [92]) for which the
near-wall operators are modified in order to handle boundary-discontinuities,
preserving the discretization accuracy up to the wall. This approach can be
enhanced using the immersed interface method (see Leveque and Li [50] and Li
and Ito [51] for more details on the method and Morgenthal and Walther [61]
for a 2D VPM application).
6.2.4 A coupling approach using an Eulerian “near-wall”

solver
The combined VPM-PFM-IB method is costly for wall-bounded external flows
at high Reynolds number. As the fundamental element of discretization is
still isotropic (h × h × h), the method is intrinsically unadapted to capture
thin boundary layers which are a main characteristic of high Reynolds number
flows. A significant gain in terms of computational cost can be achieved by
means of a Eulerian-Langragian coupling (see Daeninck and Winckelmans [23]

and Waerenburgh [86]). The coupling approach, proposed by Daeninck [23], has
been validated in 2-D using a “classical” vortex particle method (i.e., no grid),
together with a finite difference code for incompressible flows. We here suggest
to couple, on the one hand, the efficient combined VPM-PFM-IB method and,
on the other hand, a hybrid finite volume/finite element code for compressible
flows (i.e. the ARGO code developed by CENAERO). This approach aims to
make the two codes agree, at least for a highly subsonic regime (i.e. Mach
number < 0.1 · · · 0.2).
The coupling approach is summarized as follows. The Lagrangian code has

a global view of the flow; it contains the entire vorticity field. It is responsible
for capturing the flow from the wall to the wake. The boundary layers are, how-
ever, much under-resolved, yet their particle strengths are correct (see below).
The Eulerian code, allowing a body-fitted mesh, only captures the “near-wall”
region. The Eulerian vorticity field is used to correct the near-wall Lagrangian
vorticity field hence, the particle strengths in the near wall region are correct
(as this corresponds to a fine to coarse procedure). The required boundary
conditions on the outer part (i.e. velocity, pressure and temperature) are pro-
vided by the Lagrangian code. So far, only Dirichlet type boundary conditions
were investigated. The pressure field is computed using a supplementary scalar
Poisson equation (see Ploumhans [66]), whereas the temperature field is mod-
eled.
• Pressure-field evaluation
A Poisson equation (and thus also an integral equation) is solved in order

to compute the pressure field from the vorticity and the velocity fields :
∇2 B = −∇ · (ω × u),
where B = (p/ρ∞ + ½u · u) − (p∞ /ρ∞ + ½U∞ · U∞ ). The boundary con-

∂B
dition to impose at the body surface is ∂n = −ν(∇×ω)·n. This equation
is solved using the boundary element method :

Z Z
∂G ′ ′ ′
ξB(x) − (|x − x |)B(u ) dx = ν G(|x − x′ |)(∇′x × ω(x′ )) · n(x′ )dx′
S ∂n′
ZS
+ G(|x − x′ |)∇′x · (ω(x′ ) × u(x′ ))dx′ .
V
where ξ = ½ for a point x located on the surface body and ξ = 1 otherwise.

The integral equation and the Poisson equation are solved in the same
way as those involved in the immersed boundary technique.
• Temperature-field estimation
The VPM-PFM-IB solves the Navier-Stokes equation for incompressible

flows and thus it does not provide the temperature field. Various hy-
potheses can be proposed to recover an image of the temperature field
knowing the upstream conditions (p∞ , ρ∞ , T∞ , U∞ ) :
1. Assuming the temperature is such that the density is conserved:

T (x) = p(x)/(ρ∞ R), where R is the gas constant.
2. Assuming an isothermal flow T (x) = T∞ .

γ
3. Assuming a homentropic transformation: T (x) = T∞ (p∞ /p(x)) γ−1
cp
where γ = cv is the ratio of the specific heats.
m
4. Assuming a polytropic transformation T (x) = T∞ (p/p∞ ) m−1 , where
1 ≤ m ≤ γ.
The first three hypotheses have been compared for the Euler flow past a
sphere and at a Mach number M = √ U∞ = 0.1. The ARGO code was
γRT∞
used. The mesh is such that the external boundary is close to the body.
It essentially consists in a sphere inside a cube, where the side c = 32 D (D
being the sphere diameter) using ∼ 220 000 tetrahedrons, see Fig. 6.1.(a).
“Incompressible” boundary conditions (described before) are imposed. It
turns out that the hypotheses 2 and 3 provide satisfactory results (and
thus also likely hypothesis 4, not tested); the pressure field corresponding
to the case 3 is presented in Fig. 6.1.(b).
The complete coupling approach (VPM-PFM-IB-ARGO), that is with vis-

cosity effects, has been implemented but not yet validated. The impulsively
(a) (b)
Figure 6.1: (a) View of the mesh of a sphere included in a box. (b) Slice of the
pressure field (p − p∞ )/p∞ ∗ 100 just after impulsive start assuming a isothermal flow.
started flow past a sphere at Re = 300 and M = 0.1 has been computed. The
mesh used for each code is presented in Fig. 6.2.(a). Fig. 6.2.(b) shows a slice of
the velocity norm for both domains at the early age tU∞ /D = 1.5. We clearly
see that velocity field is continuous through each domain.
These are preliminary results which seem to indicate that the methodology
is coherent and could lead to a very efficient tool. Much work yet remains
in order to “get there”. Wind farms computations constitute a good example
of problems for which it is crucial to capture aerodynamics coupled to wake
dynamics and for which this methodology should be investigated.
(a)
(b)
Figure 6.2: (a) Quarter of the mesh used by the compressible ARGO code (top) and
by the incompressible VIC-PFM-IB code (bottom). (b) 2-D cut of the velocity norm
|u|/U∞ at t UD∞ = 1.5. Eulerian field inside the circle and Lagrangian field outside
the circle.
Bibliography
[1] J. Barnes and P. Hut. A hierarchical O(N Log N) force-calculation algo-

rithm. Nature, 324:446–449, 1986.
[2] M. Bergdorf, G.-H. Cottet, and P. Koumoutsakos. Multilevel adaptive

particle methods for convection-diffusion equations. Multiscale Modeling
and Simulation: A SIAM Interdisciplinary Journal, 4(1):28–357, 2005.
[3] M. Bergdorf and P. Koumoutsakos. A lagrangian particle-wavelet method.

Multiscale Modeling and Simulation: A SIAM Interdisciplinary Journal,
5(3):980–995, 2006.
[4] M. J. Berger and J. Oliger. Adaptive mesh refinement for hyperbolic

partial differential equations. J. Comput. Phys., 53(3):484–512, 1984.
[5] A. Brandt. Multi-level adaptive solutions to boundary-value problems.

Math. Comp., 31:333–390, 1977.
[6] D. Brown, W. Henshaw, and D. Quinlan. Overture: An object-oriented

framework for solving partial differential equations. In Yutaka Ishikawa,
Rodney Oldehoeft, John Reynders, and Marydell Tholburn, editors, Sci-
entific Computing in Object-Oriented Parallel Environments, volume 1343
of Lecture Notes in Computer Science, pages 177–184. Springer Berlin /
Heidelberg, 1997.
[7] D. Brown and L. Lowe. Multigrid elliptic equation solver with adaptive
mesh refinement. J. Comput. Phys., 209(2):582–598, 2005.
[8] P. Chatelain, A. Curioni, M. Bergdorf, D. Rossinelli, W. Andreoni, and

P. Koumoutsakos. Billion vortex particle direct numerical simulations of
172 BIBLIOGRAPHY
aircraft wakes. Computer Methods in Applied Mechanics and Engineering,

197(13-16):1296–1304, 2008.
[9] P. Chatelain and P. Koumoutsakos. A fourier-based elliptic solver for

vortical flows with periodic and unbounded directions. J. Comput. Phys.,
229(7):2425–2431, 4 2010.
[10] J. P. Christiansen. Numerical solution of hydrodynamics by the method

of point vortices. J. Comput. Phys., 13:363–379, 1973.
[11] R. Cocle. Combining the vortex-in-cell and parallel fast multipole meth-
ods for efficient domain decomposition simulations : DNS and LES ap-
proaches. PhD thesis, Université catholique de Louvain, Louvain-la-Neuve,
2007.
[12] R. Cocle, L. Dufresne, and G. Winckelmans. Investigation of multiscale

subgrid models for les of instabilities and turbulence in wake vortex sys-
tems. In Lecture Notes in Computational Science and Engineering series,
Symposium on Complex Effects in Large-Eddy Simulation (CY-LES 2005),
Limassol (Cyprus), January 2007. Springer, January 2007.
[13] R. Cocle, G. Winckelmans, and G. Daeninck. Combining the vortex-in-

cell and parallel fast multipole methods for efficient domain decomposition
simulations. J. Comput. Phys., 227:2263–2292, 2008.
[14] G. Constantinescu and K. Squires. Numerical investigations of flow

over a sphere in the subcritical and supercritical regimes. Phys. Fluids,
16(5):1449–1466, 2004.
[15] M. Coquerelle and G.-H. Cottet. A vortex level set method for the two-way
coupling of an incompressible fluid with colliding rigid bodies. Journal of
Computational Physics, 227(21):9121 – 9137, 2008. Special Issue Celebrat-
ing Tony Leonard’s 70th Birthday.
[16] G.-H. Cottet and P. Koumoutsakos. Vortex Methods: Theory and Practice.
Cambridge University Press, first edition, 2000.
[17] G.-H. Cottet and P. Poncet. Advances in direct numerical simulations

of 3d wall-bounded flows by vortex-in-cell methods. J. Comput. Phys.,
193(1):136–158, 2004.
BIBLIOGRAPHY 173
[18] C. Cottin. Experimental investigation of the dynamics of counter-rotating

vortex pairs in ground effect. Master’s thesis, Chalmers University of
Technology, 2006.
[19] C. Cottin, G. Winckelmans, T. Lonfils, and R. Cocle. Roll-up of a

temporally-evolving wing wake with velocity deficit. Far-wake research
project technical report 2.2.2-3, Université catholique de Louvain (UCL),
Louvain-la-Neuve, 2007.
[20] S. C. Crow. Stability theory for a pair of trailing vortices. AIAA Journal,
8(12):2172–2178, 1970.
[21] G. Daeninck. Developments in hybrid approaches. Vortex method with

known separation location. Vortex method with near-wall Eulerian solver.
RANS-LES coupling. PhD thesis, Université catholique de Louvain, 2006.
[22] G. Daeninck, R. Cocle, and G. Winckelmans. FAR-Wake, TR-1.1.2-1-les

of spatialy evolving wakes in ground effect. FAR-Wake technical report
TR3.1.2-4 (42 pages), Université catholique de Louvain, 2006.
[23] G. Daeninck and G. Winckelmans. A new eulerian-lagrangian vortex

method for flows with massive separation. In 6th National congress on
Theoretical and Applied Mechanics, Ghent, 1993.
[24] A.C. de Bruin. Estimation of exhaust velocities and temperatures for

various operation phases of a modern high by-pass turbofan. FAR-Wake
research project delivrable d 2-1.1-1b, 2005.
[25] A.C. de Bruin and G. Winckelmans. Cross-flow kinetic energy and core
size growth of analytically defined wake vortex pairs. Technical report,
Nationaal Lucht- en Ruimtevaartlaboratorium (NLR), Aug. 2005.
[26] P. de Vos and M. Janssens. Investigation expérimentale de tourbillons de

sillage d’aile en phase d’atterrissage. Master’s thesis, Université catholique
de Louvain (UCL), 2007.
[27] O. Desenfans and S. Gigantelli. Instabilités des tourbillons de sillage d’aile

avec décélétation brusque : investigation dans le bassin de halage. Master’s
thesis, Université catholique de Louvain (UCL), 2004.
174 BIBLIOGRAPHY
[28] J. Eldredge, C. Colonius, and A. Leonard. A vortex particle method for

two-dimensional compressible flow. J. Comput. Phys., 179(2):371–399,
2002.
[29] D. Fabre, Sipp D., and L. Jacquin. The kelvin waves and the singular
modes of the lamb-oseen vortex. J. Fluid Mech., 551:235–274, 2006.
[30] D. Fabre and L. Jacquin. Viscous instabilities in trailing vortices at large

swirl numbers. J. Fluid Mech., 500:239–262, 2003.
[31] L. Georges. Development and validation of a LES methodology for com-

plex wall-bounded flows : application to high-order structured and indus-
trial unstructured solvers. PhD thesis, Université catholique de Louvain,
Louvain-la-Neuve, 2007.
[32] C. Geuzaine and J.-F. Remacle. Gmsh: a three-dimensional finite element

mesh generator with build-in pre- and post-processing facilities. Interna-
tional Journal for Numerical Methods in Engineering, 79(11):1309–1331,
2009.
[33] L. Greengard. The rapid evaluation of potential fields in particle systems.

PhD thesis, 1987.
[34] L. Greengard and V. Rokhlin. A fast algorithm for particle simulations.

J. Comput. Phys., 73(2):325–348, 1987.
[35] B. Hejazialhosseini, D. Rossinelli, M. Bergdorf, and P. Koumoutsakos.

High order finite volume methods on wavelet-adapted grids with local
time-stepping on multicore architectures for the simulation of shock-
bubble interactions. Journal of Computational Physics, 229(22):8364 –
8383, 2010.
[36] S. E. Hieber and P. Koumoutsakos. An immersed boundary method for

smoothed particle hydrodynamics of self-propelled swimmers. J. Comput.
Phys., 227:8636–8654, October 2008.
[37] R. W. Hockney. A fast direct solution of poisson’s equation using fourier

analysis. J. ACM, 12(1):95–113, 1965.
BIBLIOGRAPHY 175
[38] J. Huang and L. Greengard. A fast direct solver for elliptic partial dif-
ferential equations on adaptively refined meshes. SIAM J. Sci. Comput.,
21(4):1551–1566, 1999.
[39] L. Jacquin, D. Fabre, D. Sipp, and E. Coustols. Unsteadiness, instability

and turbulence in trailing vortices. Comptes Rendus Physique, 6:399, 2005.
[40] L. Jacquin, D. Fabre, D. Sipp, V. Theofilis, and H/ Vollmers. Instability

and unsteadiness of aircraft wake vortices. Aerosp. Sci. Tech., 7:577–593,
2003.
[41] J. Jeong and F. Hussain. On the identification of a vortex. J. Fluid Mech.,

285(-1):69–94, 1995.
[42] T. A. Johnson and V. C. Patel. Flow past a sphere up to a Reynolds

number of 300. J. Fluid Mech., 378:19–70, 1999.
[43] G. Karypis and V. Kumar. METIS, a Software Package for Partition-

ing Unstructured Graphs and Computing Fill-Reduced Orderings of Sparse
Matrices.
[44] G. Karypis and V. Kumar. A fast and high quality multilevel scheme for
partitioning irregular graphs. SIAM J. Sci. Comput., 20(1):359–392, 1998.
[45] P. Koumoutsakos. Direct numerical simulations of unsteady separated

flows using vortex methods. PhD thesis, California Institute of Technology,
1993.
[46] P. Koumoutsakos. Multiscale flow simulations using particles. Annual

Review of Fluid Mechanics, 37(1):457–487, 2005.
[47] P. Koumoutsakos and A. Leonard. High-resolution simulations of the flow

around an impulsively started cylinder using vortex methods. J. Fluid
Mech., 296(1):1–38, 1995.
[48] S. Leibovich. Wave propagaton, instability, and breakdown of vortices.

In H. G. Hornung and Muller E.-A., editors, Vortex Motion, pages 50–67,
1982.
176 BIBLIOGRAPHY
[49] S. Leibovich. Vortex breakdown: A coherent transition trigger in con-

centrated vortices. In M. Lesieur and O. Metais, editors, Turbulence and
coherent structures, 1990.
[50] R. J. Leveque and Z. Li. The immersed interface method for elliptic equa-
tions with discontinuous coefficients and singular sources. SIAM Journal
on Numerical Analysis, 31(4):1019–1044, 1994.
[51] Z. Li and K. Ito. The Immersed Interface Method: Numerical Solutions

of PDEs Involving Interfaces and Irregular Domains (Frontiers in Applied
Mathematics). Society for Industrial and Applied Mathematics, Philadel-
phia, PA, USA, 2006.
[52] H. Liu, S. Krishnan, S. Marella, and H.S. Udaykumar. Sharp interface

Cartesian grid method II: A technique for simulating droplet interactions
with surfaces of arbitrary shape. J. Comput. Phys., 210(1):32–54, 2005.
[53] O. Lucca-Negro and T. O’Doherty. Vortex breakdown: a review. Progress

in Energy and Combustion Science, 27(4):431 – 481, 2001.
[54] S. Marella, S. Krishnan, H. Liu, and H.S. Udaykumar. Sharp interface

Cartesian grid method I: An easily implemented technique for 3D moving
boundary computations. J. Comput. Phys., 210(1):1 – 31, 2005.
[55] D. Martin, P. Colella, and D. Graves. A cell-centered adaptive projection

method for the incompressible navier-stokes equations in three dimensions.
J. Comput. Phys., 227(3):1863–1886, 2008.
[56] Z. Meglicki, S. Gray, and B. Norris. Multigrid FDTD with chombo. Com-
puter Physics Communications, 176(2):109–120, 2007.
[57] P. Meunier. FAR-Wake, TR-1.1.2-1 axial propagation of vortex perturba-

tions. Technical report, CNRS-IRPHE, 2006.
[58] R. Mittal and G. Iaccarino. Immersed boundary methods. Ann. Rev. Fluid
Mech., 37:239–261, 2005.
[59] H. Moet, F. Laporte, G. Chevalier, and T. Poinsot. Wave propagation in

vortices and vortex bursting. Phys. Fluids, 17(5):054109 (15 pages), 2005.
BIBLIOGRAPHY 177
[60] J.J Monaghan. Extrapolating b splines for interpolation. J. Comput.

Phys., 60(2):253–262, 1985.
[61] G. Morgenthal and J.H. Walther. An immersed interface method for the
vortex-in-cell algorithm. Comput. Struct., 85(11-14):712–726, 2007.
[62] F. Noca, D. Shiels, and D. Jeon. Measuring instantaneous fluid dynamic

forces on bodies, using only velocity fields and their derivatives. Journal
of Fluids and Structures, 11(3):345–350, 1997.
[63] F. Noca, D. Shiels, and D. Jeon. A comparison of methods for evaluating

time-dependent fluid dynamic forces on bodies, using only velocity fields
and their derivatives. Journal of Fluids and Structures, 13(5):551–578,
1999.
[64] M. Ould Salihi. Couplage de méthodes numériques en simulation directe

d’écoulements incompressibles. PhD thesis, Université Joseph Fourier,
Grenoble, 1998.
[65] C.S. Peskin. The immersed boundary method. Acta Numerica, pages
479–517, 2002.
[66] P. Ploumhans. Simulation of High Reynolds Number Flows Past Bluff Bod-
ies Using Vortex and Boundary Element Methods. PhD thesis, Université
catholique de Louvain, Louvain-la-Neuve, 2001.
[67] P. Ploumhans and G. Winckelmans. Vortex methods for high-resolution

simulations of viscous flow past bluff bodies of general geometry. J. Com-
put. Phys., 165:354–406, 2000.
[68] P. Ploumhans, G. Winckelmans, J. Salmon, A. Leonard, and M. Warren.

Vortex methods for direct numerical simulation of three-dimensional bluff
body flows: application to the sphere at Re = 300, 500, and 1000. J.
Comput. Phys., 178(2):427–463, 2002.
[69] P. Poncet. Vanishing of mode b in the wake behind a rotationally oscil-

lating circular cylinder. Phys. Fluids, 14(6):2021–2023, 2002.
178 BIBLIOGRAPHY
[70] P. Poncet. Analysis of an immersed boundary method for three-

dimensional flows in vorticity formulation. J. Comput. Phys., 228:7268–
7288, 2009.
[71] P. M. Ricker. A direct multigrid poisson solver for oct-tree adaptive

meshes. The Astrophysical Journal Supplement Series, 176(1):293–300,
2008.
[72] D. Rossinelli, M. Bergdorf, G.-H. Cottet, and P. Koumoutsakos. Gpu

accelerated simulations of bluff body flows using vortex particle methods.
Journal of Computational Physics, 229(9):3316 – 3333, 2010.
[73] D. Rossinelli, M. Bergdorf, B. Hejazialhosseini, and P. Koumoutsakos.

Wavelet-based adaptive solvers on multi-core architectures for the simula-
tion of complex systems. 5704:721–734, 2009.
[74] J. K. Salmon and M. S. Warren. Skeletons from treecode closet. SIAM J.

Sci. Comput., 111(1):136–155, 1994.
[75] K. Shariff and A. Leonard. Vortex rings. Ann. Rev. Fluid Mech., 24:235,
1992.
[76] K. Shariff, R. Verzicco, and P. Orlandi. A numerical study of three dimen-

sional vortex ring instabilities : viscous corrections and early nonlinear
stage. J. Fluid Mech., 279:351, 1994.
[77] P. R. Spalart. Airplane trailing vortices. Ann. Rev. Fluid Mech., 30:107–
138, 1998.
[78] P.N. Swarztrauber. A direct method for the discrete solution of separable
elliptic equations. SIAM Journal of Numerical Analysis, 11:1136–1150,
1974.
[79] P.N. Swarztrauber, R. A. Sweet, and J. C. Adams. FISHPACK: Efficient

FORTRAN Subprograms for the Solution of Elliptic Partial Differential
Equations. Technical report, National Center for Atmospheric Research,
http://www.cisl.ucar.edu/css/software/fishpack/technote.ps, July 1999.
[80] P.N. Swarztrauber and R.A. Sweet. The fourier and cyclic reduction meth-
ods for solving Poisson’s equation. In J.A. Schetz and A.E. Fuhs, editors,
BIBLIOGRAPHY 179
Handbook of fluid dynamics and fluid machinery. John Wiley & Sons, New
York, 1996.
[81] W. Sweldens and P. Schröder. Building your own wavelets at home. In

Roland Klees and Roger Haagmans, editors, Wavelets in the Geosciences,
volume 90 of Lecture Notes in Earth Sciences, pages 72–107. Springer
Berlin / Heidelberg, 2000. 10.1007/BFb0011093.
[82] F. Thirifay and G. Winckelmans. Development of a lagrangian method for

combustion and application to the planar methane air jet diffusion flame.
Journal of Turbulence, 3(59):1468–5248, 2002.
[83] A. G. Tomboulides and S. A. Orszag. Numerical investigation of transi-

tional and weak turbulent flow past a sphere. J. Fluid Mech., 416:45–73,
2000.
[84] T. Uchiyama and Y. Yoshii. Improvements of vortex in cell method for

incimpressible flow simulation, submitted paper. Computers and Fluids,
2011.
[85] H.S. Udaykumar, R. Mittal, P. Rampunggoon, and A. Khanna. A Sharp

Interface Cartesian Grid Method for Simulating Flows with Complex Mov-
ing Boundaries. J. Comput. Phys., 174(1):345–380, 2001.
[86] X. Waerenburgh. Simulation numérique d’écoulements 2d incompressibles

autour de corps : Couplage d’une méthdode eulérienne et d’une méthode
langragienne. Travail de fin d’étude, Université catholique de Louvain
(UCL), Louvain-la-Neuve, 2001.
[87] G. Winckelmans. Vortex Methods, volume 3 (Fluids) of Encyclopedia of

Computational Mechanics, chapter 5. John Wiley & Sons, 2004.
[88] G. Winckelmans, R. Cocle, L. Dufresne, and R. Capart. Vortex methods

and their application to trailing wake vortex simulations. Comptes Rendus
Physique, 6(4-5):467–486, 2005. Aircraft trailing vortices.
[89] G. Winckelmans, L. Dufresne, and L. Bricteux. LES investigation of air-

craft wake two-vortex system in low atmospheric turbulence. Bull. Amer.
Phys. Soc., 50(10), 2005.
180 BIBLIOGRAPHY
[90] G. Winckelmans and A. Leonard. Contributions to vortex particle methods

for the computation of three-dimensional incompressible unsteady flows.
J. Comput. Phys., 2(109):247–273, 1993.
[91] G. Winckemans. Topics in Vortex Methods for the Computation of Three-

and Two-Dimensional Incompressible Unsteady Flows. PhD thesis, Cali-
fornia Institute of Technology, 1989.
[92] Y. Yang and H.S. Udaykumar. Sharp interface Cartesian grid method III:
Solidification of pure materials and binary solutions. J. Comput. Phys.,
210(1):55–74, 2005.
[93] L. Zhang and J. Eldredge. A viscous vortex particle method for deform-
ing bodies with application to biolocomotion. International Journal for
Numerical Methods in Fluids, 59(12):1299–1320, 2009.
Appendix A
Space-developing
simulation of a near-wake
rollup using a modeled
inflow vorticity field
This appendix presents the results of the space-developing simulation of the

vortex sheet rollup with velocity deficit using an inflow. It is shown that this
simulation provided non-physical results as the circulation of each vortex of the
counter-rotating vortex-pair increase.
The straightforward space-developing extension of the time-developing ve-
locity deficit model is here used as an inflow. It basically consists in imposing
the vorticity field (such as described in Chapter 4; in the paragraph entitled “
Description of the vorticity-field inflow-model ”)
A global 3-D perspective view is presented in Fig. A.1. Fig. A.2 provides a
zoom in the near-field region. Fig. A.4 presents various slices of the averaged
streamwise vorticity. At the first sight, one could be tempted to conclude that
the topology of the flow is similar to the one presented in Fig. 4.6 and Fig. 4.7.
Yet, this is not the case. It turns out that the circulation of the two rolled-
up vortices increases, see Fig. A.3. This also leads to a significant increase
182Appendix A. S-D simulation of a near-wake rollup using a modeled inflow
of the linear impulse. This behavior is definitely not physically acceptable.

That could be explained by analyzing the spanwise vorticity (i.e. the vorticity
component due the modeling of the friction drag) flux at the inflow. Because of
the presence of the “lifting-line” vorticity, the streamcise velocity component
is not up-down anti-symmetric (which was the case in the time-developing
simulation). Hence, the total spanwise vorticity flux (i.e. the integral over a
cross-plane) is none zero. The resulting difference is gathered in both vortices.
To circumvent this issue, another approach is proposed for which the lift-
force is directly modeled (exhaustively presented in Chapter 4).
Figure A.1: Space-developing simulation of the vortex-sheet rollup with velocity deficit using an inflow: 3D perspective view of two isocontours
183
of the vorticity norm kωkt0 = [200 400] colored with the axial vorticity component.
Figure A.2: Space-developing simulation of the vortex-sheet rollup with velocity

deficit using an inflow: zoom in the near-wake region of the 3D perspective view of
two isocontours of the vorticity norm kωkt0 = [200 400] colored with the axial vorticity
component.
0.8
[−]
0.6
0.4
0.2
0
0 0.05 0.1 0.15 0.2 0.25
τ

deficit using an inflow: evolution of the linear-impulse yc Γ+ /(Γ0 b) (solid) and of
the total circulation of the half-plane Γ+ /Γ0 (dash).
185

deficit using an inflow: 2D averaged vorticity field ωx t0 ; from top to bottom, time
t = [0.0, 0.07, 0.1, 0.3].
Appendix B
3D filament simulations of
Crow-like instabilities in
ground effect
This chapter was written within the FAR-Wake project as a technical report.
European project AST4-CT-2005-012238, Report TR3.1.1-1, February, 2006.
The authors are T. Lonfils, G. Daeninck and G. Winckelmans.
B.1 Introduction
The purpose of this work is to investigate two vortex systems in ground effect
(IGE) using a vortex filament method. The method being inviscid, the scope
of this study is however reduced: 2-D and 3-D (DNS and LES) results clearly
show that a very important effect for wake vortices (WV) IGE is the viscous
interaction with the ground through the creation of secondary vortices. We
therefore restrict the study to the analysis of long wave (Crow-type) instabilities
for WV IGE.
After briefly describing the numerical method, we investigate two cases: (i)
a configuration as close as possible to the one studied experimentally by CNRS-
IRPHE (i.e., a thick vortex pair IGE with a forced long wave instability) and
(ii) the case of aircraft like vortices IGE with a forced Crow instability (at
188 Appendix B. 3D filament simulations of Crow-like instabilities IGE
different initial perturbation levels).
B.2 Vortex filament method
In 3-D inviscid flows, vortex lines move as material lines: this constitutes the
basis for the method of vortex filaments. For more details on the vortex filament
method, and vortex methods in general, refer to [88, 87]. Each filament p
corresponds to a vortex tube of circulation Γp : the strength associated to
that filament. The method must be regularized, as singular filaments have
a logarithmically infinite self-induced velocity everywhere their curvature is
nonzero. For filaments using radially symmetric regularization functions, one
has Z
1 X gσ (|x − xp |)
uσ (x) = − Γp (x − xp ) × dxp
4π p Cp |x − x|3
△ r

where the regularization function gσ (r) = g σ is such that g(ρ) → 1 for ρ
3
large and g(ρ) ∝ ρ for ρ small. The corresponding streamfunction is obtained
as Z
1 X
Ψσ (x) = Γp Gσ (|x − xp |) dxp (B.1)
4π p Cp
△ 1 r
g(ρ)
with Gσ (r) = σG σ and ρ3 = − ρ1 dG
dρ (ρ). This streamfunction also corre-
2
sponds to the solution of ∇ ψσ = −ωσ with, as vorticity field,
X Z
ωσ (x) = Γp ζσ (|x − xp |) dxp (B.2)
p Cp
△ 1 r
Rρ
where ζσ (|x|) = 4πσ3 ζ σ . One has that g(ρ) = 0 ζ(s)s2 ds.
Typically, one uses parametric splines to represent numerically the fila-
ments. The regularized method converges for regular vorticity fields when the
number of filaments is increased, provided that the overlapping condition is
R
satisfied. The regularization function has normalization ζσ (x) dx = 1. It is
R
of order r when it satisfies the moment properties xp11 xp22 xp33 ζσ (x) dx = 0
R
for 1 ≤ p1 + p2 + p3 < r, and |x|r |ζσ (x)| dx < ∞. For radially symmetric
R∞
functions, these become 0 ζ(ρ) ρ2+s dρ = 0 for s even and 2 ≤ s < r, and
R∞ 2+r
0 |ζ(ρ)| ρ dρ < ∞. Order higher than r = 2 calls for functions that are
B.2. Vortex filament method 189
not strictly positive. See, Table B.1 for usual functions. The low order alge-
braic (LOA) function has r = 0: it can only be used for global vortex tube
modelling (as done in the present work), as it does not converge for detailed
field discretizations (using many filaments per vortex tube).
Figure B.1: Schematic of a infinite periodic filament and its image with respect to
the vertical plane and to the horizontal plane (ground plane). The peak plane and the
trough plane are also shown.
The filament method is here also “filtered”: indeed, if nothing is done, a high
wave number spurious mode will pollute the simulation (and eventually make
it blow up). Here, we use a discrete filter acting on the splines representation
of the filaments (as in [88]).
Figure B.1 describes the geometry of the problem. The problem being
periodic in the x-direction, we use a periodic version of the filament method.
ˆ to take into account the inviscid
Planes of symmetry are forced following Oxy,
ˆ to ensure the symmetrical evolution of the vortices.
wall, and following Oxz
We thus discretize and compute the evolution of one vortex filament only.
Table B.1: Examples of 3-D (g(ρ) and ζ(ρ)) regularization functions. The circulation
Γ(ρ) and vorticity ω(ρ) distributions for an equivalent 2-D vortex (i.e. straight 3-D
filament) are also provided.
g(ρ) ζ(ρ) Γ(ρ)/Γ0 πσ 2 Γ0 ω(ρ) r

3 2
ρ 3 ρ 1
LOA (ρ2 +1)3/2 (ρ2 +1)5/2 (ρ2 +1) (ρ2 +1)2
0
ρ3 (ρ2 +5/2) 15/2 ρ2 (ρ2 +2) 2
HOA (ρ2 +1)5/2 (ρ2 +1)7/2 2 2
(ρ2 +1)3
2
2 2
(ρ +1) 2 2
2 4
Gaussian erf (ρ) − ρ π1/2 e−ρ π 1/2
e−ρ 1 − e−ρ e−ρ 2
B.3 Definition
ˆ
One defines the Fourier transforms of the filament position in the plane Oxy
ˆ respectively:
and Oxz
N
1 X
ŷm = y(xi ) exp−2Iπmi/N
N i=1
N
1 X
ẑm = z(xi ) exp−2Iπmi/N
N i=1
where N the number of discrete points, m = 1, 2, . . . , N/2 and I 2 = −1. The

amplitude of each mode is:
||âm ||2 = ||ŷm ||2 + ||ŷm ||2
The growth rate is then defined by:
d
σ ∗ , σt0 = (log(||âm ||))
dt∗
where t∗ = t/t0 and t0 is the characteristic time scale. This corresponds to the
time for a pair of vortices characterized by an initial circulation Γ0 to descend
Γ0
OGE a distance b0 with the self-induced velocity V0 = 2πb0 .
B.4. Comparison: numerical simulation versus experiment 191
B.4 Comparison: numerical simulation versus

experiment
In the framework of FAR-wake Task 3.1.1, experiments have been performed
by CNRS-IRPHE [18]. One of the aims of this work is to compare numerical
results obtained using the filament method with experimental data. We first
compare visualizations of vortex locations at different times. We then analyze
the trajectories and temporal evolution of the vortex positions.
Experimental setup
The vortices are generated by two plates which are fixed to a common base and
are able to rotate. The motion of the plates is properly chosen so as to have
a spacing distance between the two vortices equals on average b0 = 2.4 [cm].
In order to force a longwave instability, the edges of the plates describe a
sinusoid. Here, the wavelength and the amplitude of this sinusoid is respectively
λf = 12 [cm] = 5.0 b0 and a0 = 1 [mm] = 0.0417 b0. The initial moment is
here defined when the vortex spacing no longer varies. This corresponds to
h0 = 6.2b0 . The experimental results show that the vortex circulation profile
is relatively well fitted by a Gaussian vortex characterized by σg = 0.17 b0 .
Numerical setup
In order to compare the numerical and experimental results, the numerical

parameters must be carefully chosen in order to match the experimental con-
figuration. The two vortices are separated by a distance b0 and are initialized
at an altitude of h0 = 6.2b0 . The perturbation plane is put at 45◦ . In this
plane, the perturbation is taken as:

x
a(x) = a0 sin 2π
λf
where a0 = 0.0417 b0 (it turns out that this initial amplitude produces results
that best match those of the experiment, see later) and λf = 5.0 b0 . Notice that
this wavelength is significantly shorter than the natural Crow instability (with
λC ≃ 8 − 9 b0 ); it is nevertheless the one that was forced in the CNRS-IRPHE
experiment and we here only aim at comparing with the experiment.

As mentioned above, the experimental vortex circulation profile can be
fitted by a Gaussian vortex. However, the 3-D Gaussian regularization function
is not implemented in the current filament code; this is because the Biot-Savart
evaluation is expensive. To be able to compare two equivalent vortex evolutions,
one uses vortices which have the same long wavelength dynamics.
In other words, in vortex filament methods, one can use vortices with an-
other distribution function and require to obtain: the same translational ve-
locity for thin vortex rings (Eq.B.3) and the same rotational velocity for long
wavelength perturbations of a vortex tube (Eq.(B.4)). As shown below, those
requirements are equivalent: the proper scaling fulfills both requirements.
• The self-induced translational velocity of a thin vortex ring (σ/R ≪ 1)

is given, for a regularized vortex by

Γ 8R
U= log −C (B.3)
4πR σ
with C = 1 for the low order algebraic function, C = 1/2 for the high
order algebraic function and C = (1 − γ/2) = 0.7114 for the Gaussian
function(where γ = 0.5772 is the Euler constant).
• The dispersion relation (self-induced rotational velocity) of long wave-

length perturbations (kσ ≪ 1) is given, for a regularized vortex by

Γ 2 kσ 1
Ω=± k log + γ− +C (B.4)
4π 2 2
with the same C values as above.
We here use the high order algebraic distribution and we scale σ to obtain the
same dynamics as those of a Gaussian vortex.
Another error comes from the fact that we here discretize a vortex tube
with only one vortex filament: because the velocity is computed at the center
of the filament, it is not equal to the true translational velocity of the vortex
tube. It can be corrected by using, again, an equivalent σ. Winckelmans [91]
showed that using C = 0.5580 instead of 0.7114 in Eq.(B.3) and Eq.(B.4) leads
to the correct dynamics of a Gaussian vortex characterized by σg .
B.4. Comparison: numerical simulation versus experiment 193
Thus, in order to obtain the same dynamics as those in the experiment, one
simply uses the filament method with a high order algebraic regularization and
with σhoa such that:
log σhoa + Choa = log σg + Cg
with Cg = 0.5580, Choa = 0.5 and here with σg = 0.17 b0. It leads to
σhoa = 1.060σg = 0.180 b0
The vortex tube is represented by one vortex filament and is discretized by 320
points over λf , and cubic splines. The time step is ∆t∗ = 10−2 /(2π). We use
the same filter parameters as those detailed in [88].
Results
Comparison with the visualizations The qualitative vortex evolution is

very similar to the experimental one. Figures B.2 and B.3 show the top and
side views at t∗ = 2.7 (during the linear growth phase) and t∗ = 3.8 (at
reconnection) respectively. Even if the Reynolds number of the experiment is
not high (Re ≃ 3500 and thus there is some diffusion of the vortex core), the
inviscid vortex filament method is in good agreement with the experimental
results up to the reconnection phase: the position of the vortices and the
amplitude of the instability are very similar.
Temporal evolution up to reconnection Experimental results are given

for characteristic sections: the peak and the trough sections: i.e., the sections
where the altitude are maximal and minimal respectively (Figure B.1). Up to
the reconnection, trajectory, altitude and vortex spacing evolution are qualita-
tively similar as shown in Figure B.4.(a) and Figure B.5. The growth rate of
the long wavelength instability is obtained as σ ∗ = 0.66, as shown in Figure
B.4.(b).
(a)
(b)
Figure B.2: Comparison of the vortex wake evolution: numerical case (left) at t∗ =
2.7 and experimental case (right) at t∗ ≃ 2.7. (a) top view and (b) side view, at same
scale. The tube radius shown in the simulation is that for which the vortex induced
velocity is maximal: r = rc = 1.121 σg .
B.5 Analysis of the Crow instability in ground

effect
Numerical simulations have been performed to investigate the behavior of the

natural Crow instability in ground effect (before reconnection). Crow [20] has
shown that the most unstable long wavelength mode of a pair of wake vortices
(σ/b0 ≪ 1) has a wavelength: λ = 8.6 b0 . This mode grows in a plane inclined
at about 48◦ to the horizontal. The growth rate is σ ∗ = 0.83. The purpose of
this section is to investigate how the inviscid ground affects the Crow instability
in terms of the deformation plane’s angle and of the growth rate.
Numerical setup The following simulations are performed using the low
order algebraic regularization function (see Table B.1) as it is a good fit usual
for aircraft wake vortices. The value of σ is fixed to match a typical wake
vortex core σ = rc = 0.05 b0. The vortices are initialized at an altitude out
B.5. Analysis of the Crow instability in ground effect 195
(a)
(b)
Figure B.3: Comparison of the vortex evolution: numerical case (left) at t∗ = 3.8
and experimental case (right) at t∗ ≃ 3.9. (a) top view and (b) side view, at same
scale. In the simulation, the tube shown is that for which r = rc .
(a) (b)
7
0
10
6
5
||â||
4
z/b
0 σ*=0.66
−1
3 10
−2
0 10
−3 −2 −1 0 1 2 3 0 1 2 3 * 4 5 6
y/b0 t
Figure B.4: Trajectory of the vortices (a) and evolution of the perturbation (b). The
peak (solid line) and the trough (dash-dotted line).
of ground effect: h0 = 6 b0 . The perturbation is put in a plane inclined at

48.3◦ (accurate value determined by simulation carried OGE) to the horizontal
and is initialized with a wavelength of λC = 8.6 b0 . Different perturbation
(a) (b)
4 6
5
3
4
b/b0 z/b
0
2 3
2
1
1
0 0
0 1 2 3 4 5 6 0 1 2 3 4 5 6
t* t*
Figure B.5: Evolution of the spacing between the two vortices (a) and of the vortex
altitudes (b). The peak (solid line), the trough (dash-dotted line) and the average
(dotted line).
levels were investigated: a0 /b0 = [10−1 , 10−2 , 10−3 , 10−4 , 10−5 , 10−6 ]. The
filament is discretized using 512 points over λC , and cubic splines. The time
step is ∆t∗ = 5 10−3 /(2π) and the filter parameters is again the same as those
detailed in [88].
Instability analysis Figure B.6 and Figure B.7 show the visualization of
the filament evolutions for two initial perturbation levels: a0 /b0 = 10−2 and
a0 /b0 = 10−3 respectively. One can see that there are essentially two cases:
(i) when the initial perturbation is high enough, the vortex filaments reconnect
before entering IGE, as in Figure B.6 or (ii) for low initial perturbation level, the
vortices enter IGE before reconnection and interact (and eventually reconnect)
with their respective image.
Figures B.8 and B.9 describe the evolution of the amplitude and angle of
the Crow mode for the different initial perturbation levels. Out of ground effect
(t∗ < 4), the growth rate of the Crow mode is σ ∗ = 0.82 and the perturbation
angle is 48.3◦ . When the vortices are influenced by their respective images only
(t∗ > 9), the growth rate of the Crow mode is again σ ∗ ≃ 0.82; the perturbation
angle is then −41.7◦ (thus 48.3◦ with respect to the vertical). Between these
two phases, the angle of the perturbation plane changes (Figure B.9). Figure
B.10 further shows this phenomenon: ||ẑm || changes sign while ||ŷm || slows
down, stops, and then grows again.
As expected, if the initial perturbation level is low enough, the results show
that the evolution of the Crow instability IGE is universal and independent
of the initial perturbation amplitude. The rotation of the perturbation plane
also happens at the same altitude z/b0 ≃ 0.53 above the ground (here at the
same time t∗ ≃ 7.13 as all simulations were initially at h0 /b0 = 6). Figure B.11
shows the trajectories of the vortices and the orientation of the perturbation
planes. As expected, out of ground effect, and completely in ground effect, the
instability characteristics (growth rate and angle) are those of the usual Crow
instability.
(a) t∗ = 0 (c) t∗ = 3.6
(b) t∗ = 2.0 (d) t∗ = 4.7
Figure B.6: Visualization of the vortex wake evolution. The shadow is also shown.
The initial perturbation level is: a0 /b0 = 10−2 .
(a) t∗ = 0 (d) t∗ = 6.4
(b) t∗ = 2.0 (e) t∗ = 8.0
(c) t∗ = 3.6 (f) t∗ = 9.8
Figure B.7: Visualization of the vortex wake evolution. The shadow is also shown.
The initial perturbation level is: a0 /b0 = 10−3 .
0
10
−1
10
σ*=0.82
−2
10
||â||
−3
10
*
−4 σ =0.82
10
−5
10
−6
10
0 2 4 6 8 10 12 14 16 18 20
t*
Figure B.8: Temporal evolution of the amplitude for the Crow mode λC /b0 = 8.6
for different initial perturbation levels.
60
10−1 10−2
40
20
[°]
0
−20
10−3 10
−4
10−5 10−6
−40
0 2 4 6 8 10 12 14 16 18
t*
Figure B.9: Temporal evolution of the perturbation plane angle for each initial per-
turbation level a0 /b0 . The time at which the orientation changes sign is t∗ ≃ 7.13,
which corresponds to z/b0 = 0.53.
−2
10
||â||
−4
10
−6
10
0 2 4 6 t* 8 10 12 14
Figure B.10: Temporal evolution of the amplitude for the Crow mode λC /b0 = 8.6.
The initial perturbation level is: a0 /b0 = 10−6 . ||â|| (solid line), ||ŷ|| (dashed line)
and ||ẑ|| (dotted line).
7
b
0
y/b0
3
b
0
0
−1
−4 −3 −2 −1 0 1 2 3 4
x/b
0
Figure B.11: Trajectories of the vortices and orientation of the perturbation planes
(solid lines). Planes of symmetry (dash-dotted line). The time between each orienta-
tion shown is ∆t∗ = 0.4.

Thesis

Uploaded by

Copyright:

Available Formats

Thesis

Uploaded by

Document Information

Original Title

Copyright

Available Formats

Share this document

Share or Embed Document

Sharing Options

Did you find this document useful?

Is this content inappropriate?

Copyright:

Available Formats

Thesis

Uploaded by

Copyright:

Available Formats

"Vortex particle-mesh with combined immersed

boundary and mesh refinement techniques :

Document type : Thèse (Dissertation)

Vortex particle-mesh method with combined immersed

Composition du jury: Thèse présentée en vue de

Je voudrais avant tout remercier mon promoteur, Grégoire Winckelmans, qui

fructueuse. Il a été, de facto, mon troisième promoteur!

2 An IB technique within the combined VPM-PFM method 9

3 A mesh refinement technique for the VPM method 37

4 Investigation of velocity deficit effects on aircraft wake rollup 69

5 Numerical investigations of end-effect 113

6 Conclusions & perspectives 161

A S-D simulation of a near-wake rollup using a modeled inflow 181

B 3D filament simulations of Crow-like instabilities IGE 187

1.2 Some preliminaries

Formulation equations The Navier-Stokes equations for an incompressible

where u is the velocity field and ω , ∇ × u is the vorticity field, thus ∇ · ω = 0.

where Ψ is the streamfunction and φ is the scalar potential. We take φ such

Discretization The vorticity field is discretized using vortex particles which

• Combining a fast grid-based Poisson solver together with the handling of

• Efficiently dealing with boundary layers. Vortex particle methods are

ber and ii) including a multi-resolution approach in the VPM-PFM method.

• Chapter 2 presents the numerical algorithm developed to compute the

• Chapter 3 deals with the integration of a multiple space-resolution tech-

The second part of the manuscript is devoted to the numerical investigation of

• Chapter 5 investigates the “end-effect” phenomena. The end-effect oc-

Finally, the effect of an inviscid ground on the Crow(-like) instability growth is

Abstract This paper presents a method for the simulation of incompress-

Boundary (IB) approach. The Navier-Stokes equations in the vorticity-velocity

means of a panel method, thus in a step reminiscent of a sharp interface method

Vortex methods solve the Navier-Stokes equations in the vorticity-velocity

The velocity is then computed using the relation

This chapter is structured as followed. First, the algorithm is described in

2.2 An Immersed Boundary Method Within a

Figure 2.1: General setup of the IB method.

The algorithm reads: At time tn , the vorticity field is discretized using a

1. Mesh overlay & Random grid shifting. A Cartesian grid with a

uniform resolution h, is laid over the support of vorticity, i.e. it covers

3. Evaluation of the vorticity field on the grid. We construct a set P2

is here done using one-side second-degree polynomial functions P:

4. No-slip condition. The vortex panel vector intensity is computed such

5. Vorticity flux. A Lighthill’s model of vorticity flux at the wall is

ξ and η are the local integration coordinates. hl /2 is equal to zk if

In practice, if the triangular panels are equilateral enough, the panels

6. Poisson equation. Solving the Poisson equation ∇2 Ψnijk = −(ωijk

7. Evolution equation. The RHS of Eq.(2.1) and (2.2) are computed

The vorticity field thus has to be regularly corrected in order to project it

The vorticity field ω can be corrected to a divergence free field ω̃:

by taking D such that

9. Interpolation. Grid quantities (i.e. velocity, time-derivative of strength

After particle redistribution (typically every 5 time steps), the AB in-

x∗p = xnp + ∆tunp

and the second step (corrector) is

2.3 Validation: flow past a sphere at Re = 300

Figure 2.3: General view of the flow past an arbitrarily shape.

Figure 2.4: Sphere discretization using recursive subdivision of an icosahedron.

One defines the Drag coefficient (CD )

The time-evolution of the diagnostics (drag/lift coefficients) are presented in

CD ∆CD CL ∆CL Str

Figure 3.11: Multi-resolution