Lagrangian particle tracking

Lagrangian particle tracking (LPT) is a method used in fluid mechanics to analyze the motion of particles when subjected to a flow field. It provides a Lagrangian perspective, in which the flow is described by tracking 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 introduced 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 large scales makes it suitable for studying the dispersion of atmospheric pollutants. 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]

Experimental field

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]

Particles behavior at low and high Stokes numbers

Extracting data about the flow field is done by using particle tracking velocimetry (PTV) methods. PTV methods use a lower particle concentration compared to particle image velocimetry (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) that are 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 are 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), illumination technology (High frequency lasers and scalable LED illumination^[13]), calibration (Volume self-Calibration method ^[14] ) and in post-processing algorithms (Shake-the-Box^[15] and iterative particle reconstruction techniques ^[16] ) make 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]

Example of iterative particle reconstruction technique^[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.

Computational field

In the field of computational fluid dynamics, Lagrangian particle tracking is typically employed in coupled CFD-DEM simulations.^[4] Unlike 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]

Uniform grid

For simple Cartesian grids, the tracking algorithm is straightforward, as the problem is decoupled along the three directions. 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]

Rectilinear grid

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 directions 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 grid

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]

Visual representation of a spatial transformation which transforms a curvilinear structured grid in a rectilinear structured 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 grid

Example of particle path in an unstructured grid^[19]

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 is effectively inside a cell.^[19] Using cells of simple shapes, such as triangles or convex quadrilaterals 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 that have curved faces, wall reflection, algorithm robustness checks, and improved efficiency by limiting the cell search only to neighboring ones.^{[19][20][21][22][23]}

See also

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