<figure><img src="image_1.png" /></figure>

HAL
open science

Generic energy function to predict particle positions in
4D-PTV

Ali Rahimi Khojasteh, Dominique Heitz

► To cite this version:

Ali Rahimi Khojasteh, Dominique Heitz. Generic energy function to predict particle positions in 4D-PTV. 20th International Symposium on Application of Laser and Imaging Techniques to Fluid Mechanics, Jul 2022, Lisbon, Portugal. 7 p. hal-03864727

HAL Id: hal-03864727

https://hal.inrae.fr/hal-03864727

Submitted on 21 Nov 2022

**HAL** is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

Public Domain
# Generic energy function to predict particle positions in 4D-PTV

AR. Khojasteh¹*, D. Heitz¹

1:INRAE, OPAALE, 17 avenue de Cucillé, 35044, Rennes, France

*Corresponding author: ali.rahimi-khojasteh@inrae.fr

**Keywords:** Particle tracking, Position estimation, Coherent motions

## ABSTRACT

Recent developments in time-resolved Particle Tracking Velocimetry (4D-PTV) consistently improved tracking accuracy and robustness by introducing predictive algorithms such as Shake-The-Box (STB). We propose a generic non-dimensional energy function to predict particle positions by leveraging additional information from coherent neighbour motions. The proposed energy function is a combination of position history, neighbour velocity and acceleration terms. Coherent neighbours are quantified by locally utilising Finite Time Lyapunov Exponent (FTLE) as an objective Lagrangian Coherent Structure (LCS) diagnostic method. Synthetic analysis of the wake behind a smooth cylinder at Reynolds number equal to 3900 showed that the optimal solution of the minimised energy function could be modelled as a function of the measurement uncertainties. The model was assessed with the 2D homogeneous isotropic turbulent flow (HIT) as well and was found to be case and Reynolds independent. Results of the 4D-PTV experimental study of the same wake flow configurations are reported. We compared predicted positions with the optimised final positions of STB. It was found that the generic energy function succeeded in estimating particle positions with minimum deviation to the optimised positions in comparison with other techniques.

## 1. Generic energy function minimisation

This paper discusses a parametric approach in time-resolved Particle Tracking Velocimetry (4D-PTV) to obtain a generic particle position predictor function. In the classic 4D-PTV algorithms such as Shake-The-Box (STB) (Schanz et al., 2016), predicted particle positions are given to the optimisation process for further corrections. The optimisation can deal with slight deviations between the predicted position and the true position. However, the optimisation fails to find the true position if the deviation is large enough to have multi-candidates for a single particle at $t_{n+1}$. This shows the importance of having an appropriate prediction applicable in different flow complexities. The proposed idea is initiated by arguing that predictions in PTV techniques such as polynomial and Wiener filter (Schanz et al., 2016) focus on a single particle individually, while this single particle is not acting alone (Khojasteh, Yang, et al., 2021). The present study was designed as a complementary function for 4D-PTV algorithms such as STB (Schanz et al., 2016) and KLPT

<figure><img src="image_2.png" /></figure>
<figure><img src="image_3.png" /><figcaption>Figure 1. Contour plots of the estimation error with respect to $\alpha'_1$ and $\alpha'_2$ weighting terms. The synthetic trajectories were generated based on the 1st LPT challenge estimation errors. a) the 3D wake behind a cylinder at a Reynolds number equal to 3900, b) 2D homogeneous isotropic turbulent flow at a Reynolds number equal to 3000</figcaption></figure>

(Yang & Heitz, 2021). Khojasteh, Heitz, et al. (2021) showed that an energy function obtained from three terms of position, velocity, and acceleration could significantly reduce the bias error caused by the flow complexities compared with polynomial and Wiener filter predictors. However, an estimated position in a real experiment needs to be optimised for both bias error and measurement uncertainties. The derived energy function suggested by Khojasteh, Heitz, et al. (2021) has dimensional terms (position, velocity, and acceleration) with equal weights to predict the particle positions. In the present study, we defined weighting parameters to the energy function to reach the minimum estimation bias error. The weighted energy function can be written as,

$$J = \underbrace{\left(\frac{1}{n}\sum_{i=1}^{n}(X_i - y_i)^2\right)}_{\text{History}} + \underbrace{\alpha_1(\dot{X}_n - \ddot{y}_n)^2}_{\text{Coherent velocity}} + \underbrace{\alpha_2(\ddot{X}_n - \ddot{y}_n)^2}_{\text{Coherent acceleration}}, \quad (1)$$

where $X$ is the estimated positions based on known $y$ observations in time. $\alpha_1$ and $\alpha_2$ are dimensional velocity and acceleration weights, respectively. The energy function contains three terms. The first is the least mean square minimisation problem of the polynomial predictor based on the history of particles. The second and the third terms are neighbouring coherent velocity and accelerations. Khojasteh, Heitz, et al. (2021) showed that coherent velocity and acceleration terms could be computed by using Lagrangian Coherent Structures (LCS) to determine. Equal weights of $\alpha_{12}$ mean that the impact of the neighbouring coherent velocity is at the same level as the history of the target particle. However, the particle history should be the most significant signal in the estimation process. These terms can be non-dimensionalised based on turbulent integral scales as,
$$x' = \frac{x}{D}, y' = \frac{y}{D} \qquad (2)$$

$$\dot{x}' = \frac{\dot{x}}{U}, \dot{y}' = \frac{\dot{y}}{U} \qquad (3)$$

$$\ddot{x}' = \ddot{x} \frac{D}{U^2}, \ddot{y}' = \ddot{y} \frac{D}{U^2}, \qquad (4)$$

where D is an integral length scale, and U a velocity reference. Therefore, the energy function will become

$$J = D^2 \left( \frac{1}{n} \sum_{i=1}^{n} (X_i' - y_i')^2 \right) + \alpha_1 U^2 (\dot{X}_n' - \dot{y}_n')^2 + \alpha_2 \frac{U^4}{D^2} (\ddot{X}_n' - \ddot{y}_n')^2, \qquad (5)$$

which can be simplified as,

$$J' = \frac{1}{n} \sum_{i=1}^{n} (X_i' - y_i')^2 + \alpha_1 (\dot{X}_n' - \dot{y}_n')^2 + \alpha_2 (\ddot{X}_n' - \ddot{y}_n')^2. \qquad (6)$$

$$\alpha_1' = \alpha_1(\frac{U}{D})^2, \alpha_2' = \alpha_2(\frac{U}{D})^4 \qquad (7)$$

Two non-dimensional $\alpha_1'$ and $\alpha_2'$ weights determine how much velocity and acceleration signals can constrain the overall energy function minimisation process. In equation 6, weight of the target position history is directly linked with two $\alpha_1'$ and $\alpha_2'$ values. If both terms are set for any range below one, the history is the most significant signal. We need to minimise the energy function for finding the optimal solution for the predictor function. The minimisation process of $J'$ starts by solving $\frac{\partial J'}{\partial a} = 0$. The optimal solution of the minimised energy function leads to a minimum estimation error of the predictor function. We plot the RMS error in a range of $\alpha'$ weights to find the optimal solution.

## 2. Synthetic evaluation

Since a real experiment is associated with uncertainties and inaccuracies, it is crucial to have an appropriate uncertainty estimation of each term in the energy function (see equation (6)). To this end, we employed data reported by the 1st LPT challenge (Sciacchitano et al., 2021) as a reference starting point to estimate inaccuracies that might be introduced into a predictor function. The 1st LPT challenge assessed position estimation accuracy of six time-resolved tracking algorithms, including the coherency based tracking (Khojasteh, Yang, et al., 2021) for particle densities from 0.05 to 0.2 ppp. The non-dimensionalized position $\epsilon_X/D$, velocity $\epsilon_{\dot{X}}/U_{\infty}$, and acceleration $\epsilon_{\ddot{X}}(U_{\infty}^2/D)$ errors were employed to create the synthetic data at each particle density. As a result of the LPT challenge, the averaged RMS position error was 0.005 mm where the integral scale was $D = 10$ mm at the density of 0.12 ppp. The averaged velocity and acceleration error at the same density were found to be 0.01 $U_{\infty}$ and 0.3 $U_{\infty}^2/D$, respectively. This shows that the acceleration estimation has
<figure><img src="image_4.png" /><figcaption>**Figure 2.** Contour plots of the estimation error as a function of uncertainty levels. White arrows show the changing direction of the optimal solution (i.e. minimum estimation error) by increasing each uncertainty term.</figcaption></figure>

at least an order of magnitude higher estimation error than other terms. We introduced these errors into the Lagrangian trajectories of 2D homogeneous isotropic turbulent and 3D cylinder wake flow and then computed the correlation of the final estimation error and $\alpha'$ weights. For such a scenario, we can expect the optimal solution should have less $\alpha'_2$ weight than $\alpha'_1$. Figure 1 shows the RMS prediction error based on two weighing parameters. Since the acceleration error is higher than other terms in both synthetic data, the minimum error happens with relatively lower $\alpha'_2$ values (i.e., $\alpha'_1/\alpha'_2 \gg 1$). Both weighting parameters are found to be significantly smaller than $\ll 1$. This means that the history of the target particle with unit weight is the most valuable signal in such a scenario. Comparing the error behaviour in 2D HIT and 3D cylinder wake flow cases shows that both cases need similar weighting parameters and require $\alpha'_1/\alpha'_2 \gg 1$ with small weighting magnitudes $\alpha' \ll 1$.

The result of the minimum energy function from two synthetic test cases suggests that the optimal solution can be directly linked with uncertainty levels. Therefore, we designed three further parametric test cases to determine the mentioned statement. As shown in Figure 2, we increased the uncertainty level of each parameter solely while other terms are fixed to assess how the optimal solution behaves with respect to three position, velocity, and acceleration uncertainties. The optimal solution tends to linearly move away toward higher weighting magnitude with increased position uncertainty (see Figure 2 a). If we fix the position uncertainty, the optimal solution rotates with the same distance around the coordinate centre, depending on which parameter is increased (see Figure 2 bc). As a result of this assessment, an appropriate estimation of three uncertainty levels would provide enough information to set weighting parameters. On the other hand, we can model the optimal solution by roughly estimating position, velocity, and acceleration uncertainties by fitting a surface over the minimum $\alpha'_{12}$ values. Therefore, we can set the weighting system depending on the uncertainty level of each term.

<figure><img src="image_5.png" /></figure>
3. Experimental evaluation

An experimental study of the cylinder wake flow at Reynolds number equal to 3900 (same value as the synthetic data) was performed. We designed an experimental setup with four cameras. Four CMOS SpeedSense DANTEC cameras with a resolution of 1280 × 800 pixels and the maximum frequency of 3 kHz are empowered. Cameras are equipped with Nikon 105 mm lenses. The first two cameras are positioned in backward light scattering, while the second two cameras receive maximum intensity signal in forward scattering. The calibration error was lower than 0.06 pixel and reduced to 0.04 after the volume self calibration. The volume of interest was 280 mm × 160 mm × 46 mm starting from roughly 4D downstream of the cylinder, knowing that the vortex formation zone ends at 4D. The aperture was set at 11 to achieve 46 mm depth of focus. We used an LED system to illuminate this large volume. The seeding particles were Helium Filled Soap Bubbles (HFSB) Scarano et al. (2015) resulting desired intensity signal with appropriate particle size. However, bubbles are limited by three main factors in the wind tunnel experiments, including generation rate, lifetime, and image glare points. We placed 50 bubble generator nozzles with airfoil-shaped structures inside the wind tunnel chamber. The nozzles were far upstream of the measurement section to ensure a sufficient number of bubbles were created, and the main flow field was not disturbed by the existence of nozzles. The bubble lifetime is very short (less than 2 – 3 minutes) inside the wind tunnel, mainly because they explode by passing through honeycomb layers. To overcome this issue, we injected bubbles inside the chamber for up to 5 minutes when the wind tunnel is off before starting the acquisition. We found that particles larger than three pixels create two glare points on two sides of the bubble. This requires more image treatments before running the 4D-PTV algorithm to avoid false particle reconstruction. However, the intensity of two glare points can diffuse and merge if the particle size is around two pixels. Therefore, we adjusted the camera magnification to reach two particle pixel sizes on average to surpass the glare point issue. One snapshot of the experiment is shown in Figure 3.b. Trajectory results of the current experiment with superimposed vorticity iso-surfaces are shown in Figure 3.b.

To quantify the results of different schemes, we compared predictions with optimised positions obtained from STB Davis. As a result of the experiment, STB managed to successfully build nearly 12000 particles as shown in Figure 3.b. Noisy particle reconstruction of four time steps was used as an input of the prediction functions. We compared three techniques, Polynomial, Wiener filter, and coherent predictors, with final optimised positions. The deviation of position estimated of each technique is shown in Figure 3. The distribution shows that the coherent predictor has more accurate estimations within 1 pixel deviation from the optimised positions. Position estimations of Wiener filter and coherent predictors stay below 2.5 pixels deviation for nearly all particles. On the contrary, the Polynomial predictor has maximum deviation with STB Davis.

<figure><img src="image_6.png" /></figure>
<figure><img src="image_7.png" /><figcaption>**Figure 3.** Experiment study of the cylinder wake flow at Reynolds 3900 with four cameras: a) Experiment normal pdf results of particle position error in x direction of three predictors. Each predictor is compared with final optimised positions of STB Davis; b) Side view of particle trajectories superimposed by vorticity iso-surfaces.</figcaption></figure>

**4. Conclusion**

In the present study, we obtained a generic position predictor as a function of the measurement uncertainties. We performed DNS synthetic analyses of the 3D wake behind a cylinder at Reynolds number of 3900, and 2D homogeneous isotropic turbulent flow at a Reynolds number 3000. It was found that the optimal solution of the predictor function is a case and Reynolds independent model. We also used 4D-PTV experiment data to demonstrate that the generic energy function has better estimation errors than recently available techniques.

**References**

Khojasteh, A. R., Heitz, D., Yang, Y., & Fiabane, L. (2021). Particle position prediction based on Lagrangian coherency for flow over a cylinder in 4D-PTV. In *14th international symposium on particle image velocimetry* (pp. 1–9). Illinois, USA.

Khojasteh, A. R., Yang, Y., Heitz, D., & Laizet, S. (2021, 9). Lagrangian coherent track initialization. *Physics of Fluids, 33*(9), 095113. doi: 10.1063/5.0060644

Scarano, F., Ghaemi, S., Caridi, G. C. A., Bosbach, J., Dierksheide, U., & Sciacchitano, A. (2015, 2). On the use of helium-filled soap bubbles for large-scale tomographic PIV in wind tunnel experiments. *Experiments in Fluids, 56*, 1–12. doi: 10.1007/s00348-015-1909-7

Schanz, D., Gesemann, S., & Schröder, A. (2016, 5). Shake-The-Box: Lagrangian particle tracking at high particle image densities. *Experiments in Fluids, 57*, 1–27. doi: 10.1007/s00348-016-2157-1

<figure><img src="image_8.png" /></figure>
Sciacchitano, A., Leclaire, B., & Schroeder, A. (2021, 8). Main results of the first Lagrangian Particle Tracking Challenge. *14th International Symposium on Particle Image Velocimetry, 1(1)*. doi: 10.18409/ispiv.v1i1.197

Yang, Y., & Heitz, D. (2021, 4). Kernelized Lagrangian particle tracking. *Experiments in Fluids, 62*(12), 1-21. doi: 10.1007/s00348-021-03340-2

<figure><img src="image_9.png" /></figure>