Lagrangian particle tracking

Lagrangian particle tracking (LPT) is a method used in fluid mechanics to analyse particles' motion when subjected to a flow field. It provides a Lagrangian perspective, in which the flow is described by following fluid parcels or tracers over time, rather than observing changes at fixed locations as in the Eulerian frame^[1].

In experimental studies, LPT is typically performed using three-dimensional Particle Tracking Velocimetry (3D-PTV)^[1]. Neutrally buoyant tracer particles are seeded into the flow, and their positions are recorded using high-speed cameras and stereo reconstruction techniques. The resulting particle paths allow for the study of turbulent structures, transport phenomena, and time-resolved Lagrangian statistics^[2][3].

In computational fluid dynamics, LPT refers to the numerical simulation of discrete particles embedded in a continuous flow field. The fluid phase is typically solved in an Eulerian framework, while the particle phase is resolved using Lagrangian mechanics^[4]. This approach, also termed Discrete Particle Simulation (DPS), is particularly suited to dilute multiphase flows, such as aerosols, particle deposition in the human airways^[5] and environmental particle transport^[6], where particle–fluid coupling is weak. Applications of LPT also include cases where coupling is not negligible, and a more advanced method such as the discrete element method (DEM) is needed. Examples of this case are industrial mixing^[7], combustion modelling^[8], sprays and fluidized beds^[9].

Beyond engineering and turbulence research, LPT has been widely adopted in environmental modelling. Its capacity to resolve particle motion over complex terrain and at large scales makes it suitable for studying atmospheric pollutant dispersion. In regional air quality assessments, LPT methods have been used for both forward simulations (predicting particle transport from known sources) and inverse modelling (inferring sources from observed concentrations). These techniques have proven effective in identifying transboundary pollution pathways and assessing exposure risks.^[10]

The main goal of Lagrangian particle tracking in the experimental field is to extract data about the flow field, such as the flow velocity, acceleration (the material derivative) and pressure fields in the Lagrangian frame^[1].

This is done by using particle tracking velocimetry (PTV) methods. Differently from particle image velocimetry (PIV), PTV methods use a lower particle concentration with respect to PIV methods, allowing to follow individual particles in the flow field. Typical particle image densities for PTV methods are between ${\displaystyle \sim 0.005}$ and ${\displaystyle \sim 0.02}$ ppp (particles per pixel) while PIV experiments can be performed with higher particle densities^[1].

In both PTV and PIV methods, the flow is seeded with neutrally buoyant particles (typically fluid ones in gaseous media and solid ones in fluid media) sufficiently small to follow the flow streamlines. This is verified when the Stokes number is sufficiently small, the typical condition being: ${\displaystyle St<0.1}$^[4].

The particles are illuminated twice at given intervals and a high definition camera is used to capture an image each time particles are illuminated. Each pair of images can be elaborated to extract the velocity field. Images are then post-processed to obtain the particles' velocity and positions^[11].

Several PTV methods have been proposed and their particular setup changes based on the specific needs of the experiment. The number of cameras used varies between one for the standard 2D-PTV method, up to 6 when state-of-the-art setups are used ^[12]. Recent advancement in camera (high frequency CMOS and CCD sensors) and illumination technology (High frequency lasers and scalable LED illumination^[13]), calibration (Volume self-Calibration method ^[14] ) as well as in post-processing algorithms (Shake-the-Box^[15] and iterative particle reconstruction techniques ^[16] ) made it possible to develop complex setups like time resolved 3D-PTV with particle image densities which reach PIV-level density (${\displaystyle \sim 0.2}$ ppp)^[1].

When performing PTV experiments, typical difficulties include^[1]:

• irregular scattering behavior caused by particle's different sizes, shape, rotation velocity and position with respect to the illumination source
• background intensity and light reflection caused by the model which can diminish the signal-to-noise ratio (SNR)
• image distortions due to improper calibration or camera astigmatism which can lead to non-Gaussian particle imaging, deteriorating the quality of post-processing results.

In the field of computational fluid dynamics, Lagrangian particle tracking is typically employed in coupled CFD-DEM simulations^[4]. Differently from what happens in experiments, during a CFD simulation the velocity of the particle is already known, as is its position at a given time instant ${\displaystyle t^{n}}$, and the velocity of the particle can be interpolated from the velocity of the fluid cell where the particle lies. Lagrangian particle tracking is instead needed to forecast in which cell the particle will end at the new time instant ${\displaystyle t^{n+1}}$. Specifically, reconstructing the particle's trajectory involves the numerical integration of the ordinary differential equation:

${\displaystyle {\frac {d\mathbf {x} }{dt}}=\mathbf {v} (\mathbf {x} )}$

with initial condition ${\displaystyle \mathbf {x} (t_{n})=\mathbf {x} _{p}^{n}}$. Assuming the equation is integrated in time with an explicit integration scheme, the particle velocity can be interpolated from the velocity of the cell where the particle is and the equation can be integrated, returning ${\displaystyle \mathbf {x} _{p}^{n+1}}$. The new position of the particle ${\displaystyle \mathbf {x} _{p}^{n+1}}$ is then used at the new timestep to advance the equation in time again. In order to do that, a tracking algorithm is needed to find the (possibly new) cell where the particle has moved at ${\displaystyle t^{n+1}}$^[17].

For simple Cartesian grids, the tracking algorithm is straightforward, as the problem is decoupled along the three directions and the grid coordinates can be directly computed knowing the cell index. Knowing the grid first point in space ${\displaystyle \mathbf {x_{0}} =(x_{0},y_{0},z_{0})}$ and the grid spacing ${\displaystyle \mathbf {\Delta } =(\Delta _{x},\Delta _{y},\Delta _{z})}$ then the coordinate of the center of cell ${\displaystyle (i,j,k)}$ is uniquely identified as:

${\displaystyle {\begin{cases}x_{i}=x_{0}+i\cdot \Delta _{x}\qquad \qquad {\text{for }}i=0:n_{x}-1\\y_{i}=y_{0}+j\cdot \Delta _{y}\qquad \qquad {\text{for }}j=0:n_{y}-1\\z_{i}=z_{0}+k\cdot \Delta _{z}\qquad \qquad {\text{for }}k=0:n_{z}-1\end{cases}}}$

The index of the new surrounding cell of the particle at time ${\displaystyle t^{n+1}}$ can be computed as:

${\displaystyle {\begin{cases}i={\Big \lfloor }{\frac {x_{p}^{n+1}-x_{0}}{\Delta _{x}}}{\Big \rceil }\\j={\Big \lfloor }{\frac {y_{p}^{n+1}-y_{0}}{\Delta _{y}}}{\Big \rceil }\\k={\Big \lfloor }{\frac {z_{p}^{n+1}-z_{0}}{\Delta _{z}}}{\Big \rceil }\end{cases}}}$

where ${\displaystyle \lfloor \cdot \rceil }$ is the rounding operator^[17].

In rectilinear grids, the coordinates of mesh elements at a given index can't be computed directly and has to be stored explicitly. In rectilinear grids, the coordinates in a given direction only depend on the index in the same directions:

${\displaystyle x_{i}=f(i,j,k)=f_{x}(i)\quad {\text{for }}i=0:n_{x}-1}$

This allows to decouple the problem in the three direction as for the Cartesian grid. Assuming that ${\displaystyle f_{x}(i)}$ is a monotonically increasing function, by saving the values ${\displaystyle f_{x}(i)}$ in an array ${\displaystyle X(i)}$ and knowing ${\displaystyle x_{p}^{n+1}}$ the cell index is efficiently found performing a binary search algorithm, whose cost is ${\displaystyle O(\log {n})}$^[18].

Curvilinear grids have the same topological structure of rectilinear grids (i.e. they are "structured") but the location of a cell center in one direction depends on all the three indices of the cell:

${\displaystyle {\begin{cases}x_{i,j,k}=f(i,j,k)\quad {\text{for }}i=0:n_{x}-1\\y_{i,j,k}=g(i,j,k)\quad {\text{for }}j=0:n_{y}-1\\z_{i,j,k}=h(i,j,k)\quad {\text{for }}k=0:n_{z}-1\end{cases}}}$

As interpolations and point location algorithms are more complex in curvilinear grids, a common procedure is to transform the physical space ${\displaystyle (P)}$ in a computational space ${\displaystyle (C)}$ where the curvilinear grid is mapped to a Cartesian grid^[17].

By doing so, the problem falls back to the rectilinear case. To solve the problem in the ${\displaystyle C}$-space, the following equation must be solved:

${\displaystyle {\frac {d{\boldsymbol {\xi }}}{dt}}={\boldsymbol {u}}({\boldsymbol {\xi }})}$

where ${\displaystyle {\boldsymbol {\xi }}=(\xi ,\eta ,\zeta )}$ denotes a position in the computational space. The velocity in the ${\displaystyle C}$-space is computed by transforming the velocity in the ${\displaystyle P}$-space:

${\displaystyle \mathbf {u} =\mathbf {J} ^{-1}\cdot \mathbf {v} }$

where ${\displaystyle \mathbf {J} }$ is the Jacobian matrix of the transformation that maps the two spaces:

${\displaystyle \mathbf {J} ={\begin{bmatrix}x_{\xi }&x_{\eta }&x_{\zeta }\\y_{\xi }&y_{\eta }&y_{\zeta }\\z_{\xi }&z_{\eta }&z_{\zeta }\end{bmatrix}}}$

where ${\displaystyle x_{\xi }}$ stands for ${\displaystyle {\frac {dx}{d\xi }}}$. Typically the velocity is not saved to memory in both spaces but it's just transformed to ${\displaystyle C}$-space when needed to avoid excessive memory usage^[17]. The Jacobian is computed locally using finite differences, and the number of different Jacobians used for each cell (e.g. one for each cell or one/multiple for each cell node) significantly affects the operation's precision ^[17].

Unstructured grids pose significant challenges to LPT algorithms, as the grid topology is irregular and needs to be explicitly stored at runtime. Moreover, unstructured grids can be composed of cells of arbitrary shapes, making it complex to develop an efficient location algorithm to check if a position effectively is inside a cell^[19]. Using cells of simple shapes such as triangles or convex quadrilateral in 2D and tetrahedra in 3D simplifies the checks on point locations. In this case, the cells are convex and their faces are planar surfaces, therefore a simplified in-cell-check can be performed by verifying that:

${\displaystyle (\mathbf {r} _{c}-\mathbf {r} _{p})\cdot \mathbf {n} \geqslant 0}$

where ${\displaystyle \mathbf {r} _{c}}$ is the face centroid and ${\displaystyle \mathbf {n} }$ is the outward normal of the face^[19].

Once the new particle position has been computed, the intersection between the particle trajectory and the cell's faces must be computed to deduce which cell the particle is moving to. If a particle's trajectory crosses a face, the particle is assigned to the neighboring cell. During a single time step, the particle can cross multiple cell faces, therefore this step can be iterated several times^[19]. Several methods have been proposed to effectively implement Lagrangian particle tracking in unstructured grids. Refined methods include cases with cells with curved faces, wall reflection, algorithm robustness checks and improved efficiency by limiting the cell search only to neighboring ones^{[19][20][21][22][23]}.

• Computational fluid dynamics
• Lagrangian and Eulerian specification of the flow field
• Particle image velocimetry
• Regular grid
• Unstructured grid

• Norouzi, H. R.; Zarghami, R.; Sotudeh-Gharebagh, R.; Mostoufi, N. (2016). Coupled CFD‐DEM Modeling. John Wiley & Sons, Ltd. ISBN 9781119005315.
• Schröder, A.; Schanz, D. (2023-01-19). "3D Lagrangian Particle Tracking in Fluid Mechanics". Annual Review of Fluid Mechanics. 55 (1): 511–540. doi:10.1146/annurev-fluid-031822-041721.
• Springer handbook of experimental fluid mechanics. Tropea, C.; Yarin, A. L.; Foss, J.F. (eds.). Berlin: Springer. 2007. ISBN 978-3-540-25141-5.{{cite book}}: CS1 maint: others (link)

• openPTV, an OSS PTV post-processing library
• TracTrac, an OSS PTV post-processing library developed by University of Rennes
• CFDEM, an OSS CFD-DEM software
• MFIX-DEM, an OSS CFD-DEM software developed by NETL