# Reviving the Vortex Particle Method:
A Stable Formulation for Meshless Large Eddy Simulation

Eduardo J. Alvarez<EMAIL_ADDRESS>Andrew Ning<EMAIL_ADDRESS>Ira A.
Fulton College of Engineering, Brigham Young University, Utah, USA

###### Abstract

The vortex particle method (VPM) is a mesh-free approach to computational
fluid dynamics (CFD) solving the Navier-Stokes equations in their velocity-
vorticity form. The VPM uses a Lagrangian scheme, which not only avoids the
hurdles of mesh generation, but it also conserves vortical structures over
long distances with minimal numerical dissipation while being orders of
magnitude faster than conventional mesh-based CFD. However, VPM is known to be
numerically unstable when vortical structures break down close to the
turbulent regime. In this study, we reformulate the VPM as a large eddy
simulation (LES) in a scheme that is numerically stable, without increasing
its computational cost. A new set of VPM governing equations are derived from
the LES-filtered Navier-Stokes equations. The new equations reinforce
conservation of mass and angular momentum by reshaping the vortex elements
subject to vortex stretching. In addition to the VPM reformulation, a new
anisotropic dynamic model of subfilter-scale (SFS) vortex stretching is
developed. This SFS model is well suited for turbulent flows with coherent
vortical structures where the predominant cascade mechanism is vortex
stretching. Extensive validation is presented, asserting the scheme comprised
of the reformulated VPM and SFS model as a meshless LES that accurately
resolves large-scale features of turbulent flow. Advection, viscous diffusion,
and vortex stretching are validated through simulation of isolated and
leapfrogging vortex rings. Mean and fluctuating components of turbulent flow
are validated through simulation of a turbulent round jet, where Reynolds
stresses are resolved directly and compared to experimental measurements.
Finally, the computational efficiency of the scheme is showcased in the
simulation of an aircraft rotor in hover, showing our meshless LES to be 100x
faster than a mesh-based LES with similar fidelity, while being 10x faster
than a low-fidelity unsteady Reynolds-average Navier-Stokes simulation and
1000x faster than a high-fidelity detached-eddy simulation. The implementation
of our meshless LES scheme is hereby released as an open-source software,
called FLOWVPM.

1.125(-2.25cm,-13.9cm)  Preprint of E. J. Alvarez and A. Ning, “Stable Vortex
Particle Method Formulation for Meshless Large-Eddy Simulation,” AIAA Journal,
Nov 2023. doi: 10.2514/1.J063045
The content of this paper may differ from the final publisher version.

##  I Introduction

Large eddy simulation (LES) is a class of computational fluid dynamics (CFD)
that filters the Navier-Stokes equations to decompose small-scale fluctuations
of the flow from large ones. The large scales are then resolved directly,
while the effects of the smaller scales are modeled. The most common LES
methods use a mesh or grid to discretize the space and calculate fluxes and
derivatives, classified as finite volume/element/difference methods.
Significant user effort is spent in the generation and manipulation of the
mesh, with studies showing that about 67% of engineering time in mesh-based
CFD is spent in these efforts.Hardwick  _et al._ (2005)

Vortex methods Leonard (1980); Cottet and Koumoutsakos (2000a); Winckelmans
(2004); Mimeau and Mortazavi (2021) are a class of meshless CFD solving the
Navier-Stokes equations in their velocity-vorticity form. This form is solved
in a Lagrangian scheme, which not only avoids the hurdles of mesh generation,
but also conserves vortical structures over long distances with minimal
numerical dissipation.

The vortex particle method Winckelmans and Leonard (1993) (VPM) is a vortex
method that uses particles to discretize the Navier-Stokes equations, with the
particles representing radial basis functions that construct a continuous
vorticity field Barba (2004). This meshless CFD has several advantages over
conventional mesh-based CFD. In the absence of a mesh, the VPM (1) does not
suffer from the conventional Courant–Friedrichs–Lewy (CFL) condition, (2) does
not suffer from the numerical dissipation introduced by a mesh, (3)
derivatives are calculated exactly rather than approximated through a stencil,
and (4) it has been shown to be 100x to 1000x faster than mesh-based CFD with
comparable fidelity Alvarez  _et al._ (2020). Furthermore, the VPM is
spatially second-order accurate Yokota and Barba (2013), it is highly
efficient since elements are placed only where vorticity exists, the spatial
discretization can be automatically adapted Barba, Leonard, and Allen (2005a);
Lakkis and Ghoniem (2009); Stock and Gharakhani (2021) in the fashion of an
adaptive mesh refinement, and simulations are highly-parallelizable in
heterogeneous CPU and GPU architectures Yokota  _et al._ (2013); Hu  _et al._
(2014).

Figure 1: Meshless LES of eVTOL aircraft using reformulated VPM: (left)
computational elements (vortex particles and strength), and (right) volume
rendering of vorticity field.

VPM has gained popularity in recent years due to a growing need to predict
complex aerodynamic interactions in modern electric aircraft Kasliwal  _et
al._ (2019); nas (2020); Schwab  _et al._ (2021); McDonald  _et al._ (2021) at
a computational cost fit for conceptual design Silva  _et al._ (2018); Brelje
and Martins (2018). For example, VPM has been used for high-fidelity
simulation of rotorcraft forward flight Stock, Gharakhani, and Stone (2010)
and multirotor interactions Alvarez and Ning (2020); Alvarez  _et al._ (2020),
and as a mid and low-fidelity tool for ground effect Tan  _et al._ (2020,
2021), electrical vertical takeoff and landing (eVTOL) Arda Yücekayali (2020),
distributed electric propulsion Alvarez and Ning (2018); Teixeira and Cesnik
(2019), aeromechanics of unconventional rotorcraft Singh and Friedmann (2021);
Jacobellis  _et al._ (2021); Corle, Floros, and Schmitz (2021), rotor-rotor
interactions Alvarez and Ning (2019); Lee and Lee (2020), and wind energy Lee
and Lee (2019a, b); Mehr, Alvarez, and Ning (2020). Furthermore, a couple of
opensource VPM codes have been recently released Kuzmina, Marchevsky, and
Ryatina (2019); Stock and Gharakhani (2020); Tugnoli  _et al._ (2021). Fig. 1
shows the simulation of an eVTOL aircraft using VPM.

In spite of its growing popularity, VPM is known to be numerically unstable
when vortical structures break down close to the turbulent regime. This has
limited its range of applications in the aforementioned studies to mostly
benign cases with well-behaved numerics (e.g., coarse simulations of rotors
with axial inflow). The instability is triggered when vortex stretching gives
raise to a rapid increase of local vorticity, which the low numerical
dissipation of the method fails to damp out. Such numerical instability is
believed to be caused by a combination of Lagrangian distortion, a vorticity
field that is not divergence-free, and the absence of subfilter-scale (SFS)
turbulent diffusion, among other reasons.

Multiple meshless schemes have been developed over the years to address
Lagrangian distortion Barba, Leonard, and Allen (2005b); Lakkis and Ghoniem
(2009); Kirchhart and Rieger (2021), divergence of the vorticity field
Pedrizzetti (1992); Gharakhani (1997), and SFS effects Gharakhani (2002,
2005). While some of these schemes are effective in two-dimensional cases,
they have not succeeded at making the three-dimensional VPM numerically
stable. The need to circumvent this challenge motivated further development of
the vortex-in-cell (VIC) method Christiansen (1997); Cottet and Poncet (2004);
Koumoutsakos (2005), also known as vortex particle-mesh method. In this
method, the particles are projected onto a background mesh at every time step,
and vortex stretching, viscous diffusion, and the Biot-Savart law are computed
in a mesh-based scheme. This approach has shown to be numerically stable,
recently enabling the study of wake dynamics with unprecedented
fidelity.Chatelain  _et al._ (2008, 2017); Caprace, Chatelain, and Winckelmans
(2020); Nguyen (2019); Nguyen and Duong (2021); Colognesi, Ronsse, and
Chatelain (2021); Ramos-García  _et al._ (2021) However, the introduction of
this mesh also forfeits a few of the aforementioned benefits of a purely
Lagrangian (meshless) scheme.

In this study, we propose an LES formulation of the VPM that is both
numerically stable and meshless. A new set of VPM governing equations are
derived from the LES-filtered Navier-Stokes equations in Section II. In
Section III, the classic VPM is shown to be only one out of multiple possible
formulations of these governing equations. Such formulation is shown to
locally violate conservation of mass and angular momentum under vortex
stretching, explaining why the classic VPM tends to be numerically unstable.
In Section IV, we consider multiple possible formulations using different
element shapes and conservation laws. One formulation is selected, which uses
the particle size to reinforce local conservation of both mass and angular
momentum. This formulation, referred to as the reformulated VPM, is
implemented and hereby released as the open-source software FLOWVPM.***The
code is available at https://github.com/byuflowlab/FLOWVPM.jl

In addition to the reformulated VPM, a novel anisotropic structural model of
subfilter-scale (SFS) vortex stretching is developed in Section V. The means
for backscatter control are also provided, along with a dynamic procedure for
the automatic computation of the model coefficient. This SFS model is apt for
both meshless and mesh-based CFD, and is well suited for turbulent flows with
coherent vortical structures where the predominant cascade mechanism is vortex
stretching.

The scheme comprised of the reformulated VPM and SFS model is extensively
tested in Section VI, validating the scheme as an LES method. Advection,
viscous diffusion, and vortex stretching are validated through simulation of
isolated and leapfrogging vortex rings. Large-scale turbulent dynamics are
validated through simulation of a turbulent round jet, where predicted
fluctuations and Reynolds stress are compared to experimental measurements.
Finally, the computational efficiency of our meshless LES is showcased in an
engineering application through the simulation of an aircraft rotor in hover.

##  II Fundamentals of
the Vortex Particle Method

A new set of VPM governing equations are hereby derived from the LES-filtered
Navier-Stokes equations in their vorticity form. Along the way, the derivation
here presented introduces all the fundamental concepts necessary to bring the
reader up to speed with LES and vortex methods.

###  II.1 Vorticity Form of the Navier-Stokes Equations

In a Newtonian, incompressible fluid with constant viscosity, the linear
momentum of a differential fluid element is governed by the Navier-Stokes
equation

$\displaystyle\frac{\partial{\textbf{u}}}{\partial
t}+({\textbf{u}}\cdot\nabla){\textbf{u}}=-\frac{1}{\rho}\nabla
p+\nu\nabla^{2}{\textbf{u}},$ (1)

where ${\textbf{u}}({\textbf{x}},t)$ is the velocity field,
$p({\textbf{x}},t)$ is the pressure field, and $\rho$ and $\nu$ are the
density and kinematic viscosity of the fluid, respectively.

Taking the curl over Eq. 1, the pressure dependence disappears obtaining

$\displaystyle\frac{\text{D}}{\text{D}t}\boldsymbol{\omega}=(\boldsymbol{\omega}\cdot\nabla){\textbf{u}}+\nu\nabla^{2}\boldsymbol{\omega},$
(2)

where $\frac{\text{D}}{\text{D}t}$ denotes the material derivative operator,
$\frac{\text{D}}{\text{D}t}()\equiv\frac{\partial}{\partial
t}()+{({\textbf{u}}\leavevmode\nobreak\ \cdot\leavevmode\nobreak\ \nabla)()}$,
and
$\boldsymbol{\omega}({\textbf{x}},t)=\nabla\times{\textbf{u}}({\textbf{x}},t)$
is the vorticity field.

Eq. 2 is the vorticity form of the Navier-Stokes linear-momentum equation.
This equation depends on $\boldsymbol{\omega}$ alone since u can be calculated
from $\boldsymbol{\omega}$ through the Biot-Savart law. While Eq. 1 stems from
conservation of linear momentum, Eq. 2 can be interpreted as some form of
conservation of angular momentumChatwin (1973), which is further discussed in
Section IV.1. From the right-hand side, we see that the evolution of the
vorticity field is governed by vortex stretching (first term) and viscous
diffusion (second term).

###  II.2 Large Eddy Simulation Equation

We now derive the filtered vorticity Navier-Stokes equation—here referred to
as large-eddy simulation (LES) equation—and define the stress tensor that
transfers energy between resolved scales and subfilter scales.

Let $\phi$ be a field and $\zeta_{\sigma}$ a filter kernel, the filtered field
is denoted by a bar and defined as

$\displaystyle\overline{\phi}\left(\mathbf{x}\right)$
$\displaystyle\equiv\int\limits_{-\infty}^{\infty}\phi(\mathbf{y})\zeta_{\sigma}(\mathbf{x}-\mathbf{y})\,\mathrm{d}\mathbf{y},$

where the filter $\zeta_{\sigma}$ is associated to a certain cutoff length
scale $\sigma$. In this study, $\zeta_{\sigma}$ is defined as
$\zeta_{\sigma}(\mathbf{x})\equiv\frac{1}{\sigma^{3}}\zeta(\frac{\|\mathbf{x}\|}{\sigma})$
where $\zeta$ is a radial basis function, and is required to have a volume
integral of unity,
$\int\limits_{-\infty}^{\infty}\zeta_{\sigma}(\mathbf{y})\,\mathrm{d}\mathbf{y}=1$.

In order to derive the LES version of the vorticity equation, Eq. 2 is written
in tensor notation and filtered as

$\displaystyle\overline{\frac{\partial\omega_{i}}{\partial
t}}+\overline{u_{j}\frac{\partial\omega_{i}}{\partial
x_{j}}}=\overline{\omega_{j}\frac{\partial u_{i}}{\partial
x_{j}}}+\nu\overline{\nabla^{2}\omega_{i}}.$ (3)

In this equation, $\overline{u_{j}\frac{\partial\omega_{i}}{\partial x_{j}}}$
and $\overline{\omega_{j}\frac{\partial u_{i}}{\partial x_{j}}}$ are non-
linear terms that cannot be calculated from resolved quantities, but are
rather approximated through a tensor $T_{ij}$ that encapsulates the error
between $\overline{u_{i}\omega_{j}}$ and
$\overline{u_{i}}\,\overline{\omega_{j}}$ as

$\displaystyle\overline{u_{i}\omega_{j}}=\overline{u_{i}}\,\overline{\omega_{j}}+T_{ij}.$

The gradient of $T_{ij}$ and its transpose result in

$\displaystyle\frac{\partial T^{\prime}_{ij}}{\partial
x_{j}}=\overline{u_{j}\frac{\partial\omega_{i}}{\partial
x_{j}}}-\overline{u_{j}}\frac{\partial\overline{\omega_{i}}}{\partial x_{j}}$
$\displaystyle\frac{\partial T_{ij}}{\partial
x_{j}}=\overline{\omega_{j}\frac{\partial u_{i}}{\partial
x_{j}}}-\overline{\omega_{j}}\frac{\partial\overline{u_{i}}}{\partial x_{j}}$

Replacing this into Eq. 3 and using the derivative-filter commutation
property, we obtain the LES vorticity equation:

$\displaystyle\frac{\partial\overline{\omega_{i}}}{\partial
t}+\overline{u_{j}}\frac{\partial\overline{\omega_{i}}}{\partial x_{j}}$
$\displaystyle=$ (4)
$\displaystyle\overline{\omega_{j}}\frac{\partial\overline{u_{i}}}{\partial
x_{j}}$ $\displaystyle+\nu\nabla^{2}\overline{\omega_{i}}-\frac{\partial
T^{\prime}_{ij}}{\partial x_{j}}+\frac{\partial T_{ij}}{\partial x_{j}}.$

The term $\frac{\partial T^{\prime}_{ij}}{\partial x_{j}}$ represents the
subfilter-scale (SFS) contributions arising from the advective term (vorticity
advection), while $\frac{\partial T_{ij}}{\partial x_{j}}$ represents the
contributions arising from vortex stretching. $T_{ij}$ is associated to the
SFS vorticity stress that encapsulates the interactions between large scale
dynamics and SFS dynamics and has to be modeled in terms of resolved
quantities. The accuracy of LES hinges on the modeling of this tensor. Its
divergence represents the rate at which enstrophy—a measure of rotational
kinetic energy—is transferred from resolved scales to subfilter scales
(diffusion) and from subfilter scales to resolved scales (backscatter). In
Section V we will develop a model that approximates the vortex-stretching
component of this tensor, and provide the means for backscatter control.

For convenience, we write Eq. 4 in vector notation as

$\displaystyle\frac{\text{d}}{\text{d}t}\overline{\boldsymbol{\omega}}=\left(\overline{\boldsymbol{\omega}}\cdot\nabla\right)\overline{\mathbf{u}}+\nu\nabla^{2}\overline{\boldsymbol{\omega}}-\mathbf{E}_{\mathrm{adv}}-\mathbf{E}_{\mathrm{str}},$
(5)

where $\left(\mathbf{E}_{\mathrm{adv}}\right)_{i}\equiv\frac{\partial
T^{\prime}_{ij}}{\partial x_{j}}$ is the SFS vorticity advection,
$\left(\mathbf{E}_{\mathrm{str}}\right)_{i}\equiv-\frac{\partial
T_{ij}}{\partial x_{j}}$ is the SFS vortex stretching, and the
$\frac{\text{d}}{\text{d}t}$ operator is a linearized version of the filtered
material derivative,
$\frac{\text{d}}{\text{d}t}()\equiv\frac{\partial}{\partial
t}()+{(\overline{{\mathbf{u}}}\cdot\nabla)()}$.

###  II.3 Viscous Diffusion

One common practice for solving non-linear partial differential equations
(PDE), like Eq. 5, is that of splitting the PDE operator into a linear sum of
its non-linear pieces. This permits discretizing and solving each non-linear
piece in a separate numerical scheme. In this study we will split the PDE in
Eq. 5 into inviscid and viscous pieces as

$\displaystyle\frac{\text{d}}{\text{d}t}\overline{\boldsymbol{\omega}}$
$\displaystyle=\left(\frac{\text{d}}{\text{d}t}\overline{\boldsymbol{\omega}}\right)_{\mathrm{inviscid}}+\left(\frac{\text{d}}{\text{d}t}\overline{\boldsymbol{\omega}}\right)_{\mathrm{viscous}},$

where

$\displaystyle\left(\frac{\text{d}}{\text{d}t}\overline{\boldsymbol{\omega}}\right)_{\mathrm{inviscid}}$
$\displaystyle=\left(\overline{\boldsymbol{\omega}}\cdot\nabla\right)\overline{\mathbf{u}}-\mathbf{E}_{\mathrm{adv}}-\mathbf{E}_{\mathrm{str}}$
$\displaystyle\left(\frac{\text{d}}{\text{d}t}\overline{\boldsymbol{\omega}}\right)_{\mathrm{viscous}}$
$\displaystyle=\nu\nabla^{2}\overline{\boldsymbol{\omega}}.$

Over the years, multiple Lagrangian schemes have been developed that
accurately resolve the viscous component, like the vortex redistribution
method Shankar and Dommelen (1996); Gharakhani (2001), particle strength
exchange Degond and Mas-Gallic (1989), and core spreading Rossi (1996), to
name a few. In this study, viscous diffusion will be solved through the core
spreading method coupled with the radial basis function interpolation approach
for spatial adaptation developed by Barba Barba, Leonard, and Allen (2005b);
Barba (2004); Torres and Barba (2009). This viscous scheme diffuses the
vorticity by thickening each particle’s core size $\sigma$ over time, while
using an RBF interpolation to reset core sizes when they have overgrown. As
shown by Rossi Rossi (1996), the core spreading method has second-order
spatial convergence, while showing linear convergence when coupled with
spatial adaptation.

In the following sections we will focus on the inviscid part of the PDE,
developing a scheme for its numerical solution.

###  II.4 Lagrangian Discretization: The Vortex Particle

Figure 2: Vortex particle with unit size and strength: (left) core size
$\sigma$ and vortex strength $\Gamma$, (middle) vorticity field, and (right)
contours of velocity field and streamlines (represented as curved arrows).

The material derivative in Eq. 2 and the material-conservative nature of the
vorticity makes the $\boldsymbol{\omega}$ field especially well fit for a
Lagrangian description. We now discretize the vorticity equation with
Lagrangian elements, termed vortex particles. Each particle represents a
volume of fluid that is convected by the velocity field carrying an integral
quantity of vorticity.

The unfiltered $\boldsymbol{\omega}$ field is discretized with singular vortex
particles of positions ${\textbf{x}}_{p}$ and coefficients
$\boldsymbol{\Gamma}_{p}$, approximating $\boldsymbol{\omega}$ as

$\displaystyle\boldsymbol{\omega}({\textbf{x}},t)\approx\sum\limits_{p}\boldsymbol{\Gamma}_{p}(t)\delta({\textbf{x}}-{\textbf{x}}_{p}(t)),$
(6)

where $\delta$ is the Dirac delta. Each particle travels with the local
velocity as

$\displaystyle\frac{\text{d}}{\text{d}t}{\textbf{x}}_{p}={\textbf{u}}({\textbf{x}}_{p}),$

where ${\textbf{x}}_{p}$ is the position of the $p$-th particle. Thus, each
coefficient $\boldsymbol{\Gamma}_{p}$ (termed vortex strength) represents the
average vorticity that is carried in the volume of each particle since

$\displaystyle\int\limits_{-\infty}^{\infty}\boldsymbol{\omega}({\textbf{x}},t)\,\mathrm{d}\mathbf{x}\approx\sum\limits_{p}\boldsymbol{\Gamma}_{p}(t).$

Using the singular particle approximation in the filtered vorticity
$\overline{\boldsymbol{\omega}}$,

$\displaystyle\overline{\boldsymbol{\omega}}\left(\mathbf{x}\right)$
$\displaystyle=\int\limits_{-\infty}^{\infty}\boldsymbol{\omega}\left(\mathbf{y}\right)\zeta_{\sigma}(\mathbf{x}-\mathbf{y})\,\mathrm{d}\mathbf{y}$
$\displaystyle\approx\int\limits_{-\infty}^{\infty}\left(\sum\limits_{p}\boldsymbol{\Gamma}_{p}\delta({\textbf{y}}-{\textbf{x}}_{p})\right)\zeta_{\sigma}(\mathbf{x}-\mathbf{y})\,\mathrm{d}\mathbf{y},$

we obtain an approximation of the filtered vorticity field,

$\displaystyle\overline{\boldsymbol{\omega}}\left(\mathbf{x}\right)\approx\sum\limits_{p}\boldsymbol{\Gamma}_{p}\zeta_{\sigma}(\mathbf{x}-\mathbf{x}_{p}).$
(7)

As seen in Eq. (7), the filter operator has the effect of spreading the vortex
strength $\boldsymbol{\Gamma}_{p}$ in space, regularizing the singularity
originally introduced by the Dirac delta in Eq. (6). Thus, the filter kernel
takes the role of a basis function that is used to discretize and approximate
the filtered vorticity field through particles. We let the filter width
$\sigma$ change in time and space according to the evolution of each
individual particle. Here on, the filter width is referred to as smoothing
radius or core size, denoted $\sigma_{p}$, and defined as
$\sigma_{p}(t)\equiv\sigma(\mathbf{x}_{p},t)$. We approximate the filtered
vorticity as
$\overline{\boldsymbol{\omega}}\approx\boldsymbol{\omega}_{\sigma}$, with

$\displaystyle\boldsymbol{\omega}_{\sigma}\left(\mathbf{x},\,t\right)$
$\displaystyle\equiv\sum\limits_{p}\boldsymbol{\Gamma}_{p}(t)\zeta_{\sigma_{p}}(\mathbf{x}-\mathbf{x}_{p}(t))$

and

$\displaystyle\zeta_{\sigma_{p}}(\mathbf{x}-\mathbf{x}_{p}(t))$
$\displaystyle=\frac{1}{\sigma_{p}^{3}(t)}\zeta\left(\frac{\|\mathbf{x}-\mathbf{x}_{p}(t)\|}{\sigma_{p}(t)}\right).$

Replacing the original filter $\zeta_{\sigma}$ with the variable-width filter
$\zeta_{\sigma_{p}}$ introduces an error from commuting filter and
differentiation operators in the LES equations. However, this error is second
order in the filter width Ghosal and Moin (1995), hence it is assumed to be
negligible or part of the discretization error of the VPM. Also, Eq. 7
introduces numerical issues by approximating $\overline{\boldsymbol{\omega}}$
through a field that is in general not divergence-free. Cottet Cottet and
Koumoutsakos (2000a) suggested that the divergence of
$\overline{\boldsymbol{\omega}}$ arises from unphysical small scales produced
by the computation, which could be properly addressed with subfilter-scale
diffusion. The divergence of $\overline{\boldsymbol{\omega}}$ will be treated
using the relaxation proposed by Pedrizzetti Pedrizzetti (1992), which, as
shown in Section VI, will be sufficient to attain both numerical stability and
physical accuracy with the LES formulation and subfilter-scale model developed
in this study.

The rotational part of the velocity field is calculated from the particles
using the Helmholtz decomposition
$\overline{\mathbf{u}}=\nabla\times\boldsymbol{\psi}$ where
$\boldsymbol{\psi}$ is some vector potential. This is done by analytically
solving the Poisson equation
$\nabla^{2}\boldsymbol{\psi}=-\overline{\boldsymbol{\omega}}$, which solution
is a regularized Biot-Savart law. Fig. 2 shows a vortex particle with unit
size and strength, the spreading of the vortex strength by the filter kernel,
and the resulting velocity field.

###  II.5 General VPM Governing Equations

We will now use vortex particles to discretize the LES vorticity equation, Eq.
5, and derive the equations governing the evolution of the Lagrangian
elements. For ease of notation, here on we denote the filtered velocity field
$\overline{\mathbf{u}}$ simply $\mathbf{u}$ and use
$\frac{\mathrm{D}}{\mathrm{D}t}$ and $\frac{\mathrm{d}}{\mathrm{d}t}$
interchangeably. Also, time dependence is no longer explicitly indicated, but
$\mathbf{x}_{p}$, $\boldsymbol{\Gamma}_{p}$, and $\sigma_{p}$ are time-
dependent variables.

Starting from the inviscid part of the LES-filtered vorticity equation,

$\displaystyle\frac{\text{d}}{\text{d}t}\overline{\boldsymbol{\omega}}$
$\displaystyle=\left(\overline{\boldsymbol{\omega}}\cdot\nabla\right)\mathbf{u}-\mathbf{E}_{\mathrm{adv}}-\mathbf{E}_{\mathrm{str}},$

we write the filter operator explicitly,

$\displaystyle\frac{\text{d}}{\text{d}t}\left(\int\limits_{-\infty}^{\infty}\boldsymbol{\omega}\left(\mathbf{y}\right)\zeta_{\sigma}(\mathbf{x}-\mathbf{y})\,\mathrm{d}\mathbf{y}\right)=$
$\displaystyle\qquad\qquad\qquad\left[\left(\int\limits_{-\infty}^{\infty}\boldsymbol{\omega}\left(\mathbf{y}\right)\zeta_{\sigma}(\mathbf{x}-\mathbf{y})\,\mathrm{d}\mathbf{y}\right)\cdot\nabla\right]\mathbf{u}\left(\mathbf{x}\right)$
$\displaystyle\qquad\qquad\qquad\qquad\qquad\qquad\qquad-\mathbf{E}_{\mathrm{adv}}\left(\mathbf{x}\right)-\mathbf{E}_{\mathrm{str}}\left(\mathbf{x}\right).$

Using the singular particle approximation,

$\displaystyle\boldsymbol{\omega}({\textbf{y}})\approx\sum\limits_{q}\boldsymbol{\Gamma}_{q}\delta({\textbf{y}}-{\textbf{x}}_{q}),$

both integrals collapse resulting in

$\displaystyle\frac{\mathrm{d}}{\mathrm{d}t}\left(\sum\limits_{q}\boldsymbol{\Gamma}_{q}\zeta_{\sigma_{q}}({\mathbf{x}}-{\mathbf{x}}_{q})\right)=$
(8)
$\displaystyle\left[\left(\sum\limits_{q}\boldsymbol{\Gamma}_{q}\zeta_{\sigma_{q}}({\mathbf{x}}-{\mathbf{x}}_{q})\right)\cdot\nabla\right]{\mathbf{u}}({\mathbf{x}})-\mathbf{E}_{\mathrm{adv}}({\mathbf{x}})-\mathbf{E}_{\mathrm{str}}({\mathbf{x}}).$

With the help of Appendix A, the left-hand side unfolds as

$\displaystyle\frac{\mathrm{d}}{\mathrm{d}t}\left(\sum\limits_{q}\boldsymbol{\Gamma}_{q}\zeta_{\sigma_{q}}({\mathbf{x}}-{\mathbf{x}}_{q})\right)$
$\displaystyle=\sum\limits_{q}\left[\frac{\mathrm{d}\boldsymbol{\Gamma}_{q}}{\mathrm{d}t}\zeta_{\sigma_{q}}({\mathbf{x}}-{\mathbf{x}}_{q})+\boldsymbol{\Gamma}_{q}\frac{\mathrm{d}}{\mathrm{d}t}\left(\zeta_{\sigma_{q}}({\mathbf{x}}-{\mathbf{x}}_{q})\right)\right]$
$\displaystyle=\sum\limits_{q}\left\\{\frac{\mathrm{d}\boldsymbol{\Gamma}_{q}}{\mathrm{d}t}\zeta_{\sigma_{q}}({\mathbf{x}}-{\mathbf{x}}_{q})+\boldsymbol{\Gamma}_{q}\frac{\partial\zeta_{\sigma_{q}}}{\partial
t}({\mathbf{x}}-{\mathbf{x}}_{q})\right.$
$\displaystyle\qquad\qquad\qquad\left.+\boldsymbol{\Gamma}_{q}\left[\left(\mathbf{u}({\mathbf{x}})-\mathbf{u}({\mathbf{x}}_{q})\right)\cdot\nabla\right]\zeta_{\sigma_{q}}({\mathbf{x}}-{\mathbf{x}}_{q})\right\\}$

while vortex stretching is rearranged as

$\displaystyle\left[\left(\sum\limits_{q}\boldsymbol{\Gamma}_{q}\zeta_{\sigma_{q}}({\mathbf{x}}-{\mathbf{x}}_{q})\right)\cdot\nabla\right]{\mathbf{u}}({\mathbf{x}})=$
$\displaystyle\qquad\qquad\qquad\qquad\qquad\sum\limits_{q}\zeta_{\sigma_{q}}({\mathbf{x}}-{\mathbf{x}}_{q})\left(\boldsymbol{\Gamma}_{q}\cdot\nabla\right){\mathbf{u}}({\mathbf{x}}).$

Evaluating at the position of the $p$-th particle,
$\mathbf{x}=\mathbf{x}_{p}$, and pulling the $p$-th term out of each sum,

$\displaystyle\frac{\mathrm{d}}{\mathrm{d}t}\left(\sum\limits_{q}\boldsymbol{\Gamma}_{q}\zeta_{\sigma_{q}}({\mathbf{x}}_{p}-{\mathbf{x}}_{q})\right)$
$\displaystyle=\frac{\mathrm{d}\boldsymbol{\Gamma}_{p}}{\mathrm{d}t}\zeta_{\sigma_{p}}({\mathbf{0}})+\boldsymbol{\Gamma}_{p}\frac{\partial\zeta_{\sigma_{p}}}{\partial
t}({\mathbf{0}})$ $\displaystyle\quad+\sum\limits_{q\neq
p}\left\\{\frac{\mathrm{d}\boldsymbol{\Gamma}_{q}}{\mathrm{d}t}\zeta_{\sigma_{q}}({\mathbf{x}}_{p}-{\mathbf{x}}_{q})+\boldsymbol{\Gamma}_{q}\frac{\partial\zeta_{\sigma_{q}}}{\partial
t}({\mathbf{x}}_{p}-{\mathbf{x}}_{q})\right.$
$\displaystyle\qquad\qquad\qquad\left.+\,\boldsymbol{\Gamma}_{q}\left[\left(\mathbf{u}({\mathbf{x}}_{p})-\mathbf{u}({\mathbf{x}}_{q})\right)\cdot\nabla\right]\zeta_{\sigma_{q}}({\mathbf{x}}_{p}-{\mathbf{x}}_{q})\right\\}$

and

$\displaystyle\left[\left(\sum\limits_{q}\boldsymbol{\Gamma}_{q}\zeta_{\sigma_{q}}({\mathbf{x}}_{p}-{\mathbf{x}}_{q})\right)\cdot\nabla\right]{\mathbf{u}}({\mathbf{x}}_{p})$
$\displaystyle=\zeta_{\sigma_{p}}({\mathbf{0}})\left(\boldsymbol{\Gamma}_{p}\cdot\nabla\right){\mathbf{u}}({\mathbf{x}}_{p})+\sum\limits_{q\neq
p}\zeta_{\sigma_{q}}({\mathbf{x}}_{p}-{\mathbf{x}}_{q})\left(\boldsymbol{\Gamma}_{q}\cdot\nabla\right){\mathbf{u}}({\mathbf{x}}_{p}).$

Thus, Eq. 8 evaluated at $\mathbf{x}=\mathbf{x}_{p}$ becomes

$\displaystyle\frac{\mathrm{d}\boldsymbol{\Gamma}_{p}}{\mathrm{d}t}\zeta_{\sigma_{p}}({\mathbf{0}})$
$\displaystyle+\boldsymbol{\Gamma}_{p}\frac{\partial\zeta_{\sigma_{p}}}{\partial
t}({\mathbf{0}})+\mathbf{M}_{p}^{0}$ (9)
$\displaystyle=\zeta_{\sigma_{p}}({\mathbf{0}})\left(\boldsymbol{\Gamma}_{p}\cdot\nabla\right)\mathbf{u}({\mathbf{x}}_{p})$
$\displaystyle\quad+(\mathbf{M}_{p}^{1}+\mathbf{M}_{p}^{2})-\left(\mathbf{E}_{\mathrm{adv}}({\mathbf{x}}_{p})+\mathbf{E}_{\mathrm{str}}({\mathbf{x}}_{p})\right),$

where

$\displaystyle\bullet\quad$
$\displaystyle\mathbf{M}_{p}^{0}\equiv\sum\limits_{q\neq
p}\left(\frac{\mathrm{d}\boldsymbol{\Gamma}_{q}}{\mathrm{d}t}\zeta_{\sigma_{q}}({\mathbf{x}}_{p}-{\mathbf{x}}_{q})+\boldsymbol{\Gamma}_{q}\frac{\partial\zeta_{\sigma_{q}}}{\partial
t}({\mathbf{x}}_{p}-{\mathbf{x}}_{q})\right)$ (10) $\displaystyle\bullet\quad$
$\displaystyle\mathbf{M}_{p}^{1}\equiv-\sum\limits_{q\neq
p}\boldsymbol{\Gamma}_{q}\left(\mathbf{u}({\mathbf{x}}_{p})-\mathbf{u}({\mathbf{x}}_{q})\right)\cdot\nabla\zeta_{\sigma_{q}}({\mathbf{x}}_{p}-{\mathbf{x}}_{q})$
(11) $\displaystyle\bullet\quad$
$\displaystyle\mathbf{M}_{p}^{2}\equiv\sum\limits_{q\neq
p}\zeta_{\sigma_{q}}({\mathbf{x}}_{p}-{\mathbf{x}}_{q})\left(\boldsymbol{\Gamma}_{q}\cdot\nabla\right){\mathbf{u}}({\mathbf{x}}_{p})$
(12)

These $\mathbf{M}$-terms pose an interdependence between the $p$-th particle
and neighboring particles, which arises from filtering the vorticity equation.
In this study we will neglect this interdependence as explained in Section
III. However, we will see in Section V that $\mathbf{M}_{p}^{1}$ and
$\mathbf{M}_{p}^{2}$ are closely related to advection and vortex stretching at
the subfilter scale.

Recalling that $\zeta_{\sigma}$ is defined as
$\zeta_{\sigma}(\mathbf{x})=\frac{1}{\sigma^{3}}\zeta(\frac{\|\mathbf{x}\|}{\sigma})$,
its time derivative is calculated as

$\displaystyle\frac{\partial\zeta_{\sigma}}{\partial t}({\mathbf{x}})$
$\displaystyle=-3\frac{1}{\sigma^{4}}\frac{\partial\sigma}{\partial
t}\zeta\left(\frac{\|\mathbf{x}\|}{\sigma}\right)+\frac{1}{\sigma^{3}}\frac{\partial}{\partial
t}\left(\zeta\left(\frac{\|\mathbf{x}\|}{\sigma}\right)\right)$
$\displaystyle=-3\frac{1}{\sigma}\frac{\partial\sigma}{\partial
t}\zeta_{\sigma}(\mathbf{x})-\frac{1}{\sigma^{3}}\frac{\partial\zeta}{\partial
r}\left(\frac{\|\mathbf{x}\|}{\sigma}\right)\frac{\|\mathbf{x}\|}{\sigma^{2}}\frac{\partial\sigma}{\partial
t}.$

Assuming that $\zeta(r)$ reaches a maximum at $r=0$ (i.e.,
$\frac{\partial\zeta}{\partial r}(0)=0$) and evaluating at
$\mathbf{x}=\mathbf{0}$, we get

$\displaystyle\frac{\partial\zeta_{\sigma}}{\partial t}({\mathbf{0}})$
$\displaystyle=-3\frac{1}{\sigma}\frac{\partial\sigma}{\partial
t}\zeta_{\sigma}(\mathbf{0}).$ (13)

Since viscous effects have been set aside through operator splitting,
$\frac{\partial\sigma}{\partial t}$ in this derivation only accounts for
inviscid effects. For clarity, this means that core spreading due to viscous
diffusion must not be included in Eq. 13.

Finally, substituting Eq. 13 into Eq. 9 and assuming
$\zeta_{\sigma_{p}}({\mathbf{0}})\neq 0$, we arrive to the equation governing
the evolution of vortex strength,

$\displaystyle\frac{\mathrm{d}\boldsymbol{\Gamma}_{p}}{\mathrm{d}t}=\left(\boldsymbol{\Gamma}_{p}\cdot\nabla\right)\mathbf{u}({\mathbf{x}}_{p})$
$\displaystyle+3\boldsymbol{\Gamma}_{p}\frac{1}{\sigma_{p}}\frac{\partial\sigma_{p}}{\partial
t}$ (14)
$\displaystyle+\frac{1}{\zeta_{\sigma_{p}}({\mathbf{0}})}(-\mathbf{M}_{p}^{0}+\mathbf{M}_{p}^{1}+\mathbf{M}_{p}^{2})$
$\displaystyle-\frac{1}{\zeta_{\sigma_{p}}({\mathbf{0}})}\left(\mathbf{E}_{\mathrm{adv}}({\mathbf{x}}_{p})+\mathbf{E}_{\mathrm{str}}({\mathbf{x}}_{p})\right).$

Thus, particle convection as in

$\displaystyle\frac{\text{d}}{\text{d}t}{\mathbf{x}}_{p}={\mathbf{u}}({\mathbf{x}}_{p}),$
(15)

strength evolution as in Eq. 14, some expression for
$\frac{\partial\sigma_{p}}{\partial t}$, and the viscous diffusion equation
make up the general governing equations of the vortex particle method that
solve the LES-filtered Navier-Stokes vorticity equation.

We see in Eq. 14 that the evolution of the vortex strength is dictated by
vortex stretching (first term), inviscid expansion/contraction of the particle
size $\sigma_{p}$ (second term), a dependence on other particles through the
$\mathbf{M}$ terms, and SFS contributions through the $\mathbf{E}$ terms. One
obtains different formulations of the VPM depending on how
$\frac{\partial\sigma}{\partial t}$, the $\mathbf{M}$ terms, and the SFS
$\mathbf{E}$ terms are handled. In particular, we will see in Section III that
the classic VPM is equivalent to assuming†††The only instance that the classic
VPM uses $\frac{\partial\sigma}{\partial t}\neq 0$ is in solving the viscous
diffusion equation through the core spreading scheme, but even in this case
the effects of $\frac{\partial\sigma}{\partial t}$ on vortex strength are
neglected. $\frac{\partial\sigma}{\partial t}=0$ while neglecting all
$\mathbf{M}$ terms. In Section IV we propose a new formulation that uses
$\frac{\partial\sigma}{\partial t}\neq 0$ to reinforce conservations laws in
spherical elements. In Section V we discuss existing SFS models for meshless
schemes and develop a new anisotropic dynamic model of SFS vortex stretching,
$\mathbf{E}_{\mathrm{str}}$.

###  II.6 Numerical Schemes

In FLOWVPM—the solver used for this study—vortex stretching is solved in the
transposed scheme Winckelmans (1989); Winckelmans and Leonard (1993) and the
divergence of the vorticity field is treated through the relaxation scheme
developed by Pedrizzeti Pedrizzetti (1992). The time integration of the
governing equations is done through a low-storage third-order Runge-Kutta
scheme Williamson (1980). A Gaussian kernel is used as the LES filter
$\zeta_{\sigma}$ (or VPM radial basis function). The fast multipole method
Greengard (1987); Cheng, Greengard, and Rokhlin (1999) (FMM) is used for the
computation of the regularized Biot-Savart law, approximating the velocity
field and vortex stretching through spherical harmonics with computational
complexity $\mathcal{O}(N)$, where $N$ is the number of particles. The FMM
computation of vortex stretching is performed through an efficient complex-
step derivative approximation Alvarez and Ning (2020), implemented in a
modified version of the open-source, parallelized code ExaFMMYokota and Barba
(2011); Wang, Yokota, and Barba (2021). FLOWVPM is implemented in the Julia
language Bezanson  _et al._ (2017), which is a modern, high-level, dynamic
programming language for high-performance computing.

##  III Classic Vortex Particle Method

Before further exploring the general equation governing vortex stretching, Eq.
14, we pause to analyze the equation that has been used extensively in the
literature throughout the years:

$\displaystyle\frac{\mathrm{d}\boldsymbol{\Gamma}_{p}}{\mathrm{d}t}=\left(\boldsymbol{\Gamma}_{p}\cdot\nabla\right)\mathbf{u}(\mathbf{x}_{p}).$
(16)

We denote this equation as the Classic VPM.

Multiple variations of Eq. 16 have been used over the years. For instance,
Gharakhani Gharakhani (1997) introduced a new term accounting for the
divergence of the approximated vorticity field, Winckelmans and Leonard
Winckelmans and Leonard (1993) replaced the differential operator in vortex
stretching with its transpose, and Mansfield et al. Mansfield, Knio, and
Meneveau (1998, 1999) and Cottet Cottet (1996) reintroduced the SFS
contributions. However, all these variations and the classic equation have in
common that they all neglect $\frac{\partial\sigma_{p}}{\partial t}$ and are
free of the interdependence posed by the $\mathbf{M}$-terms in Eq. 14.

We now delineate the assumption underlying the classic VPM that will justify
neglecting the $\mathbf{M}$-terms, termed localized-vorticity assumption, as
follows.

Given a vorticity field $\boldsymbol{\omega}$ that is compact in a small
volume $\mathrm{Vol}$, the average vorticity
$\left<\boldsymbol{\omega}\right>$ is calculated as

$\displaystyle\left<\boldsymbol{\omega}\right>=\frac{\int\limits_{\mathrm{Vol}}\boldsymbol{\omega}(\mathbf{y})\,\mathrm{d}\mathbf{y}}{\mathrm{Vol}}.$

The field can be approximated through a radially-symmetric field
$\tilde{\boldsymbol{\omega}}$ defined as

$\displaystyle\tilde{\boldsymbol{\omega}}(\mathbf{x})\equiv\left<\boldsymbol{\omega}\right>\mathrm{Vol}\zeta_{\sigma}(\mathbf{x}-\mathbf{x}_{0}),$

where $\zeta_{\sigma}$ is a radial basis function of spread $\sigma$ and
center $\mathbf{x}_{0}$ approximating the vorticity distribution of the
original field $\boldsymbol{\omega}$. If $\zeta_{\sigma}$ is normalized such
as to have a volume integral of unity, $\tilde{\boldsymbol{\omega}}$
approximates $\boldsymbol{\omega}$ in an average sense since

$\displaystyle\int\limits_{-\infty}^{\infty}\tilde{\boldsymbol{\omega}}(\mathbf{y})\,\mathrm{d}\mathbf{y}$
$\displaystyle=\int\limits_{\mathrm{Vol}}\boldsymbol{\omega}(\mathbf{y})\,\mathrm{d}\mathbf{y}.$

Defining $\boldsymbol{\Gamma}=\left<\boldsymbol{\omega}\right>\mathrm{Vol}$,
the localized-vorticity field $\boldsymbol{\omega}$ can be approximated with a
single particle as

$\displaystyle\boldsymbol{\omega}\left(\mathbf{x}\right)\approx\boldsymbol{\Gamma}\zeta_{\sigma}(\mathbf{x}-\mathbf{x}_{0}).$

Replacing this in the inviscid part of the vorticity equation,

$\displaystyle\frac{\text{d}}{\text{d}t}\boldsymbol{\omega}=(\boldsymbol{\omega}\cdot\nabla){\mathbf{u}},$

evaluating at $\mathbf{x}=\mathbf{x}_{0}$, and following similar steps as in
Section II.5, we get

$\displaystyle\frac{\mathrm{d}\boldsymbol{\Gamma}}{\mathrm{d}t}=\left(\boldsymbol{\Gamma}\cdot\nabla\right){\mathbf{u}}(\mathbf{x}_{0})+3\boldsymbol{\Gamma}\frac{1}{\sigma}\frac{\partial\sigma}{\partial
t}.$ (17)

This is the evolution equation of the vortex strength approximating the field
of localized vorticity, discretized with only one particle. In contrast, in
Section II.5 we showed that a more general vorticity field that is discretized
using multiple particles leads to a governing equation that is slightly
different:

$\displaystyle\frac{\mathrm{d}\boldsymbol{\Gamma}_{p}}{\mathrm{d}t}=\left(\boldsymbol{\Gamma}_{p}\cdot\nabla\right)\mathbf{u}(\mathbf{x}_{p})$
$\displaystyle+3\boldsymbol{\Gamma}_{p}\frac{1}{\sigma_{p}}\frac{\partial\sigma_{p}}{\partial
t}$ (18)
$\displaystyle+\frac{1}{\zeta_{\sigma_{p}}({\mathbf{0}})}(-\mathbf{M}_{p}^{0}+\mathbf{M}_{p}^{1}+\mathbf{M}_{p}^{2}),$

after ignoring the SFS contributions. Here, the $\mathbf{M}$-terms, as defined
in Eqs. 10, 11 and 12, pose a dependence on neighboring particles that arises
from having filtered the vorticity equation.

Discarding these $\mathbf{M}$-terms is equivalent to assuming that the
vorticity field can be approximated by the superposition of blobs of fluid
with compact vorticity that evolve somewhat independently from each
other.‡‡‡The interactions neglected here are only the ones pertaining to the
evolution of vortex strength, while other interactions are still accounted for
in the convection and viscous diffusion of the particles. This assumption,
which we call localized-vorticity assumption, reduces Eq. 18 into Eq. 17 and
is valid as long as there is no significant particle overlap. In both
derivations of the classic VPM by Winckelmans and Leonard Winckelmans and
Leonard (1993) and Cottet and Koumutsakos Cottet and Koumoutsakos (2000b),
these $\mathbf{M}$-terms are not present due to the construction of the
method: In the classic derivation, the unfiltered vorticity equation (Eq. 2)
is discretized with singular particles and only the velocity field is filtered
to obtain a regularized field, while in our derivation we discretized the LES-
filtered vorticity equation (Eq. 5) which leads to Eq. 18. In order to bring
both approaches into agreement, in Section IV we will use the localized-
vorticity assumption in the derivation of our reformulated VPM to neglect
these $\mathbf{M}$-terms.

Interestingly, discarding the $\mathbf{M}$-terms resembles the LES
decomposition approach of truncated basis functions Meneveau and Katz (2000).
In such an approach, the flow field is expanded using orthonormal basis
functions. The summation of bases is then truncated to define the large-scale
field, and the discarded modes represent the range of subfilter scales. In the
localized-vorticity assumption, the sum over all the particles is truncated
after the leading term, $p$, while neglecting the contributions of neighboring
particles. Hence, the localized-vorticity assumption can be regarded as a
secondary LES filter with the neglected $\mathbf{M}$-terms becoming part of
the subfilter-scale contributions.

While it is justifiable to neglect the $\mathbf{M}$-terms through the
localized-vorticity assumption, it is unclear to us what the basis is for the
classic VPM to assume $\frac{\partial\sigma_{p}}{\partial t}=0$. In fact, we
hypothesize that this last assumption is the cause of the numerical
instabilities that pervade the classic VPM.

The classic VPM simply regards $\sigma$ as a numerical parameter with no
physical significance. However, Leonard Leonard (1980) suggested that $\sigma$
should change according to conservation of mass, and Nakanishi, Ojima, and
Maremoto Nakanishi and Kamemoto (1993); Ojima and Kamemoto (2000) suggested
that $\sigma$ should change according to Kelvin’s theorem, which was more
recently implemented by Kornev et al. Kornev, Denev, and Samarbakhsh (2020)
Even though these authors let $\sigma$ evolve in time, they did not include
these effects back in the equation that governs $\boldsymbol{\Gamma}$,
effectively assuming $\frac{\partial\sigma_{p}}{\partial t}=0$ in the
evolution of $\boldsymbol{\Gamma}$. In the following chapter we propose a
formulation that uses $\frac{\partial\sigma_{p}}{\partial t}\neq 0$ in the
governing equation of $\boldsymbol{\Gamma}$, while letting $\sigma$ change as
to reinforce conservation of both mass and angular momentum, which will be
shown to lead to remarkable numerical stability.

##  IV Reformulated Vortex Particle Method

The general governing equations derived in Section II.5 call for an expression
for $\frac{\partial\sigma_{p}}{\partial t}$. While the classic VPM simply
assumes $\frac{\partial\sigma_{p}}{\partial t}=0$, we will now go back to
first principles to find some plausible expressions for
$\frac{\partial\sigma_{p}}{\partial t}$. As a preamble, in Section IV.1 we put
forth two physical implications that distill from the vorticity equation,
which will then guide our search for possible candidate formulations of
$\frac{\partial\sigma_{p}}{\partial t}$ in Section IV.2. The process for
constructing candidate formulation is then generalized in Section IV.3, and
the entire space of possible formulation is analyzed in Section IV.4. This
will then lead us to one formulation, termed reformulated VPM and summarized
in Section IV.5, which reinforces conservation of both mass and angular
momentum by reshaping the vortex elements as they are subject to vortex
stretching.

###  IV.1 Physical Implications of the Vorticity Equation

In Section II.1 we showed that the linear-momentum Navier-Stokes equation can
be transformed into an expression that only depends on vorticity, namely Eq.
2, here repeated:

$\displaystyle\frac{\text{d}}{\text{d}t}\boldsymbol{\omega}=(\boldsymbol{\omega}\cdot\nabla)\mathbf{u}+\nu\nabla^{2}\boldsymbol{\omega}.$
(2)

We now point our attention to two laws that follow from Eq. 2: conservation of
mass and angular momentum. The discussion that follows is inspired by the
writing of P. A. Davidson Davidson (2001, 2019).

####  IV.1.1 Conservation of Angular Momentum

Consider a spherical differential fluid element carrying a mean vorticity
$\boldsymbol{\omega}$ with moment of inertia $I$. Due to
$\boldsymbol{\omega}$, the element is then rotating at an angular velocity of
$\boldsymbol{\omega}/2$ and its angular momentum $\mathbf{L}$ is calculated as
(Batchelor, 1967, p. 82)

$\displaystyle\mathbf{L}=\frac{I\boldsymbol{\omega}}{2}.$ (19)

Given that the element is spherical (and before strain distorts the element),
the pressure field exerts no torque on the element and the only torque-
producing forces are due to viscous effects, namely
$\boldsymbol{\tau}_{\mathrm{viscous}}$. The change of angular momentum is then
calculated as

$\displaystyle\frac{\mathrm{d}}{\mathrm{d}t}\mathbf{L}=\boldsymbol{\tau}_{\mathrm{viscous}}.$
(20)

Replacing Eq. 19 in Eq. 20, we arrive to the following expression of the total
derivative of vorticity§§§More generally, this equation should be expressed
with an inertia tensor, but for a spherical body this tensor reduces to a
scalar due to the symmetry of the body.,

$\displaystyle\frac{\mathrm{d}}{\mathrm{d}t}\boldsymbol{\omega}=-\frac{\boldsymbol{\omega}}{I}\frac{\mathrm{d}}{\mathrm{d}t}I+\boldsymbol{\tau}_{\mathrm{viscous}}^{*},$
(21)

where
$\boldsymbol{\tau}_{\mathrm{viscous}}^{*}=2\boldsymbol{\tau}_{\mathrm{viscous}}/I$.
When viscous effects are ignored, Eq. 21 implies that, in order to conserve
angular momentum, the angular velocity (or vorticity) must decrease whenever
the moment of inertia increases, and vice versa.

Eqs. 2 and 21 suggest that the vorticity Navier-Stokes equation is simply an
expression of the conservation of angular momentum in a spherical fluid
elementChatwin (1973). Furthermore, the first term in the right hand side of
Eq. 2 accounts for the change of moment of inertia as

$\displaystyle(\boldsymbol{\omega}\cdot\nabla)\mathbf{u}=-\frac{\boldsymbol{\omega}}{I}\frac{\mathrm{d}}{\mathrm{d}t}I,$
(22)

which leads to an increase/decrease of vorticity to conserve momentum as the
moment of inertia decreases/increases. Hence,
$(\boldsymbol{\omega}\cdot\nabla)\mathbf{u}$ is referred to as vortex
stretching as it accounts for the deformation exerted by the velocity field on
the fluid element, intensifying the vorticity in the direction that the
element is stretched, as illustrated in Fig. 3.

Figure 3: Stretching of a spherical fluid element and increase in vorticity
(or angular velocity).

####  IV.1.2 Conservation of Mass

Another physical implication of Eq. 2 is the existence of vortex lines and
vortex tubes, described as follows. We define a short material line
$\boldsymbol{\ell}$ that moves with the fluid, with start and end points
$\mathbf{x}$ and $\mathbf{x}+\boldsymbol{\ell}$, respectively. If
$\boldsymbol{\ell}$ is infinitesimally small, the evolution of
$\boldsymbol{\ell}$ is given by

$\displaystyle\frac{\mathrm{d}}{\mathrm{d}t}\boldsymbol{\ell}=\mathbf{u}\left(\mathbf{x}+\boldsymbol{\ell}\right)-\mathbf{u}\left(\mathbf{x}\right)=\left(\boldsymbol{\ell}\cdot\nabla\right)\mathbf{u}\left(\mathbf{x}\right).$
(23)

Defining $\boldsymbol{\ell}$ as being tangent to the local vorticity, and
comparing Eq. 23 to the inviscid part of Eq. 2, we see that vorticity evolves
as the material line, being identically stretched and reoriented by the
velocity field. This led Helmholtz to conclude the existence of lines of
vorticity, termed vortex lines, that move with the fluid. We then define a
vortex tube as the surface formed by all the vortex lines passing through a
closed curve, depicted in Fig. 4.

The fact that vortex lines move with the fluid implies that the volume of a
vortex tube segment must be conserved in incompressible flow. For instance,
let us discretize and approximate a vortex tube through Lagrangian cylindrical
elements of length $\ell$ and cross section $r$. As the length of the element
is stretched, the cross section $r$ must shrink to conserve the same volume
$\mathrm{Vol}=\pi r^{2}\ell$. This is,

$\displaystyle\frac{\mathrm{d}}{\mathrm{d}t}\mathrm{Vol}=0,$

leading to

$\displaystyle\frac{\mathrm{d}}{\mathrm{d}t}r=-\frac{r}{2\ell}\frac{\mathrm{d}}{\mathrm{d}t}\ell.$
(24)

Figure 4: Vortex tube formed by the vortex lines passing through a closed
curve.

####  IV.1.3 Implications for the VPM

The aforementioned laws of conservation point our attention to two major
pitfalls of the classic VPM, if particles are considered as material elements.
First, Eq. 22 establishes that vortex stretching must lead to a change in the
shape of the fluid element. Otherwise, if the vorticity (angular velocity) has
changed due to vortex stretching while the shape of the element (moment of
inertia) is kept constant, the conservation of angular momentum has been
locally violated. This is exactly the case when the classic VPM assumes a
constant core size, $\frac{\mathrm{d}}{\mathrm{d}t}\sigma=0$.

Second, in an incompressible tube of vorticity, Eqs. 23 and 2 establish that
an increase in vorticity due to vortex stretching is associated to a
lengthening of the tube, which must cause the cross section of the tube to
shrink according to Eq. 24. Otherwise, the tube has kept a constant cross
section while being stretched, thus locally violating conservation of mass.
Again, this is exactly the case when the classic VPM assumes
$\frac{\mathrm{d}}{\mathrm{d}t}\sigma=0$.

In summary, we believe that the classic VPM violates conservation of both
angular momentum and mass when it assumes
$\frac{\mathrm{d}}{\mathrm{d}t}\sigma=0$. We also hypothesize that this is the
reason for its poor numerical stability. In the following sections we will
consider some possible formulations that use
$\frac{\mathrm{d}}{\mathrm{d}t}\sigma$ as a degree of freedom to locally
reinforce these conservation laws.

###  IV.2 Candidate Formulations

Starting from the general equation governing vortex strength (Eq. 14),
neglecting SFS terms, and using the localized-vorticity assumption, we obtain

$\displaystyle\frac{\mathrm{d}}{\mathrm{d}t}\boldsymbol{\Gamma}_{p}=\left(\boldsymbol{\Gamma}_{p}\cdot\nabla\right)\mathbf{u}(\mathbf{x}_{p})+3\boldsymbol{\Gamma}_{p}\frac{1}{\sigma_{p}}\frac{\partial\sigma_{p}}{\partial
t}.$ (25)

We will now explore multiple combinations of element shapes and conservation
laws to derive possible expressions for $\frac{\partial\sigma_{p}}{\partial
t}$.

####  IV.2.1 Sphere with Conservation of Momentum

Given that the moment of inertia of a solid sphere is $I=\frac{2}{5}mr^{2}$,
where $m=\frac{4}{3}\rho\pi r^{3}$ is the mass of the sphere and $r$ its
radius, the inertial term in Eq. 22 becomes

$\displaystyle\frac{1}{I}\frac{\mathrm{d}}{\mathrm{d}t}I=\frac{1}{r^{5}}\frac{\mathrm{d}}{\mathrm{d}t}\left(r^{5}\right)=\frac{5}{r}\frac{\mathrm{d}}{\mathrm{d}t}r.$

Replacing this into Eq. 22 and taking the dot product with
$\boldsymbol{\omega}$,

$\displaystyle-\left(\boldsymbol{\omega}\cdot\boldsymbol{\omega}\right)\frac{5}{r}\frac{\mathrm{d}}{\mathrm{d}t}r=\left[(\boldsymbol{\omega}\cdot\nabla)\mathbf{u}\right]\cdot\boldsymbol{\omega},$

results in

$\displaystyle\frac{1}{r}\frac{\mathrm{d}}{\mathrm{d}t}r$
$\displaystyle=-\frac{1}{5}\frac{1}{\|\boldsymbol{\omega}\|^{2}}\left[(\boldsymbol{\omega}\cdot\nabla)\mathbf{u}\right]\cdot\boldsymbol{\omega}$
$\displaystyle=-\frac{1}{5}\left[(\hat{\boldsymbol{\omega}}\cdot\nabla)\mathbf{u}\right]\cdot\hat{\boldsymbol{\omega}}.$

When $\boldsymbol{\omega}$ is filtered and discretized with vortex particles,
and the localized-vorticity assumption is applied, the spherical element
becomes the vortex particle itself. Vortex strengths are roughly aligned with
the vorticity field¶¶¶This is especially true when using Pedrizzetti’s
relaxation scheme.Pedrizzetti (1992), meaning
$\hat{\boldsymbol{\Gamma}}_{p}\approx\hat{\boldsymbol{\omega}}\left(\mathbf{x}_{p}\right)$,
and the size of the particle can be expressed as $\sigma_{p}=\alpha r$ with
$\alpha$ some scaling factor, obtaining

$\displaystyle\frac{1}{\sigma_{p}}\frac{\mathrm{d}}{\mathrm{d}t}\sigma_{p}$
$\displaystyle=-\frac{1}{5}\left[(\hat{\boldsymbol{\Gamma}}_{p}\cdot\nabla)\mathbf{u}\left(\mathbf{x}_{p}\right)\right]\cdot\hat{\boldsymbol{\Gamma}}_{p},$

or, equivalently

$\displaystyle\frac{\mathrm{d}}{\mathrm{d}t}\sigma_{p}$
$\displaystyle=-\frac{1}{5}\frac{\sigma_{p}}{\|\boldsymbol{\Gamma}_{p}\|}\left[(\boldsymbol{\Gamma}_{p}\cdot\nabla)\mathbf{u}\left(\mathbf{x}_{p}\right)\right]\cdot\hat{\boldsymbol{\Gamma}}_{p}.$
(26)

Noticing that, in a Lagrangian scheme, $\sigma_{p}$ is only a function of time
(i.e.,
$\frac{\mathrm{d}}{\mathrm{d}t}\sigma_{p}=\frac{\partial\sigma_{p}}{\partial
t}$), we replace Eq. 26 into Eq. 25 to obtain the equation governing
$\boldsymbol{\Gamma}_{p}$:

$\displaystyle\frac{\mathrm{d}}{\mathrm{d}t}\boldsymbol{\Gamma}_{p}=\left(\boldsymbol{\Gamma}_{p}\cdot\nabla\right)\mathbf{u}(\mathbf{x}_{p})-\frac{3}{5}\left\\{\left[(\boldsymbol{\Gamma}_{p}\cdot\nabla)\mathbf{u}\left(\mathbf{x}_{p}\right)\right]\cdot\hat{\boldsymbol{\Gamma}}_{p}\right\\}\hat{\boldsymbol{\Gamma}}_{p}.$
(27)

Thus, Eqs. 26 and 27 are the governing equations required for the particles to
preserve angular momentum.

####  IV.2.2 Tube with Conservation of Mass

Assuming that the field $\boldsymbol{\omega}$ encompasses a tube of vorticity
of radius $r$, we now use conservation of mass to derive an expression for
$\frac{\mathrm{d}}{\mathrm{d}t}\sigma$. When $\boldsymbol{\omega}$ is
discretized through cylindrical vortex elements of length $\ell$ and radius
$\sigma_{t}=\alpha_{t}r$, with $\alpha_{t}$ some scaling factor, Eq. 24
becomes

$\displaystyle\frac{1}{\sigma_{t}}\frac{\mathrm{d}}{\mathrm{d}t}\sigma_{t}=-\frac{1}{2\ell}\frac{\mathrm{d}}{\mathrm{d}t}\ell.$
(28)

The length $\ell$ of the cylindrical element relates to the evolution of its
vorticity $\boldsymbol{\omega_{t}}$ through vortex stretching as

$\displaystyle\frac{1}{\ell}\frac{\mathrm{d}}{\mathrm{d}t}\ell=\frac{1}{\|\boldsymbol{\omega_{t}}\|}\frac{\mathrm{d}\boldsymbol{\omega_{t}}}{\mathrm{d}t}\cdot\hat{\boldsymbol{\omega_{t}}}.$

Thus, Ojima and Maremoto Ojima and Kamemoto (2000) proposed that this
lengthening of $\ell$ and contraction of $\sigma_{t}$ caused by vortex
stretching in a cylindrical element translates into the following relationship
in order for a vortex particle to conserve mass:

$\displaystyle\frac{1}{\sigma_{p}}\frac{\mathrm{d}}{\mathrm{d}t}\sigma_{p}=-\frac{1}{2\|\boldsymbol{\Gamma}_{p}\|}\frac{\mathrm{d}\boldsymbol{\Gamma}_{p}}{\mathrm{d}t}\cdot\hat{\boldsymbol{\Gamma}}_{p}.$
(29)

Now, using Eq. 29 in the vortex strength equation, Eq. 25, we obtain

$\displaystyle\frac{\mathrm{d}}{\mathrm{d}t}\boldsymbol{\Gamma}_{p}=\left(\boldsymbol{\Gamma}_{p}\cdot\nabla\right)\mathbf{u}(\mathbf{x}_{p})-\frac{3}{2}\left(\frac{\mathrm{d}\boldsymbol{\Gamma}_{p}}{\mathrm{d}t}\cdot\hat{\boldsymbol{\Gamma}}_{p}\right)\hat{\boldsymbol{\Gamma}}_{p},$
(30)

and from its dot product with $\hat{\boldsymbol{\Gamma}}_{p}$ we get

$\displaystyle\frac{\mathrm{d}\boldsymbol{\Gamma}_{p}}{\mathrm{d}t}\cdot\hat{\boldsymbol{\Gamma}}_{p}=\frac{2}{5}\left[\left(\boldsymbol{\Gamma}_{p}\cdot\nabla\right)\mathbf{u}(\mathbf{x}_{p})\right]\cdot\hat{\boldsymbol{\Gamma}}_{p}.$

Finally, substituting this back into Eqs. 29 and 30, we obtain the governing
equations required for the particles to preserve mass in a tube of vorticity:

$\displaystyle\frac{\mathrm{d}}{\mathrm{d}t}\sigma_{p}$
$\displaystyle=-\frac{1}{5}\frac{\sigma_{p}}{\|\boldsymbol{\Gamma}_{p}\|}\left[\left(\boldsymbol{\Gamma}_{p}\cdot\nabla\right)\mathbf{u}(\mathbf{x}_{p})\right]\cdot\hat{\boldsymbol{\Gamma}}_{p}$
$\displaystyle\frac{\mathrm{d}}{\mathrm{d}t}\boldsymbol{\Gamma}_{p}$
$\displaystyle=\left(\boldsymbol{\Gamma}_{p}\cdot\nabla\right)\mathbf{u}(\mathbf{x}_{p})-\frac{3}{5}\left\\{\left[\left(\boldsymbol{\Gamma}_{p}\cdot\nabla\right)\mathbf{u}(\mathbf{x}_{p})\right]\cdot\hat{\boldsymbol{\Gamma}}_{p}\right\\}\hat{\boldsymbol{\Gamma}}_{p}$

Surprisingly, we have arrived to the same governing equations of the momentum-
conserving sphere formulation∥∥∥With the caveat that this is so only because
we have neglected the SFS term. For this reason, in Section IV.3 we will
revert back to Eq. 29., Eqs. 26 and 27. This stems from the fact that both
conservation laws are physical consequences of the same equation: the
vorticity Navier-Stokes equation. These new governing equations, Eqs. 26 and
27, ensure local conservation of both angular momentum and mass, thus
overcoming the pitfalls of the classic VPM.

####  IV.2.3 Other Formulations

Combining the two conservations laws with the two element shapes, one can
devise up to six different formulations: Each permutation between sphere/tube
elements, momentum/mass conservation, and the assumption that mass
conservation has or has not been already ensured by the convection of the
elements. Two of them have already been discussed and turn out to be the same.
For completeness, we now derive and briefly discuss the four remaining
formulations.

Mass-Conserving Sphere

The mass of a spherical fluid element is conserved in incompressible flow as
long as its volume remains constant. This translates into vortex particles of
constant radii,

$\displaystyle\frac{\mathrm{d}}{\mathrm{d}t}\sigma_{p}=0,$ (31)

which reduces the vortex strength equation to

$\displaystyle\frac{\mathrm{d}}{\mathrm{d}t}\boldsymbol{\Gamma}_{p}$
$\displaystyle=\left(\boldsymbol{\Gamma}_{p}\cdot\nabla\right)\mathbf{u}(\mathbf{x}_{p}).$

Notice that these are the equations of the classic VPM. Even though this
formulation preserves mass of spherical chunks of the flow, it is deemed
unphysical as it violates local conservation of mass when the vorticity field
encompasses a tube of vorticity that is being stretched, as discussed in
Section IV.1.

Momentum-Conserving Tube

Consider a cylindrical element that rotates with angular velocity
${\boldsymbol{\Omega}=\boldsymbol{\omega}/2}$. Let $\mathbb{I}$ denote its
inertia tensor and $I_{z}$ the moment of inertia about the centerline axis. If
we assume that $\boldsymbol{\Omega}$ is always aligned with its centerline
axis, the inertia tensor product $\mathbb{I}\boldsymbol{\Omega}$ becomes
$\mathbb{I}\boldsymbol{\Omega}=I_{z}\boldsymbol{\Omega}$ since
$\boldsymbol{\Omega}$ is an eigenvector of $\mathbb{I}$, with $I_{z}$ the
corresponding eigenvalue. In order to conserve angular momentum, neglecting
pressure torque and following the discussion in Section IV.1.1, the moment of
inertia must change according to

$\displaystyle-\frac{\boldsymbol{\omega}}{I_{z}}\frac{\mathrm{d}}{\mathrm{d}t}I_{z}=(\boldsymbol{\omega}\cdot\nabla)\mathbf{u}.$
(32)

Using $I_{z}=\frac{1}{2}mr^{2}$ and without assuming that the mass $m=\rho\pi
r^{2}\ell$ is automatically conserved, the dot product of Eq. 32 with
$\boldsymbol{\omega}$ leads to

$\displaystyle\frac{\mathrm{d}}{\mathrm{d}t}r=-\frac{r}{4\ell}\frac{\mathrm{d}}{\mathrm{d}t}\ell-\frac{r}{4\|\boldsymbol{\omega}\|^{2}}\left[(\boldsymbol{\omega}\cdot\nabla)\mathbf{u}\right]\cdot\boldsymbol{\omega}.$

In the vortex particle scheme this becomes

$\displaystyle\frac{\mathrm{d}}{\mathrm{d}t}\sigma_{p}=-\frac{1}{4}\frac{\sigma_{p}}{\|\boldsymbol{\Gamma}_{p}\|}\frac{\mathrm{d}\boldsymbol{\Gamma}_{p}}{\mathrm{d}t}\cdot\hat{\boldsymbol{\Gamma}}_{p}-\frac{1}{4}\frac{\sigma_{p}}{\|\boldsymbol{\Gamma}_{p}\|}\left[(\boldsymbol{\Gamma}_{p}\cdot\nabla)\mathbf{u}\right]\cdot\hat{\boldsymbol{\Gamma}}_{p},$
(33)

which in conjunction to the vortex strength equation, Eq. 25, we obtain the
following governing equations:

$\displaystyle\frac{\mathrm{d}}{\mathrm{d}t}\sigma_{p}$
$\displaystyle=-\frac{2}{7}\frac{\sigma_{p}}{\|\boldsymbol{\Gamma}_{p}\|}\left[(\boldsymbol{\Gamma}_{p}\cdot\nabla)\mathbf{u}\right]\cdot\hat{\boldsymbol{\Gamma}}_{p}$
$\displaystyle\frac{\mathrm{d}}{\mathrm{d}t}\boldsymbol{\Gamma}_{p}$
$\displaystyle=\left(\boldsymbol{\Gamma}_{p}\cdot\nabla\right)\mathbf{u}(\mathbf{x}_{p})-\frac{6}{7}\left\\{\left[(\boldsymbol{\Gamma}_{p}\cdot\nabla)\mathbf{u}\right]\cdot\hat{\boldsymbol{\Gamma}}_{p}\right\\}\hat{\boldsymbol{\Gamma}}_{p}.$

Assumed Mass Conservation

Vortex methods typically assume that the continuity equation is implicitly
solved by convecting the elements with the velocity field, automatically
satisfying mass conservation. Repeating the momentum-conserving tube
derivation, but this time assuming that the convection of the particles
automatically ensures $\frac{\mathrm{d}}{\mathrm{d}t}m=0$, we get the
following governing equations:

$\displaystyle\frac{\mathrm{d}}{\mathrm{d}t}\sigma_{p}$
$\displaystyle=-\frac{1}{2}\frac{\sigma_{p}}{\|\boldsymbol{\Gamma}_{p}\|}\left[(\boldsymbol{\Gamma}_{p}\cdot\nabla)\mathbf{u}\right]\cdot\hat{\boldsymbol{\Gamma}}_{p}$
(34) $\displaystyle\frac{\mathrm{d}}{\mathrm{d}t}\boldsymbol{\Gamma}_{p}$
$\displaystyle=\left(\boldsymbol{\Gamma}_{p}\cdot\nabla\right)\mathbf{u}(\mathbf{x}_{p})-\frac{3}{2}\left\\{\left[(\boldsymbol{\Gamma}_{p}\cdot\nabla)\mathbf{u}\right]\cdot\hat{\boldsymbol{\Gamma}}_{p}\right\\}\hat{\boldsymbol{\Gamma}}_{p}.$

Taking the dot product of $\hat{\boldsymbol{\Gamma}}_{p}$ over the vortex
strength equation, we get

$\displaystyle\frac{\mathrm{d}\boldsymbol{\Gamma}_{p}}{\mathrm{d}t}\cdot\hat{\boldsymbol{\Gamma}}_{p}$
$\displaystyle=-\frac{1}{2}\left[(\boldsymbol{\Gamma}_{p}\cdot\nabla)\mathbf{u}\right]\cdot\hat{\boldsymbol{\Gamma}}_{p}.$

In this last equation, notice that the derivative of the strength projected on
its own direction is negative when vortex stretching is positive. This means
that the magnitude of vortex strength (and associated vorticity) decreases
when the particle is stretched, which is an unphysical result.

Along the same lines, we can repeat the momentum-conserving sphere derivation
from Section IV.2.1, but this time assuming automatic mass conservation,
$\frac{\mathrm{d}}{\mathrm{d}t}m=0$. This turns out to lead to the same
governing equations in Eq. 34. Hence, we conclude that it is not possible to
impose momentum conservation with either tube or sphere elements while
assuming automatic mass conservation.

###  IV.3 Generalized Formulation

In the preceding section, the conservation laws were imposed in six different
ways to derive various expressions of
$\frac{\mathrm{d}}{\mathrm{d}t}\sigma_{p}$, namely Eqs. 26, 29, 31, 33 and 34.
In turn, each version of $\frac{\mathrm{d}}{\mathrm{d}t}\sigma_{p}$ was used
with the vortex strength equation, Eq. 25, to obtain a distinct VPM
formulation. We can generalize this formulation procedure noticing that all
versions of $\frac{\mathrm{d}}{\mathrm{d}t}\sigma_{p}$ take the functional
form

$\displaystyle\frac{\mathrm{d}}{\mathrm{d}t}\sigma_{p}$
$\displaystyle=-{\color[rgb]{0.30859375,0.578125,0.80859375}\definecolor[named]{pgfstrokecolor}{rgb}{0.30859375,0.578125,0.80859375}{f}}\sigma_{p}\frac{1}{\|\boldsymbol{\Gamma}_{p}\|}\frac{\mathrm{d}\boldsymbol{\Gamma}_{p}}{\mathrm{d}t}\cdot\hat{\boldsymbol{\Gamma}}_{p}$
(35)
$\displaystyle\qquad-{\color[rgb]{0.62890625,0.55078125,0.4296875}\definecolor[named]{pgfstrokecolor}{rgb}{0.62890625,0.55078125,0.4296875}{g}}\sigma_{p}\frac{1}{\|\boldsymbol{\Gamma}_{p}\|}\left[\left(\boldsymbol{\Gamma}_{p}\cdot\nabla\right)\mathbf{u}{\small(\mathbf{x}_{p})}\right]\cdot\hat{\boldsymbol{\Gamma}}_{p},$

where the parameters
${\color[rgb]{0.30859375,0.578125,0.80859375}\definecolor[named]{pgfstrokecolor}{rgb}{0.30859375,0.578125,0.80859375}{f}},{\color[rgb]{0.62890625,0.55078125,0.4296875}\definecolor[named]{pgfstrokecolor}{rgb}{0.62890625,0.55078125,0.4296875}{g}}\in\mathbb{R}$
are derived applying the conservation laws.

We will now use the
${\color[rgb]{0.30859375,0.578125,0.80859375}\definecolor[named]{pgfstrokecolor}{rgb}{0.30859375,0.578125,0.80859375}{f}}$–${\color[rgb]{0.62890625,0.55078125,0.4296875}\definecolor[named]{pgfstrokecolor}{rgb}{0.62890625,0.55078125,0.4296875}{g}}$
form of $\frac{\mathrm{d}}{\mathrm{d}t}\sigma_{p}$, Eq. 35, to derive a new
set of governing equation of the VPM, also in the
${\color[rgb]{0.30859375,0.578125,0.80859375}\definecolor[named]{pgfstrokecolor}{rgb}{0.30859375,0.578125,0.80859375}{f}}$–${\color[rgb]{0.62890625,0.55078125,0.4296875}\definecolor[named]{pgfstrokecolor}{rgb}{0.62890625,0.55078125,0.4296875}{g}}$
form. We start from the general equation governing vortex strength, Eq. 14,
here repeated

$\displaystyle\frac{\mathrm{d}\boldsymbol{\Gamma}_{p}}{\mathrm{d}t}=\left(\boldsymbol{\Gamma}_{p}\cdot\nabla\right)\mathbf{u}(\mathbf{x}_{p})+3\boldsymbol{\Gamma}_{p}\frac{1}{\sigma_{p}}\frac{\partial\sigma_{p}}{\partial
t}+\mathbf{M}_{p},$ (14)

where

$\displaystyle\mathbf{M}_{p}$
$\displaystyle\equiv\frac{1}{\zeta_{\sigma_{p}}(\mathbf{0})}(-\mathbf{M}_{p}^{0}+\mathbf{M}_{p}^{1}+\mathbf{M}_{p}^{2})$
$\displaystyle\qquad-\frac{1}{\zeta_{\sigma_{p}}(\mathbf{0})}\left(\mathbf{E}_{\mathrm{adv}}(\mathbf{x}_{p})+\mathbf{E}_{\mathrm{str}}(\mathbf{x}_{p})\right).$

Recalling that
$\frac{\mathrm{d}\sigma_{p}}{\mathrm{d}t}=\frac{\partial\sigma_{p}}{\partial
t}$ since $\sigma_{p}$ is a Lagrangian quantity, we substitute Eq. 35 into Eq.
14 and take its dot product with $\hat{\boldsymbol{\Gamma}}_{p}$ to get

$\displaystyle\frac{\mathrm{d}\boldsymbol{\Gamma}_{p}}{\mathrm{d}t}\cdot\hat{\boldsymbol{\Gamma}}_{p}=\frac{1-3{\color[rgb]{0.62890625,0.55078125,0.4296875}\definecolor[named]{pgfstrokecolor}{rgb}{0.62890625,0.55078125,0.4296875}{g}}}{1+3{\color[rgb]{0.30859375,0.578125,0.80859375}\definecolor[named]{pgfstrokecolor}{rgb}{0.30859375,0.578125,0.80859375}{f}}}\left[\left(\boldsymbol{\Gamma}_{p}\cdot\nabla\right)\mathbf{u}{\small(\mathbf{x}_{p})}\right]\cdot\hat{\boldsymbol{\Gamma}}_{p}+\frac{1}{1+3{\color[rgb]{0.30859375,0.578125,0.80859375}\definecolor[named]{pgfstrokecolor}{rgb}{0.30859375,0.578125,0.80859375}{f}}}\mathbf{M}_{p}\cdot\hat{\boldsymbol{\Gamma}}_{p}.$
(36)

Substituting Eq. 36 into Eq. 35, we arrive to the equation governing the
evolution of particle size,

$\displaystyle\frac{\mathrm{d}}{\mathrm{d}t}\sigma_{p}$
$\displaystyle=-\left(\frac{{\color[rgb]{0.62890625,0.55078125,0.4296875}\definecolor[named]{pgfstrokecolor}{rgb}{0.62890625,0.55078125,0.4296875}{g}}+{\color[rgb]{0.30859375,0.578125,0.80859375}\definecolor[named]{pgfstrokecolor}{rgb}{0.30859375,0.578125,0.80859375}{f}}}{1+3{\color[rgb]{0.30859375,0.578125,0.80859375}\definecolor[named]{pgfstrokecolor}{rgb}{0.30859375,0.578125,0.80859375}{f}}}\right)\frac{\sigma_{p}}{\|\boldsymbol{\Gamma}_{p}\|}\left[\left(\boldsymbol{\Gamma}_{p}\cdot\nabla\right)\mathbf{u}{\small(\mathbf{x}_{p})}\right]\cdot\hat{\boldsymbol{\Gamma}}_{p}$
(37)
$\displaystyle\qquad-\left(\frac{{\color[rgb]{0.30859375,0.578125,0.80859375}\definecolor[named]{pgfstrokecolor}{rgb}{0.30859375,0.578125,0.80859375}{f}}}{1+3{\color[rgb]{0.30859375,0.578125,0.80859375}\definecolor[named]{pgfstrokecolor}{rgb}{0.30859375,0.578125,0.80859375}{f}}}\right)\frac{\sigma_{p}}{\|\boldsymbol{\Gamma}_{p}\|}\mathbf{M}_{p}\cdot\hat{\boldsymbol{\Gamma}}_{p},$

and we substitute this into Eq. 14 to get the equation governing the evolution
of vortex strength,

$\displaystyle\begin{split}\frac{\mathrm{d}}{\mathrm{d}t}\boldsymbol{\Gamma}_{p}=\left(\boldsymbol{\Gamma}_{p}\cdot\nabla\right)\mathbf{u}(\mathbf{x}_{p})&-\frac{{\color[rgb]{0.62890625,0.55078125,0.4296875}\definecolor[named]{pgfstrokecolor}{rgb}{0.62890625,0.55078125,0.4296875}{g}}+{\color[rgb]{0.30859375,0.578125,0.80859375}\definecolor[named]{pgfstrokecolor}{rgb}{0.30859375,0.578125,0.80859375}{f}}}{\nicefrac{{1}}{{3}}+{\color[rgb]{0.30859375,0.578125,0.80859375}\definecolor[named]{pgfstrokecolor}{rgb}{0.30859375,0.578125,0.80859375}{f}}}\left\\{\left[\left(\boldsymbol{\Gamma}_{p}\cdot\nabla\right)\mathbf{u}{\small(\mathbf{x}_{p})}\right]\cdot\hat{\boldsymbol{\Gamma}}_{p}\right\\}\hat{\boldsymbol{\Gamma}}_{p}\\\
+\mathbf{M}_{p}&-\frac{{\color[rgb]{0.30859375,0.578125,0.80859375}\definecolor[named]{pgfstrokecolor}{rgb}{0.30859375,0.578125,0.80859375}{f}}}{\nicefrac{{1}}{{3}}+{\color[rgb]{0.30859375,0.578125,0.80859375}\definecolor[named]{pgfstrokecolor}{rgb}{0.30859375,0.578125,0.80859375}{f}}}\left(\mathbf{M}_{p}\cdot\hat{\boldsymbol{\Gamma}}_{p}\right)\hat{\boldsymbol{\Gamma}}_{p}.\end{split}$
(38)

Eqs. 37 and 38 constitute the governing equations of the reformulated VPM,
where the parameters
${\color[rgb]{0.30859375,0.578125,0.80859375}\definecolor[named]{pgfstrokecolor}{rgb}{0.30859375,0.578125,0.80859375}{f}},\,{\color[rgb]{0.62890625,0.55078125,0.4296875}\definecolor[named]{pgfstrokecolor}{rgb}{0.62890625,0.55078125,0.4296875}{g}}$
are derived from a specific implementation of the conservations laws. In
particular when
${{\color[rgb]{0.30859375,0.578125,0.80859375}\definecolor[named]{pgfstrokecolor}{rgb}{0.30859375,0.578125,0.80859375}{f}}={\color[rgb]{0.62890625,0.55078125,0.4296875}\definecolor[named]{pgfstrokecolor}{rgb}{0.62890625,0.55078125,0.4296875}{g}}=0}$,
Eqs. 37 and 38 collapse back to the classic VPM equations, making the
reformulated VPM a generalization of the classic method. Furthermore, notice
that the reformulated equations do not require more computation than the
classic method: When SFS and non-localized vorticity effects are neglected
($\mathbf{M}_{p}=0$), both $\frac{\mathrm{d}\sigma_{p}}{\mathrm{d}t}$ and
$\frac{\mathrm{d}\boldsymbol{\Gamma}_{p}}{\mathrm{d}t}$ are calculated
directly from vortex stretching,
$\left(\boldsymbol{\Gamma}_{p}\cdot\nabla\right)\mathbf{u}(\mathbf{x}_{p})$.

###  IV.4 Formulation Analysis

In Section IV.2 we have shown several plausible formulations that can be
derived from the conservation laws. We now cast each formulation into its
${\color[rgb]{0.30859375,0.578125,0.80859375}\definecolor[named]{pgfstrokecolor}{rgb}{0.30859375,0.578125,0.80859375}{f}}$–${\color[rgb]{0.62890625,0.55078125,0.4296875}\definecolor[named]{pgfstrokecolor}{rgb}{0.62890625,0.55078125,0.4296875}{g}}$
form, as summarized in Table 1, and contrast the physical implications of each
formulation to help us choose one formulation to further explore.

We wish to determine how each formulation may hinder or augment the effects of
vortex stretching on both vortex strength and vortex re-orientation. First, we
focus on vortex re-orientation. The rate of vortex re-orientation, denoted
$\frac{\mathrm{d}\hat{\boldsymbol{\Gamma}}_{p}}{\mathrm{d}t}$, is constructed
as

$\displaystyle\frac{\mathrm{d}\hat{\boldsymbol{\Gamma}}_{p}}{\mathrm{d}t}=\frac{1}{\|\boldsymbol{\Gamma}_{p}\|}\frac{\mathrm{d}\boldsymbol{\Gamma}_{p}}{\mathrm{d}t}-\frac{1}{\|\boldsymbol{\Gamma}_{p}\|}\left(\frac{\mathrm{d}\boldsymbol{\Gamma}_{p}}{\mathrm{d}t}\cdot\hat{\boldsymbol{\Gamma}}_{p}\right)\hat{\boldsymbol{\Gamma}}_{p},$

which is the change in direction of vortex strength over time, as shown in
Appendix B. Using Eq. 36 and Eq. 38 and neglecting $\mathbf{M}$, this becomes

$\displaystyle\frac{\mathrm{d}\hat{\boldsymbol{\Gamma}}_{p}}{\mathrm{d}t}=\left(\hat{\boldsymbol{\Gamma}}_{p}\cdot\nabla\right)\mathbf{u}(\mathbf{x}_{p})-\left\\{\left[\left(\hat{\boldsymbol{\Gamma}}_{p}\cdot\nabla\right)\mathbf{u}{\small(\mathbf{x}_{p})}\right]\cdot\hat{\boldsymbol{\Gamma}}_{p}\right\\}\hat{\boldsymbol{\Gamma}}_{p},$

which is independent of
${\color[rgb]{0.30859375,0.578125,0.80859375}\definecolor[named]{pgfstrokecolor}{rgb}{0.30859375,0.578125,0.80859375}{f}}$
and
${\color[rgb]{0.62890625,0.55078125,0.4296875}\definecolor[named]{pgfstrokecolor}{rgb}{0.62890625,0.55078125,0.4296875}{g}}$.
Hence, vortex re-orientation is not affected by the formulation, staying the
same between the classic VPM and any
${\color[rgb]{0.30859375,0.578125,0.80859375}\definecolor[named]{pgfstrokecolor}{rgb}{0.30859375,0.578125,0.80859375}{f}}$–${\color[rgb]{0.62890625,0.55078125,0.4296875}\definecolor[named]{pgfstrokecolor}{rgb}{0.62890625,0.55078125,0.4296875}{g}}$
formulation.

Table 1: ${\color[rgb]{0.30859375,0.578125,0.80859375}\definecolor[named]{pgfstrokecolor}{rgb}{0.30859375,0.578125,0.80859375}{f}}$–${\color[rgb]{0.62890625,0.55078125,0.4296875}\definecolor[named]{pgfstrokecolor}{rgb}{0.62890625,0.55078125,0.4296875}{g}}$ form of formulations derived in Section IV.2 Formulation | ${\color[rgb]{0.30859375,0.578125,0.80859375}\definecolor[named]{pgfstrokecolor}{rgb}{0.30859375,0.578125,0.80859375}{f}}$ | ${\color[rgb]{0.62890625,0.55078125,0.4296875}\definecolor[named]{pgfstrokecolor}{rgb}{0.62890625,0.55078125,0.4296875}{g}}$ | $h_{\Gamma}=\frac{1-3{\color[rgb]{0.62890625,0.55078125,0.4296875}\definecolor[named]{pgfstrokecolor}{rgb}{0.62890625,0.55078125,0.4296875}{g}}}{1+3{\color[rgb]{0.30859375,0.578125,0.80859375}\definecolor[named]{pgfstrokecolor}{rgb}{0.30859375,0.578125,0.80859375}{f}}}$ | $h_{\sigma}=\frac{{\color[rgb]{0.62890625,0.55078125,0.4296875}\definecolor[named]{pgfstrokecolor}{rgb}{0.62890625,0.55078125,0.4296875}{g}}+{\color[rgb]{0.30859375,0.578125,0.80859375}\definecolor[named]{pgfstrokecolor}{rgb}{0.30859375,0.578125,0.80859375}{f}}}{1+3{\color[rgb]{0.30859375,0.578125,0.80859375}\definecolor[named]{pgfstrokecolor}{rgb}{0.30859375,0.578125,0.80859375}{f}}}$
---|---|---|---|---
| Mass-conserving sphere
---
(Classic VPM)
0 | 0 | 1 | 0
Mass-conserving tube | ${\color[rgb]{0.30859375,0.578125,0.80859375}\definecolor[named]{pgfstrokecolor}{rgb}{0.30859375,0.578125,0.80859375}{\nicefrac{{1}}{{2}}}}$ | 0 | 0.4 | 0.2
| Momentum-conserving sphere
---
(Reformulated VPM)
0 | ${\color[rgb]{0.62890625,0.55078125,0.4296875}\definecolor[named]{pgfstrokecolor}{rgb}{0.62890625,0.55078125,0.4296875}{\nicefrac{{1}}{{5}}}}$ | 0.4 | 0.2
Momentum-conserving tube | ${\color[rgb]{0.30859375,0.578125,0.80859375}\definecolor[named]{pgfstrokecolor}{rgb}{0.30859375,0.578125,0.80859375}{\nicefrac{{1}}{{4}}}}$ | ${\color[rgb]{0.62890625,0.55078125,0.4296875}\definecolor[named]{pgfstrokecolor}{rgb}{0.62890625,0.55078125,0.4296875}{\nicefrac{{1}}{{4}}}}$ | 0.143 | 0.286
| Momentum-conservation
---
w/ assumed mass conservation
0 | ${\color[rgb]{0.62890625,0.55078125,0.4296875}\definecolor[named]{pgfstrokecolor}{rgb}{0.62890625,0.55078125,0.4296875}{\nicefrac{{1}}{{2}}}}$ | -0.5 | 0.5

Now, we shift our attention to the effects of vortex stretching on vortex
strength. Neglecting SFS and non-localized vorticity effects
($\mathbf{M}_{p}=0$), from Eq. 36 we have

$\displaystyle\frac{\mathrm{d}\boldsymbol{\Gamma}_{p}}{\mathrm{d}t}\cdot\hat{\boldsymbol{\Gamma}}_{p}=\frac{1-3{\color[rgb]{0.62890625,0.55078125,0.4296875}\definecolor[named]{pgfstrokecolor}{rgb}{0.62890625,0.55078125,0.4296875}{g}}}{1+3{\color[rgb]{0.30859375,0.578125,0.80859375}\definecolor[named]{pgfstrokecolor}{rgb}{0.30859375,0.578125,0.80859375}{f}}}\left[\left(\boldsymbol{\Gamma}_{p}\cdot\nabla\right)\mathbf{u}{\small(\mathbf{x}_{p})}\right]\cdot\hat{\boldsymbol{\Gamma}}_{p},$
(39)

where
$\frac{\mathrm{d}\boldsymbol{\Gamma}_{p}}{\mathrm{d}t}\cdot\hat{\boldsymbol{\Gamma}}_{p}$
is the change in the magnitude of vortex strength. We introduce the parameters

$\displaystyle
h_{\Gamma}\equiv\frac{1-3{\color[rgb]{0.62890625,0.55078125,0.4296875}\definecolor[named]{pgfstrokecolor}{rgb}{0.62890625,0.55078125,0.4296875}{g}}}{1+3{\color[rgb]{0.30859375,0.578125,0.80859375}\definecolor[named]{pgfstrokecolor}{rgb}{0.30859375,0.578125,0.80859375}{f}}},\quad
h_{\sigma}\equiv\frac{{\color[rgb]{0.62890625,0.55078125,0.4296875}\definecolor[named]{pgfstrokecolor}{rgb}{0.62890625,0.55078125,0.4296875}{g}}+{\color[rgb]{0.30859375,0.578125,0.80859375}\definecolor[named]{pgfstrokecolor}{rgb}{0.30859375,0.578125,0.80859375}{f}}}{1+3{\color[rgb]{0.30859375,0.578125,0.80859375}\definecolor[named]{pgfstrokecolor}{rgb}{0.30859375,0.578125,0.80859375}{f}}}$

and rewrite this as

$\displaystyle\frac{\mathrm{d}\boldsymbol{\Gamma}_{p}}{\mathrm{d}t}\cdot\hat{\boldsymbol{\Gamma}}_{p}=h_{\Gamma}\left[\left(\boldsymbol{\Gamma}_{p}\cdot\nabla\right)\mathbf{u}{\small(\mathbf{x}_{p})}\right]\cdot\hat{\boldsymbol{\Gamma}}_{p},$

while Eq. 37 becomes

$\displaystyle\frac{\|\boldsymbol{\Gamma}_{p}\|}{\sigma_{p}}\frac{\mathrm{d}}{\mathrm{d}t}\sigma_{p}=-h_{\sigma}\left[\left(\boldsymbol{\Gamma}_{p}\cdot\nabla\right)\mathbf{u}{\small(\mathbf{x}_{p})}\right]\cdot\hat{\boldsymbol{\Gamma}}_{p}.$

It can be shown that $h_{\Gamma}$ and $h_{\sigma}$ must satisfy
${h_{\Gamma}+3h_{\sigma}=1}$ in order for a formulation to be a solution to
the vorticity Navier-Stokes equation. Hence, any
${\color[rgb]{0.30859375,0.578125,0.80859375}\definecolor[named]{pgfstrokecolor}{rgb}{0.30859375,0.578125,0.80859375}{f}}$–${\color[rgb]{0.62890625,0.55078125,0.4296875}\definecolor[named]{pgfstrokecolor}{rgb}{0.62890625,0.55078125,0.4296875}{g}}$
formulation simply differs in how vortex stretching is distributed between the
lengthening of vortex strength $\boldsymbol{\Gamma}_{p}$, and shrinking of
particle size $\sigma_{p}$. The classic VPM results in $h_{\Gamma}=1$ and
$h_{\sigma}=0$, meaning that vortex strength takes the full amount of vortex
stretching.

The space of all possible
${\color[rgb]{0.30859375,0.578125,0.80859375}\definecolor[named]{pgfstrokecolor}{rgb}{0.30859375,0.578125,0.80859375}{f}}$–${\color[rgb]{0.62890625,0.55078125,0.4296875}\definecolor[named]{pgfstrokecolor}{rgb}{0.62890625,0.55078125,0.4296875}{g}}$
formulations is shown in Fig. 5, along with the contours of $h_{\Gamma}$ that
satisfy ${h_{\Gamma}+3h_{\sigma}=1}$. Any formulation with $h_{\Gamma}>1$ will
result in a negative $h_{\sigma}$, meaning that the lengthening of vortex
strength under stretching is so aggressive that $\sigma_{p}$ has to expand
instead of shrink, which is unphysical. This is referred to as overstretching.
On the other extreme, $h_{\sigma}>\nicefrac{{1}}{{3}}$ results in a negative
$h_{\Gamma}$, meaning that the shrinking of the particle is so aggressive that
$\boldsymbol{\Gamma}$ has to compensate by shrinking when it ought to have
lengthened with vortex stretching, which is also unphysical. This is referred
to as anti-stretching. The region between overstretching and anti-stretching
is the space of VPM formulations that are physically plausible.

Figure 5: Space of possible
${\color[rgb]{0.30859375,0.578125,0.80859375}\definecolor[named]{pgfstrokecolor}{rgb}{0.30859375,0.578125,0.80859375}{f}}$-${\color[rgb]{0.62890625,0.55078125,0.4296875}\definecolor[named]{pgfstrokecolor}{rgb}{0.62890625,0.55078125,0.4296875}{g}}$
formulations, with contour levels of
$h_{\Gamma}=\frac{1-3{\color[rgb]{0.62890625,0.55078125,0.4296875}\definecolor[named]{pgfstrokecolor}{rgb}{0.62890625,0.55078125,0.4296875}{g}}}{1+3{\color[rgb]{0.30859375,0.578125,0.80859375}\definecolor[named]{pgfstrokecolor}{rgb}{0.30859375,0.578125,0.80859375}{f}}}$
that are solutions to the vorticity Navier-Stokes equation. Markers show each
formulation in Table 1.

There are four formulations that lay in the space of physically plausible
formulations, shown in Fig. 5. First, notice that the classic VPM (mass-
conserving sphere) lies at the threshold of overstretching, which explains its
tendency to be numerically unstable. We are then left with three possible
formulations, two of which lay on the same $h_{\Gamma}$ contour line: the
momentum-conserving sphere and the mass-conserving tube. Even though these two
formulations stem from different physical principles resulting in different
${\color[rgb]{0.30859375,0.578125,0.80859375}\definecolor[named]{pgfstrokecolor}{rgb}{0.30859375,0.578125,0.80859375}{f}}$–${\color[rgb]{0.62890625,0.55078125,0.4296875}\definecolor[named]{pgfstrokecolor}{rgb}{0.62890625,0.55078125,0.4296875}{g}}$
values, they share the same $h_{\Gamma}$ and $h_{\sigma}$ values, leading to
the same governing equations******This is only true after neglecting
$\mathbf{M}$. When SFS and non-localized vorticity effects are re-
incorporated, the $\mathbf{M}$-term ends up being multiplied by a factor
${h_{M}\equiv\frac{1}{1+3{\color[rgb]{0.30859375,0.578125,0.80859375}\definecolor[named]{pgfstrokecolor}{rgb}{0.30859375,0.578125,0.80859375}{f}}}}$
as shown in Eq. 36. Hence, SFS and non-localized vorticity effects are
stronger in the momentum-conserving sphere formulation (${h_{M}=1}$) than the
mass-conserving tube formulation (${h_{M}=0.4}$)., as mentioned in Section
IV.2.2. Thus, by choosing either set of
${\color[rgb]{0.30859375,0.578125,0.80859375}\definecolor[named]{pgfstrokecolor}{rgb}{0.30859375,0.578125,0.80859375}{f}}$–${\color[rgb]{0.62890625,0.55078125,0.4296875}\definecolor[named]{pgfstrokecolor}{rgb}{0.62890625,0.55078125,0.4296875}{g}}$
values, we end up in a formulation that preserves angular momentum in
spherical chunks of the flow, while at the same time conserving mass in tubes
of vorticity. Here on, the momentum-conserving sphere formulation is referred
to as the reformulated vortex particle method, or rVPM.

###  IV.5 Reformulated Governing Equations

Summarizing the preceding sections, the governing equations of the
reformulated VPM are

$\displaystyle\frac{\text{d}}{\text{d}t}\mathbf{x}_{p}$
$\displaystyle=\mathbf{u}(\mathbf{x}_{p})$ (40)
$\displaystyle\frac{\mathrm{d}}{\mathrm{d}t}\sigma_{p}$
$\displaystyle=-\left(\frac{{\color[rgb]{0.62890625,0.55078125,0.4296875}\definecolor[named]{pgfstrokecolor}{rgb}{0.62890625,0.55078125,0.4296875}{g}}+{\color[rgb]{0.30859375,0.578125,0.80859375}\definecolor[named]{pgfstrokecolor}{rgb}{0.30859375,0.578125,0.80859375}{f}}}{1+3{\color[rgb]{0.30859375,0.578125,0.80859375}\definecolor[named]{pgfstrokecolor}{rgb}{0.30859375,0.578125,0.80859375}{f}}}\right)\frac{\sigma_{p}}{\|\boldsymbol{\Gamma}_{p}\|}\left[\left(\boldsymbol{\Gamma}_{p}\cdot\nabla\right)\mathbf{u}{\small(\mathbf{x}_{p})}\right]\cdot\hat{\boldsymbol{\Gamma}}_{p}$
(41) $\displaystyle\frac{\mathrm{d}}{\mathrm{d}t}\boldsymbol{\Gamma}_{p}$
$\displaystyle=\left(\boldsymbol{\Gamma}_{p}\cdot\nabla\right)\mathbf{u}(\mathbf{x}_{p})-\frac{{\color[rgb]{0.62890625,0.55078125,0.4296875}\definecolor[named]{pgfstrokecolor}{rgb}{0.62890625,0.55078125,0.4296875}{g}}+{\color[rgb]{0.30859375,0.578125,0.80859375}\definecolor[named]{pgfstrokecolor}{rgb}{0.30859375,0.578125,0.80859375}{f}}}{\nicefrac{{1}}{{3}}+{\color[rgb]{0.30859375,0.578125,0.80859375}\definecolor[named]{pgfstrokecolor}{rgb}{0.30859375,0.578125,0.80859375}{f}}}\left\\{\left[\left(\boldsymbol{\Gamma}_{p}\cdot\nabla\right)\mathbf{u}{\small(\mathbf{x}_{p})}\right]\cdot\hat{\boldsymbol{\Gamma}}_{p}\right\\}\hat{\boldsymbol{\Gamma}}_{p}$
(42)
$\displaystyle\left(\frac{\text{d}}{\text{d}t}\overline{\boldsymbol{\omega}}\right)_{\mathrm{viscous}}=\nu\nabla^{2}\overline{\boldsymbol{\omega}},$
(43)

where Eq. 40 resolves vorticity advection by convecting the particles, Eq. 41
governs the evolution of particle size, and Eq. 42 governs the evolution of
vortex strength. Eq. 42 in conjunction with Eqs. 40 and 41 resolve the
inviscid part of the LES-filtered vorticity Navier-Stokes equation, while the
viscous part in Eq. 43 is resolved by one of the schemes mentioned in Section
II.3. For simplicity, Eqs. 41 and 42 omit the SFS contributions, but they are
readily re-incorporated as shown in Eqs. 37 and 38.

The formulation used for the rest of this study implements
${\color[rgb]{0.30859375,0.578125,0.80859375}\definecolor[named]{pgfstrokecolor}{rgb}{0.30859375,0.578125,0.80859375}{f}}=0$
and
${\color[rgb]{0.62890625,0.55078125,0.4296875}\definecolor[named]{pgfstrokecolor}{rgb}{0.62890625,0.55078125,0.4296875}{g}}=\nicefrac{{1}}{{5}}$,
which ensures conservation of angular momentum in spherical chunks of the
flow, while at the same time conserving mass in tubes of vorticity. As
previously mentioned, the rVPM equations do not require more computation than
the classic VPM since $\frac{\mathrm{d}\sigma_{p}}{\mathrm{d}t}$ and
$\frac{\mathrm{d}\boldsymbol{\Gamma}_{p}}{\mathrm{d}t}$ are calculated
directly and solely from vortex stretching,
$\left(\boldsymbol{\Gamma}_{p}\cdot\nabla\right)\mathbf{u}(\mathbf{x}_{p})$.

##  V Anisotropic Dynamic SFS Model

In the previous section we have formulated a numerical scheme for solving the
large scales of the LES-filtered vorticity equation,

$\displaystyle\frac{\text{d}}{\text{d}t}\overline{\boldsymbol{\omega}}=\left(\overline{\boldsymbol{\omega}}\cdot\nabla\right)\overline{\mathbf{u}}+\nu\nabla^{2}\overline{\boldsymbol{\omega}}-\mathbf{E}_{\mathrm{adv}}-\mathbf{E}_{\mathrm{str}}.$
(5)

We will now focus on the subfilter-scale (SFS) stresses associated with
advection and vortex stretching, $\mathbf{E}_{\mathrm{adv}}$ and
$\mathbf{E}_{\mathrm{str}}$, respectively. In Section V.1 we briefly discuss
existing SFS models relevant to meshless vortex methods. In Section V.2 we
introduce a novel anisotropic structural model of SFS vortex stretching that
is suitable for turbulent flows where the predominant cascade mechanism is
vortex stretching. In Sections V.3 and V.4 we develop a dynamic procedure for
computing the model coefficient while also providing the means for backscatter
control. Finally, in Section V.5 we show how our model can be implemented in
conventional mesh-based CFD with a pressure-velocity solver.

###  V.1 Existing Models

Over the years, only a few models have been proposed that are suitable for
meshless vortex methods. The most popular one is a variant of the Smagorinsky
eddy-viscosity model formulated for the vorticity stress Winckelmans (1995);
Mansfield, Knio, and Meneveau (1998, 1996, 1999),

$\displaystyle\mathbf{E}=-\nabla\times\left(\nu_{\mathrm{SFS}}\nabla\times\overline{\boldsymbol{\omega}}\right),$

where $\mathbf{E}\equiv\mathbf{E}_{\mathrm{adv}}+\mathbf{E}_{\mathrm{str}}$,
$\nu_{\mathrm{SFS}}=C_{d}^{2}\sigma^{2}\sqrt{2S_{mn}S_{mn}}$, $S_{mn}$ is the
strain-rate tensor, $\sigma$ is the filter width, and $C_{d}$ is a model
coefficient which is either prescribed or computed dynamically. This
functional model and others alike were developed on the basis of homogeneous
isotropic turbulence, which makes them overly diffusive in simulations with
coherent vortical structures. In the latest developments of the vortex
particle-mesh schemeCaprace, Winckelmans, and Chatelain (2020, 2021), this
drawback has been avoided with the variational multiscale methodHughes,
Mazzei, and Jansen (2000); Vreman (2003); Jeanmart and Winckelmans (2007);
Cocle, Bricteux, and Winckelmans (2009), however its applicability to a
meshless scheme is not clear.

In a different approach, CottetCottet (1996, 1998); Cottet, Jiroveanu, and
Michaux (2003) developed an anisotropic structural model of the advective SFS
term $\mathbf{E}_{\mathrm{adv}}$ as

$\displaystyle\mathbf{E}_{\mathrm{adv}}(\mathbf{x}_{p})=$ $\displaystyle\quad
C_{d}\sum\limits_{q}\mathrm{Vol}_{q}\left(\boldsymbol{\omega}_{p}-\boldsymbol{\omega}_{q}\right)\left[\left(\mathbf{u}(\mathbf{x}_{p})-\mathbf{u}(\mathbf{x}_{q})\right)\cdot\nabla\zeta_{\sigma_{q}}(\mathbf{x}_{p}-\mathbf{x}_{q})\right],$

where the model coefficient $C_{d}$ is usually prescribed with ad hoc
calibration Mimeau, Cottet, and Mortazavi (2019). This model is
reportedlyStock, Gharakhani, and Stone (2010); Mimeau, Cottet, and Mortazavi
(2018) less dissipative than Smagorinsky-type models, however, it was
developed on the basis of 2D flow which is absent of vortex stretching. Thus,
this model neglects $\mathbf{E}_{\mathrm{str}}$ even though vortex stretching
is known to be one of the main mechanisms for enstrophy production in the
energy cascade in three dimensions.

To address the need for a low-dissipation SFS model that captures vortex
stretching as the physical mechanism for turbulence, in the following sections
we develop an anisotropic structural model of SFS vortex stretching with a
dynamic model coefficient.

###  V.2 SFS Vortex Stretching Model

Starting from the definition of SFS vortex stretching,

$\displaystyle E_{i}^{\mathrm{str}}=-\frac{\partial T_{ij}}{\partial
x_{j}}\equiv-\left(\overline{\omega_{j}\frac{\partial u_{i}}{\partial
x_{j}}}-\overline{\omega_{j}}\frac{\partial\overline{u_{i}}}{\partial
x_{j}}\right),$ (44)

we group the filter operator as

$\displaystyle E_{i}^{\mathrm{str}}\left(\mathbf{x}\right)=-$
$\displaystyle\int\limits\omega_{j}\left(\mathbf{y}\right)\left(\frac{\partial
u_{i}}{\partial
x_{j}}{\small\left(\mathbf{y}\right)}-\frac{\partial\overline{u_{i}}}{\partial
x_{j}}\left(\mathbf{x}\right)\right)\zeta_{\sigma}(\mathbf{x}-\mathbf{y})\,\mathrm{d}\mathbf{y}.$

Using the singular particle approximation,
$\boldsymbol{\omega}(\mathbf{x})\approx\sum\limits_{q}\boldsymbol{\Gamma}_{q}\delta(\mathbf{x}-\mathbf{x}_{q})$,
the integral collapses to

$\displaystyle
E_{i}^{\mathrm{str}}\left(\mathbf{x}\right)\approx-\sum\limits_{q}\Gamma_{j}^{q}\left(\frac{\partial
u_{i}}{\partial
x_{j}}\left(\mathbf{x}_{q}\right)-\frac{\partial\overline{u_{i}}}{\partial
x_{j}}\left(\mathbf{x}\right)\right)\zeta_{\sigma}(\mathbf{x}-\mathbf{x}_{q}).$

Assuming $\frac{\partial u_{i}}{\partial
x_{j}}\approx\frac{\partial\overline{u_{i}}}{\partial x_{j}}$ and writing in
vector notation, the model then becomes

$\displaystyle\mathbf{E}_{\mathrm{str}}\left(\mathbf{x}\right)\approx\sum\limits_{q}\zeta_{\sigma}(\mathbf{x}-\mathbf{x}_{q})\left(\boldsymbol{\Gamma}_{q}\cdot\nabla\right)\left(\overline{\mathbf{u}}\left(\mathbf{x}\right)-\overline{\mathbf{u}}\left(\mathbf{x}_{q}\right)\right).$

Interestingly, both Cottet’s model and our model are strikingly similar to the
non-localized vorticity terms $\mathbf{M}_{p}^{1}$ and $\mathbf{M}_{p}^{2}$,
respectively (defined in Eqs. 11 and 12), which are neglected (or filtered
out) by the localized-vorticity assumption. Thus, these models can be thought
of as soft deconvolutions of the localized-vorticity assumption. Furthermore,
introducing a model coefficient, $C_{d}$, they can approximate full
deconvolutions encompassing the entire spectrum of subfilter scales. This
model coefficient is computed dynamically as follows.

###  V.3 Dynamic Procedure

####  V.3.1 Derivative Balance

Since SFS vortex stretching needs to be modeled from resolved quantities, the
initial equality given in Eq. 44, written here in vector notation

$\displaystyle\mathbf{E}_{\mathrm{str}}=-\left[\overline{(\boldsymbol{\omega}\cdot\nabla)\mathbf{u}}-(\overline{\boldsymbol{\omega}}\cdot\nabla)\overline{\mathbf{u}}\right],$

becomes only an approximation,

$\displaystyle\mathbf{E}_{\mathrm{str}}\approx-\left[\overline{(\boldsymbol{\omega}\cdot\nabla)\mathbf{u}}-(\overline{\boldsymbol{\omega}}\cdot\nabla)\overline{\mathbf{u}}\right]$

once the $\mathbf{E}_{\mathrm{str}}$ model is introduced. To recover the
equality (or at least improve the approximation), we introduce a dynamic model
coefficient $C_{d}(\mathbf{x},\,t)$ satisfying

$\displaystyle C_{d}\mathbf{E}_{\mathrm{str}}$
$\displaystyle=-\left[\overline{(\boldsymbol{\omega}\cdot\nabla)\mathbf{u}}-(\overline{\boldsymbol{\omega}}\cdot\nabla)\overline{\mathbf{u}}\right].$
(45)

However, this equation is not useful per se since, in order to determine
$C_{d}$, Eq. 45 requires knowing the SFS quantity
${\overline{(\boldsymbol{\omega}\cdot\nabla)\mathbf{u}}-(\overline{\boldsymbol{\omega}}\cdot\nabla)\overline{\mathbf{u}}}$,
which is exactly what we are trying to model with $\mathbf{E}_{\mathrm{str}}$.
However, assuming scale similarity (meaning that $C_{d}$ is independent of the
filter width), we differentiate this equation with respect to the filter width
$\sigma$ to obtain a more useful relation:

$\displaystyle
C_{d}\frac{\partial\mathbf{E}_{\mathrm{str}}}{\partial\sigma}=-\frac{\partial}{\partial\sigma}\left[\overline{(\boldsymbol{\omega}\cdot\nabla)\mathbf{u}}-(\overline{\boldsymbol{\omega}}\cdot\nabla)\overline{\mathbf{u}}\right].$
(46)

Using the singular particle approximation, the right-hand side of Eq. 46
before differentiation becomes

$\displaystyle\left.\left[\overline{(\boldsymbol{\omega}\cdot\nabla)\mathbf{u}}-(\overline{\boldsymbol{\omega}}\cdot\nabla)\overline{\mathbf{u}}\right]\right|_{\mathbf{x}}\approx$
$\displaystyle\qquad\qquad\sum\limits_{q}\zeta_{\sigma}(\mathbf{x}-\mathbf{x}_{q})\left(\boldsymbol{\Gamma}_{q}\cdot\nabla\right)\left(\mathbf{u}\left(\mathbf{x}_{q}\right)-\overline{\mathbf{u}}\left(\mathbf{x}\right)\right),$

and after differentiation,

$\displaystyle\frac{\partial}{\partial\sigma}\left.\left[\overline{(\boldsymbol{\omega}\cdot\nabla)\mathbf{u}}-(\overline{\boldsymbol{\omega}}\cdot\nabla)\overline{\mathbf{u}}\right]\right|_{\mathbf{x}}\approx$
$\displaystyle\qquad\qquad\sum\limits_{q}\frac{\partial\zeta_{\sigma}}{\partial\sigma}(\mathbf{x}-\mathbf{x}_{q})\left(\boldsymbol{\Gamma}_{q}\cdot\nabla\right)\left(\mathbf{u}\left(\mathbf{x}_{q}\right)-\overline{\mathbf{u}}\left(\mathbf{x}\right)\right)$
$\displaystyle\qquad\qquad-\sum\limits_{q}\zeta_{\sigma}(\mathbf{x}-\mathbf{x}_{q})\left(\boldsymbol{\Gamma}_{q}\cdot\nabla\right)\frac{\partial\overline{\mathbf{u}}}{\partial\sigma}\left(\mathbf{x}\right).$

Recalling that $\zeta_{\sigma}$ is defined as
$\zeta_{\sigma}(\mathbf{x})=\frac{1}{\sigma^{3}}\zeta\left(\frac{\|\mathbf{x}\|}{\sigma}\right)$,
its width-derivative is

$\displaystyle\frac{\partial\zeta_{\sigma}}{\partial\sigma}(\mathbf{x})$
$\displaystyle=-\frac{3}{\sigma}\zeta_{\sigma}\left(\mathbf{x}\right)-\frac{\|\mathbf{x}\|}{\sigma^{5}}\frac{\partial\zeta}{\partial
r}\left(\frac{\|\mathbf{x}\|}{\sigma}\right),$

and assuming that $\zeta_{\sigma}(\mathbf{x})$ reaches a maximum at
$\mathbf{x}=\mathbf{0}$ (meaning, ${\frac{\partial\zeta}{\partial r}(0)=0}$),

$\displaystyle\frac{\partial\zeta_{\sigma}}{\partial\sigma}(\mathbf{0})=-\frac{3}{\sigma}\zeta_{\sigma}\left(\mathbf{0}\right).$

Then, evaluating at $\mathbf{x}=\mathbf{x}_{p}$, assuming
$\zeta_{\sigma}\left(\mathbf{0}\right)\neq 0$, and using the localized-
vorticity assumption to neglect all terms $q\neq p$, Eq. 46 becomes

$\displaystyle
C_{d}(\mathbf{x}_{p})\frac{1}{\zeta_{\sigma}\left(\mathbf{0}\right)}\frac{\partial\mathbf{E}_{\mathrm{str}}}{\partial\sigma}(\mathbf{x}_{p})=$
$\displaystyle\qquad\quad\frac{3}{\sigma}\left(\boldsymbol{\Gamma}_{p}\cdot\nabla\right)\left(\mathbf{u}(\mathbf{x}_{p})-\overline{\mathbf{u}}(\mathbf{x}_{p})\right)+\left(\boldsymbol{\Gamma}_{p}\cdot\nabla\right)\frac{\partial\overline{\mathbf{u}}}{\partial\sigma}(\mathbf{x}_{p}),$

or

$\displaystyle C_{d}\mathbf{m}=\mathbf{L},$ (47)

with

$\displaystyle\mathbf{m}\equiv\frac{\sigma^{3}}{\zeta(0)}\frac{\partial\mathbf{E}_{\mathrm{str}}}{\partial\sigma}(\mathbf{x}_{p})$
$\displaystyle\mathbf{L}\equiv\frac{3}{\sigma}\left(\boldsymbol{\Gamma}_{p}\cdot\nabla\right)\left(\mathbf{u}(\mathbf{x}_{p})-\overline{\mathbf{u}}(\mathbf{x}_{p})\right)+\left(\boldsymbol{\Gamma}_{p}\cdot\nabla\right)\frac{\partial\overline{\mathbf{u}}}{\partial\sigma}(\mathbf{x}_{p}).$

The relation in Eq. 47 will be the basis for our dynamic procedure. This
procedure differs from the classic dynamic procedure developed by Germano et
al.Germano  _et al._ (1991); Germano (1992) in that the former is motivated by
the balance of derivatives between true and modeled SFS contributions given in
Eq. 46, while the latter is based on the Germano identity.

####  V.3.2 Enstrophy-Production Balance

The procedure aims at calculating $C_{d}$ such as to impose the relation given
in Eq. 47, however, this is an overdetermined system as there are three
equations (one for each spatial dimension) and only one unknown, $C_{d}$.
Thus, the relation is now contracted by also imposing a balance of enstrophy
production between true and modeled SFS contributions as follows.

Enstrophy, a measure of the rotational kinetic energy of the flow, is defined
locally as $\xi\equiv\frac{1}{2}\boldsymbol{\omega}\cdot\boldsymbol{\omega}$.
The rate of local enstrophy production is then calculated as

$\displaystyle\frac{\mathrm{d}}{\mathrm{d}t}\xi=\boldsymbol{\omega}\cdot\left(\frac{\mathrm{d}}{\mathrm{d}t}\boldsymbol{\omega}\right),$

which can be decomposed between resolved and unresolved domains as

$\displaystyle\frac{\mathrm{d}}{\mathrm{d}t}\xi=\boldsymbol{\omega}\cdot\left(\frac{\mathrm{d}}{\mathrm{d}t}\overline{\boldsymbol{\omega}}\right)+\boldsymbol{\omega}\cdot\left[\frac{\mathrm{d}}{\mathrm{d}t}\left(\boldsymbol{\omega}-\overline{\boldsymbol{\omega}}\right)\right],$

The local enstrophy production in the resolved domain is then defined as

$\displaystyle\frac{\mathrm{d}}{\mathrm{d}t}\xi_{r}\equiv\boldsymbol{\omega}\cdot\left(\frac{\mathrm{d}}{\mathrm{d}t}\overline{\boldsymbol{\omega}}\right),$

and the global enstrophy production in the resolved domain is then calculated
by integration,

$\displaystyle\int\limits_{-\infty}^{\infty}\frac{\mathrm{d}}{\mathrm{d}t}\xi_{r}\,\mathrm{d}\mathbf{y}=\int\limits_{-\infty}^{\infty}\boldsymbol{\omega}\cdot\left(\frac{\mathrm{d}}{\mathrm{d}t}\overline{\boldsymbol{\omega}}\right)\,\mathrm{d}\mathbf{y}.$

The SFS contribution to the rate of enstrophy production is isolated through
operator splitting as

$\displaystyle\left(\int\limits_{-\infty}^{\infty}\frac{\mathrm{d}}{\mathrm{d}t}\xi_{r}\,\mathrm{d}\mathbf{y}\right)_{\mathrm{SFS}}\equiv\int\limits_{-\infty}^{\infty}\boldsymbol{\omega}\cdot\left(\frac{\mathrm{d}}{\mathrm{d}t}\overline{\boldsymbol{\omega}}\right)_{\mathrm{SFS}}\,\mathrm{d}\mathbf{y}.$

Now, we constrain the SFS model to match the enstrophy production of the true
SFS contribution,

$\displaystyle\int\limits_{-\infty}^{\infty}\boldsymbol{\omega}\cdot\left(\frac{\mathrm{d}}{\mathrm{d}t}\overline{\boldsymbol{\omega}}\right)_{\mathrm{SFS\,model}}\,\mathrm{d}\mathbf{y}=\int\limits_{-\infty}^{\infty}\boldsymbol{\omega}\cdot\left(\frac{\mathrm{d}}{\mathrm{d}t}\overline{\boldsymbol{\omega}}\right)_{\mathrm{true\,SFS}}\,\mathrm{d}\mathbf{y}.$

Using the singular particle approximation,
$\boldsymbol{\omega}(\mathbf{x})\approx\sum\limits_{p}\boldsymbol{\Gamma}_{p}\delta(\mathbf{x}-\mathbf{x}_{p})$,
the integrals collapse into

$\displaystyle\sum\limits_{p}\boldsymbol{\Gamma}_{p}\cdot\left(\frac{\mathrm{d}}{\mathrm{d}t}\overline{\boldsymbol{\omega}}\right)_{(\mathbf{x}_{p})}^{\mathrm{SFS\,model}}=\sum\limits_{p}\boldsymbol{\Gamma}_{p}\cdot\left(\frac{\mathrm{d}}{\mathrm{d}t}\overline{\boldsymbol{\omega}}\right)_{(\mathbf{x}_{p})}^{\mathrm{true\,SFS}}.$

One instance that this equality is satisfied is when each term in the sum
satisfies

$\displaystyle\boldsymbol{\Gamma}_{p}\cdot\left(\frac{\mathrm{d}}{\mathrm{d}t}\overline{\boldsymbol{\omega}}\right)_{(\mathbf{x}_{p})}^{\mathrm{SFS\,model}}=\boldsymbol{\Gamma}_{p}\cdot\left(\frac{\mathrm{d}}{\mathrm{d}t}\overline{\boldsymbol{\omega}}\right)_{(\mathbf{x}_{p})}^{\mathrm{true\,SFS}}.$

Differentiating with respect to the filter width, we arrive to

$\displaystyle\boldsymbol{\Gamma}_{p}\cdot\frac{\partial}{\partial\sigma}\left(\frac{\mathrm{d}}{\mathrm{d}t}\overline{\boldsymbol{\omega}}\right)_{(\mathbf{x}_{p})}^{\mathrm{SFS\,model}}=\boldsymbol{\Gamma}_{p}\cdot\frac{\partial}{\partial\sigma}\left(\frac{\mathrm{d}}{\mathrm{d}t}\overline{\boldsymbol{\omega}}\right)_{(\mathbf{x}_{p})}^{\mathrm{true\,SFS}}.$
(48)

When SFS advection is neglected, the SFS model and the true SFS become

$\displaystyle\left(\frac{\mathrm{d}}{\mathrm{d}t}\overline{\boldsymbol{\omega}}\right)^{\mathrm{SFS\,model}}=-C_{d}\mathbf{E}_{\mathrm{str}}$
$\displaystyle\left(\frac{\mathrm{d}}{\mathrm{d}t}\overline{\boldsymbol{\omega}}\right)^{\mathrm{true\,SFS}}={\overline{(\boldsymbol{\omega}\cdot\nabla)\mathbf{u}}-(\overline{\boldsymbol{\omega}}\cdot\nabla)\overline{\mathbf{u}}},$

thus Eq. 48 is simply the dot-product of Eq. 46 with
$\boldsymbol{\Gamma}_{p}$,

$\displaystyle
C_{d}\boldsymbol{\Gamma}\cdot\frac{\partial\mathbf{E}_{\mathrm{str}}}{\partial\sigma}=-\boldsymbol{\Gamma}\cdot\frac{\partial}{\partial\sigma}\left[\overline{(\boldsymbol{\omega}\cdot\nabla)\mathbf{u}}-(\overline{\boldsymbol{\omega}}\cdot\nabla)\overline{\mathbf{u}}\right].$

Finally, following the same steps that led from Eq. 46 to Eq. 47, the
enstrophy balance in Eq. 48 becomes

$\displaystyle
C_{d}\boldsymbol{\Gamma}_{p}\cdot\mathbf{m}=\boldsymbol{\Gamma}_{p}\cdot\mathbf{L},$

and $C_{d}$ is then calculated dynamically as

$\displaystyle
C_{d}=\frac{\boldsymbol{\Gamma}_{p}\cdot\mathbf{L}}{\boldsymbol{\Gamma}_{p}\cdot\mathbf{m}}.$
(49)

Thus, this $C_{d}$ calculated at the position of every particle is the
coefficient that approximates the derivative balance while satisfying a local
balance of enstrophy production between the model and the true SFS
contribution.

####  V.3.3 Lagrangian Average

Similar to the classic procedure based on the Germano identity Germano  _et
al._ (1991), Eq. 49 poses numerical issues when the denominator is close to
zero, leading to large fluctuations. In the classic procedure, Meneveau et al.
Meneveau, Lund, and Cabot (1996) addressed this issue by integrating both
numerator and denominator along Lagrangian trajectories (pathlines),
effectively building ensemble averages. Applying this technique, our dynamic
procedure becomes

$\displaystyle
C_{d}=\frac{\left<\boldsymbol{\Gamma}_{p}\cdot\mathbf{L}\right>}{\left<\boldsymbol{\Gamma}_{p}\cdot\mathbf{m}\right>},$
(50)

where $\left<\cdot\right>$ denotes the Lagrangian integration. As suggested by
Meneveau et al., the integration is performed as a relaxation process at every
time step of the form

$\displaystyle\left<\phi\right>_{\mathrm{new}}=(1-\alpha)\left<\phi\right>_{\mathrm{old}}+\alpha\phi,$

where $\alpha$ is calculated as $\alpha=\Delta t/T\leq 1$, $\Delta t$ is the
time step of the simulation, and $T$ is the time length of the ensemble
average.

###  V.4 Backscatter Control

For stability reasons, it is common in LES to control the amount of energy
being backscattered from the unresolved scales into the resolved scales. The
SFS term is purely-dissipative if its contributions to the enstrophy budget
decrease the total enstrophy. This is,

$\displaystyle\left(\int\limits_{-\infty}^{\infty}\frac{\mathrm{d}}{\mathrm{d}t}\xi_{r}\,\mathrm{d}\mathbf{y}\right)_{\mathrm{SFS}}\leq
0.$

Following the derivation in Section V.3.2, this condition is satisfied if

$\displaystyle\boldsymbol{\Gamma}_{p}\cdot\left(\frac{\mathrm{d}}{\mathrm{d}t}\overline{\boldsymbol{\omega}}\right)_{(\mathbf{x}_{p})}^{\mathrm{SFS\,model}}\leq
0$

for each particle, or equivalently,

$\displaystyle
C_{d}\boldsymbol{\Gamma}_{p}\cdot\mathbf{E}_{\mathrm{str}}(\mathbf{x}_{p})\geq
0$

Therefore, enstrophy backscatter can be filtered out at each particle by
clipping the model coefficient as $C_{d}=0$ whenever the condition shown above
is not satisfied.

###  V.5 Usage in Conventional Mesh-Based CFD

Even though our SFS model is tailored for vortex methods as it only uses the
primitive variables of the VPM, the model can also be readily applied to
conventional mesh-based CFD as follows. The vortex strength is expressed from
mesh quantities as
$\boldsymbol{\Gamma}_{q}=\left<\boldsymbol{\omega}\right>_{q}\mathrm{Vol}_{q}$,
where $\left<\boldsymbol{\omega}\right>_{q}$ and $\mathrm{Vol}_{q}$ are the
average vorticity and volume associated to each element in the grid (i.e.,
cells in a finite volume method or nodes in a finite difference method). The
model in the vorticity transport equation then becomes

$\displaystyle\frac{\partial T_{ij}}{\partial x_{j}}\left(\mathbf{x}\right)=$
$\displaystyle
C_{d}\Delta^{-3}\sum\limits_{q}\left<\omega\right>_{j}^{q}\mathrm{Vol}_{q}\left(\frac{\partial\overline{u_{i}}}{\partial
x_{j}}\left(\mathbf{x}_{q}\right)-\frac{\partial\overline{u_{i}}}{\partial
x_{j}}\left(\mathbf{x}\right)\right)\zeta\left(\frac{\|\mathbf{x}-\mathbf{x}_{q}\|}{\Delta}\right),$

where the index $q$ iterates over each element in the grid, $\Delta$ is the
width of the grid filter, and $\zeta$ is the filter kernel. Furthermore, the
model can be used in the pressure-velocity form of the momentum equation
noticing that the SFS terms in the vorticity transport equation are simply the
curl of the SFS term in the pressure-velocity equation,
$\frac{\partial\tau_{ij}}{\partial x_{j}}$, as

$\displaystyle\epsilon_{ijk}\frac{\partial}{\partial
x_{j}}\frac{\partial\tau_{kl}}{\partial x_{l}}=\frac{\partial
T^{\prime}_{ij}}{\partial x_{j}}-\frac{\partial T_{ij}}{\partial x_{j}},$

where $\epsilon_{ijk}$ is the Levi-Civita tensor. Hence, the SFS models of the
vorticity transport equation can be used to model
$\frac{\partial\tau_{ij}}{\partial x_{j}}$ by “un-curling” the equation above.

##  VI Results

In the preceding sections we have developed a scheme for numerically solving
the LES-filtered Navier-Stokes equations in their vorticity form. The proposed
scheme uses a reformulation of the VPM and a novel model of SFS vortex
stretching to achieve a meshless large eddy simulation.

We now proceed to test and validate both the VPM reformulation and the vortex-
stretching SFS model that comprise our meshless LES. Simulations are compared
to results reported in the literature from experimental work, direct numerical
simulation (DNS), the lattice Boltzmann method (LBM), and unsteady Reynolds-
average Navier-Stokes (URANS) simulation. These validation cases also serve as
examples on how to impose initial and boundary conditions in our meshless
scheme.

In Section VI.1, advection and viscous diffusion is validated on an isolated
vortex ring, which is also used to test convergence of the reformulated VPM.
In Section VI.2, vortex stretching is validated through the simulation of two
leapfrogging vortex rings. Results with classic and reformulated VPM are first
compared without SFS effects. The SFS model is then introduced to capture the
subfilter energy cascade that leads the rings into turbulent breakdown. In
Section VI.3, a turbulent round jet is simulated to test the accuracy and
numerical stability of the LES scheme. The evolution of the jet and predicted
Reynolds stress are compared to experimental measurements, validating the
scheme as an LES method accurately resolving large-scale features of turbulent
flow. Finally, in Section VI.4, our meshless LES is compared to conventional
mesh-based CFD in a real-world engineering problem: an aircraft rotor and its
wake.

###  VI.1 Vortex Ring

A vortex ring travels with a self-induced velocity that slowly decelerates as
its vorticity spreads due to viscous diffusion. This poses a good validation
case for both vorticity advection and viscous diffusion. The initial vorticity
inside the ring core is assumed to follow a Gaussian distribution of the form
${\omega_{\theta}(r)=\frac{\Gamma_{0}}{\pi a^{2}}e^{-r^{2}/a^{2}}}$, where
$\Gamma_{0}$ and $a$ are the ring’s initial circulation and core size,
respectively. The initial Reynolds number is defined as
$\mathrm{Re}=\Gamma_{0}/\nu$. In all the following simulations, the rings were
initially discretized by placing particles (evenly spaced by $\Delta x$)
everywhere that the local vorticity was larger than 5% of the peak vorticity
$\omega_{\theta}(0)$, and then performing a radial-basis function (RBF) fit to
the target vorticity $\omega_{\theta}$ to determine the initial vortex
strength of the particles. The initial particle size $\sigma$ was set as to
provide a particle overlap $\lambda=\frac{\sigma}{\Delta x}$ of 2.4. Fig. 6
shows a discretized vortex ring along with a depiction of the ring’s geometric
parameters.

A vortex ring at $\mathrm{Re}=7500$ and thickness $\frac{a}{R}=0.2$ was
simulated and compared to direct-numerical simulation (DNS) reported by Archer
et al. Archer, Thomas, and Coleman (2008) The DNS is a finite-difference
method with second-order spatial accuracy solving the incompressible pressure-
velocity form of the Navier Stokes equations, using over 100 million grid
cells.

Figure 6: Vortex ring as discretized through particles with $\Delta x=0.13a$.
Figure 7: Spatial convergence of vortex ring simulated with reformulated VPM,
compared to direct numerical simulation (DNS).

Convergence of the reformulated VPM was tested without the SFS model by
increasing the spatial resolution from $\Delta x=0.17a$ (resulting in 60k
particles) to $\Delta x=0.08a$ (resulting in 500k particles), shown in Fig. 7.
To minimize numerical noise, we suppressed the core-resetting step of the
viscous diffusion scheme Barba, Leonard, and Allen (2005b), which leads the
VPM simulation to become increasingly coarser as time goes by. For this
reason, and without the SFS model, the VPM simulation is under-resolved and
inaccurate in the turbulent regime This is observed in Fig. 7 as an abrupt
turbulent breakdown after transition. In the laminar regime (nondimensional
time between 0 and 50), the velocity shows little sensitivity to a spatial
discretization finer than $\Delta x=0.13a$, being between 4% of the DNS
velocity at time 0, as shown in the zoom-in box of Fig. 7. The DNS shows a
slow transition from laminar to turbulent starting at nondimensional time 75,
while the VPM also transitions around time 75 for $\Delta x\leq 0.13a$,
quickly leading the ring into turbulent breakdown. Hence, $\Delta x=0.13a$ was
used in all the simulations that follow, which results in 150k particles.
These 150k vortex elements, contrasted with the 100 million grid elements used
in the DNS, highlight the computational efficiency of our meshless LES, which
enabled us to run our simulations in only a matter of hours on a desktop
workstation.

Fig. 8 compares the DNS to classic and reformulated VPM simulations, along
with the analytic velocity of a viscous Gaussian vortex ring (derived in
Appendix C). Archer et al. Archer, Thomas, and Coleman (2008) identified that
the initial Gaussian distribution goes through a transient phase until
achieving a steady skewed distribution. They denoted the time that the
simulation takes to achieve this steady distribution as $t^{*}$, reporting
$t^{*}=25$ for their DNS. Similarly, the VPM simulations were also observed to
have a phase of transient vorticity distribution. In Fig. 8 both the classic
VPM simulation and the analytic solution were shifted by $t^{*}=25$, showing
good agreement with the DNS velocity at time 0. However, the classic VPM
slowly drifts away from both the analytic and DNS velocity after time 0. In
contrast, the reformulated VPM starts off with a velocity between 4% of the
DNS velocity at time 0 and maintains good agreement up to turbulent breakdown.
This confirms that the reformulated VPM accurately resolves vorticity
advection and viscous diffusion in the laminar regime.

Figure 8: Vortex ring simulated with classic VPM ($t^{*}=25$) and reformulated
VPM ($t^{*}=0$), compared to direct numerical simulation (DNS, $t^{*}=25$) and
analytic solution ($t^{*}=25$).

###  VI.2 Leapfrogging Vortex Rings

Two vortex rings that travel in the same direction can repeatedly leapfrog,
stretching and contracting as they pass one another. This poses a good
validation case for vortex stretching. Two identical vortex rings with initial
radii $R_{0}=1$, thickness $\frac{a}{R_{0}}=0.1$, and separation $\Delta
Z=R_{0}$ were simulated at $\mathrm{Re}=3000$ and compared to the simulation
through lattice Boltzmann method (LBM) reported by Cheng et al. Cheng, Lou,
and Lim (2015) The LBM simulation is a finite-difference DNS with second order
spatial accuracy solving the generalized lattice-Boltzmann equation,
reportedly using 15 billion cells.Cheng, Lou, and Luo (2010) In contrast, the
VPM simulations used only 600k particles, which again highlights the
computational efficiency of our meshless LES.

First, the simulation was performed on both classic and reformulated VPM
without the SFS model, shown in Fig. 9. Reformulated VPM and LBM show good
agreement in the first leapfrog ($Z/R_{0}\leq 3$), but they deviate in
subsequent cycles ($Z/R_{0}>3$) once turbulence starts to arise in the LBM
causing the rings to slow down and run into each other. This is not captured
in the reformulated VPM as the simulation is under-resolved (it does not
resolve all the scales of turbulent motion) and leads to a laminar behavior,
which indicates the need for an SFS model. Before incorporating the SFS model,
however, we conclude that the reformulated VPM accurately resolves vortex
stretching in the laminar regime.

Figure 9: Leapfrogging vortex rings simulated with classic VPM (cVPM) and
reformulated VPM (rVPM) without SFS model, compared to lattice Boltzmann
method (LBM) and analytic solution.

Fig. 9 also shows the classic VPM simulation and compares it to an analytic
solution derived in Appendix C. The classic VPM leads to dynamics that are
substantially different from both reformulated VPM and LBM, but it matches the
analytic solution. As described in Appendix C, the analytic solution assumes
that (1) the core size $a$ is not affected by the stretching of the ring and
that (2) the core size does not affect the interactions between rings. The
former clearly violates conservation of mass, while the latter is suspected to
violate conservation of angular momentum. Such good agreement between the
classic VPM and this unphysical analytic solution supports our previous claim
in Section IV.1.3 that the classic VPM may lead to unphysical results, which
may also cause the method to be numerically unstable.

Figure 10: Leapfrogging of vortex rings using reformulated VPM with and
without SFS model, compared to lattice Boltzmann method (LBM). Visualization
of rVPM+SFS simulation overlaid on top.

Figure 11: Volume rendering of vorticity field as the leapfrogging rings
merge. Particles shown as black points. Video available at
https://youtu.be/viE-i0KzVOw in the preprint version of this paper.

Finally, the leapfrog simulation was repeated with the SFS model added to the
reformulated VPM, shown in Fig. 10. Notice that the SFS model only affects the
dynamics after the first leapfrog ($Z/R_{0}>3$) once turbulence starts to
develop. This confirms that the dynamic procedure succeeds at auto-regulating
the SFS model, suppressing it in the laminar regime while activating it in the
turbulent regime. The diffusivity of the SFS model slows the rings down,
eventually running into each other and mixing after three leapfrogs at
$Z/R_{0}>6$ (shown in more detail in Fig. 11). As seen in Fig. 10, the process
of slowing down and running into each other happens at a different pace
between the rVPM+SFS simulation and the LBM. It is unclear whether this
discrepancy is due to our vortex-stretching SFS model lacking the full range
of turbulent diffusion (requiring the addition of a vorticity-advection SFS
model), or due to excessive numerical dissipation in the LBM. However, they
are both consistent in that the rings mix after three leapfrog cycles, showing
the ability of the SFS model to capture diffusion associated with vortex
stretching in the turbulent regime.

In order to estimate the computational cost of the VPM reformulation and the
SFS model, Table 2 compares the computational time of each simulation to the
classic VPM. The VPM reformulation runs as fast as the classic VPM, adding no
significant overhead. The SFS model adds a computational overhead of 8% when a
constant model coefficient is prescribed and 43% with the dynamic procedure.
Comparing this to a benchmark study by Chapelier et al. Chapelier, Wasistho,
and Scalo (2018), we conclude that the overhead of our SFS model is comparable
to a Smagorinsky eddy-viscosity model, while the dynamic procedure adds an
overhead comparable to a Germano-identity dynamic model. The work of Chapelier
et al. suggests that further speedup could be achieved by implementing a
sensor function, but this is left for future work.

Table 2: Computational time of the leapfrog VPM simulations. Formulation |  SFS Model |  CPU Time |  Overhead
---|---|---|---
Classic | None | $t_{\mathrm{ref}}$ | –
Reformulated | None | $1.01\,t_{\mathrm{ref}}$ |  $<$1%
Reformulated | Constant $C_{d}$ | $1.08\,t_{\mathrm{ref}}$ | +8%
Reformulated | Dynamic $C_{d}$ | $1.43\,t_{\mathrm{ref}}$ | +43%

###  VI.3 Turbulent Round Jet

Figure 12: Boundary condition of round jet simulation defined with particles.
Figure 13: Streamwise (top) and radial (bottom) velocity close to nozzle exit
in experiment by Quinn and Militzer Quinn and Militzer (1989) (markers),
compared to boundary condition of simulation (solid line) probed at $z/d=0.1$.
Figure 14: Global enstrophy of turbulent round jet with classic VPM (cVPM) and
reformulated VPM (rVPM).

A jet discharging into a quiescent environment poses a canonical case for the
study of turbulence, encompassing a laminar region near the nozzle that breaks
down into turbulence away from the nozzle. Experimental measurements on a
round jet were used to test that the reformulated VPM is able to resolve the
mean and fluctuating components of turbulent flow, while resolving Reynolds
stress directly. The simulation replicated the experiment by Quinn and
Militzer Quinn and Militzer (1989) which used a contoured nozzle with an exit
diameter $d$ of ${45.4\,\mathrm{mm}}$, discharging air into stagnant ambient
air with a centerline velocity $U_{c}$ of ${60\,\mathrm{m/s}}$. This
corresponds to a Mach number of 0.18 and a diameter-based Reynolds number
$\mathrm{Re}=U_{c}d/\nu$ of ${2\times 10^{5}}$.

The simulation assumed a top-hat velocity profile with smooth edges at the
nozzle exit given by

$\displaystyle u_{e}(r)=U_{c}\tanh\left(\frac{\frac{d}{2}-r}{\theta}\right),$

where $\theta$ is the momentum thickness of the shear layer, assumed to be
$\theta=0.025d$. As depicted in Fig. 12, this velocity profile was imposed as
a boundary condition at $z=0$ by defining a volumetric cylinder of length $h$
with particles spanning from $z=-h/2$ to $z=h/2$ and performing an RBF fit to
the vorticity profile $\omega_{\theta}=\frac{\mathrm{d}u_{e}}{\mathrm{d}r}$.
This computes the vortex strengths that induce the velocity profile $u_{e}$
inside the cylinder. The $+z$ half of the cylinder was then removed, the $-z$
half was kept in the computational domain as static particles throughout the
simulation, while the set of particles computed at $z=0$ were injected at each
time step as free particles. The resulting boundary condition showed
sensitivity to the cylinder length $h$, leading to a straight jet when $h\geq
2d$, a contracting jet when $h<2d$, and an actuator disk when $h\rightarrow
0$. Since Quinn and Militzer reported a jet that was slightly contracting, the
length of the boundary-condition cylinder was tailored to match the streamwise
and radial velocity that they measured at the nozzle exit plane shown in Fig.
13, finding sufficient agreement when $h=1.7d$.

Figure 15: Vortex particles (top), vortex strength (middle), and SFS model
coefficient (bottom) close to jet nozzle at $t=48\,\mathrm{ms}$. Figure 16:
Volume rendering of vorticity field in turbulent jet at $t=48\,\mathrm{ms}$
showing distinct flow features: Coherent structures form in the initial region
($z<1d$) that mix and break down by $z>3d$.

$t={\color[rgb]{1,1,1}\definecolor[named]{pgfstrokecolor}{rgb}{1,1,1}\pgfsys@color@gray@stroke{1}\pgfsys@color@gray@fill{1}{0}}2\,\mathrm{ms}$

$t={\color[rgb]{1,1,1}\definecolor[named]{pgfstrokecolor}{rgb}{1,1,1}\pgfsys@color@gray@stroke{1}\pgfsys@color@gray@fill{1}{0}}4\,\mathrm{ms}$

$t={\color[rgb]{1,1,1}\definecolor[named]{pgfstrokecolor}{rgb}{1,1,1}\pgfsys@color@gray@stroke{1}\pgfsys@color@gray@fill{1}{0}}8\,\mathrm{ms}$

$t=10\,\mathrm{ms}$

Figure 17: Evolution of turbulent jet simulation: (left) volume rendering and
(right) slice of vorticity field.

$t=22\,\mathrm{ms}$

$t=45\,\mathrm{ms}$

$t=50\,\mathrm{ms}$

Figure 18: Evolution of turbulent jet simulation (continued).

The simulation was first attempted with the classic VPM, but it quickly ended
in numerical blow up at the initial stage of the jet. This is shown in Fig. 14
through an abrupt jump in global enstrophy. The reformulated VPM proved to be
numerically stable in the initial and transition stages of the jet, however,
the simulation becomes unstable in the fully-developed turbulent regime as
enstrophy builds up in the absence of SFS turbulent diffusion. Introducing the
SFS model, the enstrophy production of the inlet balances out with the
forward-scatter of the fully-developed region and the rVPM simulation becomes
indefinitely stable, as shown in Fig. 14. Interestingly, the SFS model did not
stabilize the classic VPM simulation, which indicates that both the
reformulated VPM and the SFS model are needed to achieve numerical stability.

Figure 19: Average SFS model coefficient in turbulent jet simulation as
computed through the dynamic procedure. Figure 20: Development of streamwise
mean velocity profile along jet. Simulation: - rVPM. Experimental: $\star$
Quinn and Militzer Quinn and Militzer (1989); $\oplus$ Fellouah and Pollard
Fellouah and Pollard (2009); $\triangleright$ Mi et al., contraction nozzle
Mi, Nobes, and Nathan (2001); $\ast$ Iqbal and Thomas Iqbal and Thomas (2007);
$\diamond$ Romano Romano (2002); $\boxdot$ Xu and Antonia, contraction nozzle
Xu and Antonia (2002) . Figure 21: Development of fluctuating component of the
streamwise velocity along jet. Simulation: - rVPM. Experimental: $\star$ Quinn
and Militzer Quinn and Militzer (1989); $\oplus$ Fellouah and Pollard Fellouah
and Pollard (2009); $\ast$ Iqbal and Thomas Iqbal and Thomas (2007); $\boxdot$
Xu and Antonia, contraction nozzle Xu and Antonia (2002) . Figure 22:
Development of Reynolds stress along jet. Simulation: - rVPM. Experimental:
$\oplus$ Fellouah and Pollard Fellouah and Pollard (2009); $\boxdot$ Xu and
Antonia, contraction nozzle Xu and Antonia (2002) .

Once proven stable, the simulation was run for $50$ ms with time step $\Delta
t=0.02\,\mathrm{ms}$ until achieving a fully-developed region that was
statistically stationary. This was run on a single node with 128 CPU cores
(dual 64-core AMD EPYC 7702, 2.0 GHz) resulting on a wall-clock time of two
and a half days using up to $3\times 10^{6}$ particles. Fig. 15 shows the
particle field, vortex strengths, and the SFS model coefficient $C_{d}$ close
to the nozzle at $t=48\,\mathrm{ms}$, while Fig. 16 shows the vorticity of the
entire field. As seen in Fig. 16, coherent vortical structures form in the
initial region ($z<1d$), which leapfrog and mix transitioning to fully-
developed turbulent flow by $z>3d$. Notice in Fig. 15 that the model
coefficient is negligibly small in the initial region and increases to
$0.1<C_{d}<1$ in the fully-developed regime. This shows that the dynamic
procedure succeeds at automatically calibrating $C_{d}$, hindering SFS
turbulent diffusion in the laminar regime while facilitating it in turbulent
regions.

The formation and time evolution of the jet are shown in Fig. 18 and available
as a video in the supplemental content††††††Available at
https://youtu.be/V9hthE7m1d4 in the preprint version of this paper..
Initially, a vortex ring forms at the head of the jet. As the ring travels
downstream, a vortex sheet is deployed forming the shear layer, seen at
$t=2\,\mathrm{ms}$. The vortex sheet is stretched as the flow develops,
eventually rolling up and forming filaments, seen at $t=4\,\mathrm{ms}$. The
filaments pair up, leapfrog, merge, and breakdown (seen at $t=8\,\mathrm{ms}$)
eventually forming the fully-developed turbulent region at
$t>20\,\mathrm{ms}$. Fig. 19 shows the history of $C_{d}$ throughout this
process, calculated at each time step as the average $|C_{d}|$ over all the
particles where $C_{d}\neq 0$. Here we see that the average $C_{d}$ is
negligibly small in the initial development of the jet, but it ramps up and
converges to a value of 0.43 as the flow becomes fully developed, again
confirming the ability of the dynamic procedure to automatically calibrate the
SFS model at each flow regime.

To validate the dynamics predicted by the simulation, the velocity was probed
at five stations along the jet and statistical properties were compared to the
experimental measurements reported by Quinn and Militzer Quinn and Militzer
(1989). The data was also supplemented with other experiments of similar round
jets compiled by Ball et al. Ball, Fellouah, and Pollard (2012) Statistical
properties were calculated through temporal and spatial ensemble averages
after the flow became statistically stationary in the region $0\leq z\leq 5d$.
Averaging was performed over the time interval $40\,\mathrm{ms}\leq
t<50\,\mathrm{ms}$.

Fig. 20 shows the mean component of the streamwise velocity profile, $u_{z}$,
normalized by the centerline velocity at each station. The mean component
shows excellent agreement with the experiments in the initial and transition
regions ($0\leq z\leq 3d$) and reasonable agreement in the fully-developed
turbulent region ($z\geq 4d$) though slightly overexpanded. Fig. 21 shows the
fluctuating component of the streamwise velocity, $u^{\prime}_{z}$, or
standard deviation, while Fig. 22 shows the Reynolds stress between streamwise
and radial velocity, defined as the square root of the covariance,
$\sqrt{u^{\prime}_{z}u^{\prime}_{r}}$. At $z=0$, the simulation shows only a
small fluctuation and Reynolds stress as the flow is dominated by the boundary
condition. Away from the nozzle exit plane ($z>0$), fluctuations and Reynolds
stress concentrate at the shear layer ($r/R=1$) in the laminar region
($z\approx 1d$), and gradually spread as the jet breaks down in the turbulent
regime ($z\geq 4d$). Fluctuations are overpredicted (and, as a consequence,
also the Reynolds stress), however this is in within reasonable agreement with
the experiments. These predictions can be further improved in future work with
the addition of an SFS model of vorticity advection increasing the turbulent
diffusion (and damping out fluctuations), or implementing a spatial adaptation
strategy to better resolve small scales. However, the current agreement with
the experiments suffices to confirm that our scheme is an LES able to resolve
mean and fluctuating large-scale features of the flow. Also, its ability to
directly resolve Reynolds stress renders the reformulated VPM a higher
fidelity approach than Reynolds-average approaches like RANS and URANS (where
Reynolds stresses are rather modeled), while being completely meshless.

###  VI.4 Aircraft Rotor in Hover

The rotation of blades in static air drives a strong axial flow caused by the
shedding of tip vortices. This is a challenging case to simulate since, in the
absence of a freestream, the wake quickly becomes fully turbulent as tip
vortices leapfrog and mix close to the rotor. Thus, a rotor in hover is a good
engineering application to showcase the accuracy, numerical stability, and
computational efficiency of the reformulated VPM.

In this test case, we simulated the experiment by Zawodny et al. Zawodny,
Boyd, Jr., and Burley (2016) consisting of a DJI 9443 rotor in hover at 5400
RPM. This two-bladed rotor is 9.4 inches in diameter, resulting in a tip Mach
number of 0.20 and chord and diameter-based Reynolds numbers at 70% of the
blade span of $6\times 10^{4}$ and $7\times 10^{5}$, respectively. Our
simulations were also compared to unsteady Reynolds-average Navier-Stokes
(URANS) results reported by Schenk, Schenk (2020) obtained with the commercial
software STAR-CCM+. The URANS is an unsteady compressible solver with an SST
$k$–$\omega$ turbulence model, resolving the blades with an all-$y^{+}$ wall
treatment on a rotating mesh surrounding the rotor. It used an unstructured
mesh with 14 million cells, mesh refinement down to a $y^{+}$ of only 30, and
time steps equivalent to $3^{\circ}$ of rotation on a first-order time
integration scheme, which is a rather coarse simulation, but requiring very
low computational resources.

The rotating blades are computed in our VPM through an actuator line model