# Vélocimétrie par imagerie de particules 3 dimensions-3 composantes

Gauthier Delerce, Joël Sommeria, Henri Didelle, Samuel Viboud

► To cite this version:

Gauthier Delerce, Joël Sommeria, Henri Didelle, Samuel Viboud. Vélocimétrie par imagerie de particules 3 dimensions-3 composantes. FLUVISU11: Colloque de visualisation et de traitement d'images en mécanique des fluides, Jun 2005, Lyon, France. hal-00261014

HAL Id: hal-00261014

https://hal.science/hal-00261014v1

Submitted on 23 Mar 2021

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.

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

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
Velocimétrie par images de particules
3 Dimensions-3 Composantes

G. Delerce, J. Sommeria, H. Didelle, S. Viboud

GDID, 2 pl du cdt Nal, 38700 La Tronche
LEGI UMR5519, 21 avenue des Martyrs, 38000 Grenoble
gauthier@civx.org

Résumé :

Les techniques de mesures tridimensionnelles des vitesses par imagerie de particules sont en plein essor grâce notamment aux nouvelles caméras (fréquence, résolution et sensibilité) et aux sources lumineuses (puissance, fréquence). Après une présentation de la méthode de balayage de volume (Scanning PIV, (SPIV)) nous décrirons les algorithmes de traitement et de calibration. Nous finirons par une application sur des écoulements à grandes échelles (m3)

Mots-clefs :

PIV 3D

# 1 Introduction

Les méthodes de Particules Imaging Velocimetry (PIV) sont arrivées à maturité pour les mesures dans un plan (2D-PIV), les derniers développements de ces techniques (Scarano 02) sont accessibles à tout chercheur ou ingénieur via des solutions commerciales ou libres. Par constructions, ces techniques sont limitées aux mesures de la projection de la vitesse dans un plan et leur application est limitée aux écoulements dont la composantes de vitesse transverse au plan lumineux est faible. La stéréoPIV (2D3C-PIV) permet quant à elle d'obtenir les 3 composantes (3C) de la vitesse dans un plan et son extension multi-plans (Kähler 99) permet d'obtenir les dérivées. Certains écoulements peuvent être mesurés par de telles techniques, par exemple les écoulements en milieu stratifié ont une composante suivant la stratification très faible ou encore les écoulements en rotation qui ont une dérivée de la vitesse suivant l'axe de rotation très faible (stéréoPIV). Cependant dans un cas général les écoulements sont tridimensionnels et nécessitent des techniques de mesures réellement tridimensionnelles. La technique SPIV utilise un balayage du volume par une tranche laser: une unique caméra enregistre une image pour chaque position de la tranche laser; la juxtaposition des images permettant la reconstitution du volume de mesure, le temps dépendant de la position dans le volume. Le volume de mesure peut être ensemencé par des particules pour un traitement type PIV ou bien par un traceur fluorescent (LIF)

2 Méthode de balayage et résolution suivant z

2.1 Volumes synthétiques

Afin de valider le code de calcul et de déterminer les meilleurs paramètres expérimentaux nous avons simulé l'acquisition de volumes de particules dans un écoulement tridimensionnel connu (tourbillon de Burger, d'axe horizontal) puis nous avons traité ces données avec le code CIV_3D3C présenté ci dessous. La comparaison des résultats avec l'écoulement théorique nous permet de déterminer les meilleurs paramètres suivant les conditions expérimentales. Un volume de dimension
$L_x, L_y, L_z$ est balayé par une tranche lumineuse horizontale (plan $xy$) suivant l'axe vertical (Oz). Le balayage peut être continu ou séquentiel: dans le premier cas la tranche lumineuse est déplacée continûment à la vitesse $V_z$ alors que pour une acquisition séquentielle la tranche est déplacée de $\delta_z$ entre chaque prise de vue. Nous pouvons définir un profil de tranche apparent en intégrant le profil de la tranche sur le déplacement vertical.

Nous considérons dans cet article un balayage continu afin de pouvoir déterminer l'influence de la forme du profil sur la précision des mesures. Pour cela nous avons créé des volumes de particules numériques avec différents paramètres et avons également utilisé un bloc de résine contenant des particules figées que nous avons balayé à différentes vitesses.

Les paramètres que nous faisons varier sont :

- la vitesse de déplacement de la tranche laser et son épaisseur

- la densité de particules et leur diamètre.

Lors du traitement nous testons l'influence de la taille des boites de corrélation, après calcul et comparaison nous analysons les erreurs suivants les trois directions ainsi que le nombre de vecteurs faux.

L'unité de longueur utilisée est le pixel horizontal (pixel) et le diamètre des particules est le diamètre de la tache d'une particule sur une image.

<figure><img src="image_2.png" /><figcaption>Illustration 1: Erreur en fonction de l'épaisseur de la tranche et de sa vitesse de déplacement</figcaption></figure>

L'illustration 1 montre deux résultats, premièrement que le déplacement de la tranche entre deux images doit être inférieur au diamètre des particules, dans le cas représenté les particules ont un diamètre de 3 pixels et l'erreur pour un déplacement de 1 pixel est de 0.1 pixel, pour un déplacement de 3 pixels l'erreur est comprise entre 0.15 et 0.25 pixel et pour un déplacement de 3 pixels l'erreur est de 0.35 pixel. L'autre résultat est que l'augmentation de l'épaisseur de la tranche dans le cas d'un déplacement entre deux images équivalent au diamètre des particules réduit les erreurs.

Ceci montre que pour une acquisition tridimensionnelle l'augmentation de l'épaisseur de la tranche permet de réduire les erreurs sans réduire la résolution verticale, en effet la résolution est fonction du déplacement de la tranche et de la taille des boites de corrélation, or une augmentation de l'épaisseur de la tranche augmente la densité de particule par tranche et permet de réduire la longueur
en z des boites de corrélation.

Ces tests ont montré que la précision et la résolution verticale dépendent principalement de l'épaisseur de la tranche et de sa vitesse de déplacement. Ainsi l'optimum est d'avoir un déplacement de la tranche durant chaque exposition de 30 à 50% de son épaisseur et inférieur au diamètre des particules.

Nous venons de déterminer quelques plages de paramètres pour l'acquisition de volumes de fluides par balayage; après une présentation de l'algorithme de traitement nous verrons l'application de cette méthode à un cas réel d'écoulement en laboratoire.

# 3 Algorithmes

## 3.1 Inter corrélation

A l'instar de l'algorithme CIV (Fincham 97), l'algorithme CIV3D3C n'est pas basé sur l'identification individuelle de particules mais sur la conservation de la texture d'une région de fluide soumise à une petite translation et déformation. L'appariement de zones de fluides entre les deux volumes est déterminé par une fonction d'inter corrélation (1) dépendante du déplacement recherché $\vec{D}$.

Soit deux volumes $^A$ et $^B$ acquis aux temps t et t+dt constitués par la superposition d'images. Pour chaque position $^p$ où nous désirons obtenir une estimation de la vitesse nous définissons un cube $^a_p$ de dimensions $\Box b_x, \Box b_y, \Box b_z$ centré en $\Box x_p, \Box y_p, \Box z_p$ dans le volume A. Nous calculons l'inter corrélation entre $^a_p$ et tous les cubes de même dimensions présents dans un cube centré en $\Box x_p, \Box shift_{x_p}, \Box y_p, \Box shift_{y_p}, \Box z_p, \Box shift_{z_p}$ et de dimension $\Box s_x, \Box s_y, \Box s_z$:

$$C_{\hat{i},j,k} = \frac{\sum_{l=-x_p/2}^{x_p/2} \sum_{m=-y_p/2}^{y_p/2} \sum_{n=-z_p/2}^{z_p/2} \Box_a \Box_d m, n \leftarrow \hat{l}_a \Box_b d, i \leftarrow \hat{i} \Box shift_{x_p}, m \Box j \Box shift_{y_p}, n \Box k \Box shift_{z_p} \leftarrow \hat{l}_b \Box \sum_{l=-x_p/2}^{x_p/2} \sum_{m=-y_p/2}^{y_p/2} \sum_{n=-z_p/2}^{z_p/2} \Box_b \Box_d m, n \leftarrow \hat{l}_b^2 \sum_{l=-x_p/2}^{x_p/2} \sum_{m=-y_p/2}^{y_p/2} \sum_{n=-z_p/2}^{z_p/2} \Box_d l, m, n \leftarrow \hat{l}_b^2}{\hat{l}_a^2}$$

$$\text{avec } -\frac{S_x - B_x}{2} \le i \le \frac{S_x - B_x}{2}; -\frac{S_y - B_y}{2} \le j \le \frac{S_y - B_y}{2}; -\frac{S_z - B_z}{2} \le k \le \frac{S_z - B_z}{2}$$

et

$$\hat{l}_a = \frac{\sum_{l=-x_p/2}^{x_p/2} \sum_{m=-y_p/2}^{y_p/2} \sum_{n=-z_p/2}^{z_p/2} \Box_a l, m, n \Box}{b_x \cdot b_y \cdot b_z}, \quad \hat{l}_b = \frac{\sum_{l=-x_p/2}^{x_p/2} \sum_{m=-y_p/2}^{y_p/2} \sum_{n=-z_p/2}^{z_p/2} \Box_b l, m, n \Box}{b_x \cdot b_y \cdot b_z}$$

$l_a \Box x, y, z$ Étant l'intensité en niveau de gris du voxel à la position x,y,z du volume A

Le déplacement recherché du cube $^a_p$ du volume A pendant l'intervalle de temps dt est définit
comme le maximum de (1); une estimation sous pixel est obtenue par interpolation de cette fonction par un spline 3D lissé.

Nous avons ajouté un paramètre d'estimation local du déplacement dans l'équation (1) $\llap{\text{shift}}_{xp}, \text{shift}_{yp}, \text{shift}_{zp}$ □ afin de mesurer les fluctuations d'un écoulement ayant une composante de vitesse moyenne élevée (composante estimée).

A l'opposé des méthodes par FFT, la dynamique de mesure $\pm \frac{S_x - B_x}{2}; \pm \frac{S_y - B_y}{2}; \pm \frac{S_z - B_z}{2}$ (déplacement maximal autour de la position estimée) est indépendante de la dimension du volume $a_p \llap{\text{b}}_x, \text{b}_y, \text{b}_z$ □.

Le déplacement obtenu est le déplacement moyen de la texture du cube de dimension $\llap{\text{b}}_x, \text{b}_y, \text{b}_z$ □, la résolution spatiale de la mesure est donc de l'ordre de $\frac{\llap{\text{b}}_x, \text{b}_y, \text{b}_z}{2}$.

La corrélation est normalisée par la variance afin de donner un poids aussi important aux textures sombres qu'aux textures lumineuses et ainsi augmenter la sensibilité de la fonction aux variations d'intensité plutôt qu'aux points très lumineux.

## 3.2 Multiplication des fonctions de corrélation.

Suivant les conditions expérimentales la qualité des données varie entre autre par la présence importante de bruits sur les images. Il arrive que ce bruit soit assez important pour faire apparaître dans la fonction d'inter corrélation des maximums supérieurs au maximum correspondant au déplacement recherché. Hart (Hart 00) a remarqué que les maximums dus au bruit apparaissent aléatoirement dans l'espace et a proposé une technique pour améliorer la robustesse des mesures en présence de bruit. En multipliant les fonctions d'inter corrélations de domaines suffisamment voisin (pour que les déplacements correspondants soient proches), nous obtenons une fonction débruitée car les maximums dus aux bruits étant aléatoirement distribués dans le temps et l'espace et donc non corrélés, alors que les maximums correspondants aux déplacements recherchés le sont. Grâce à cette technique nous ne considérerons que les maximums qui sont dans le voisinage du maximum de la fonction obtenue par multiplication. Sans détériorer la précision ni le temps de calcul cette technique permet de supprimer nombre de vecteurs faux sans faire appels à des critères statistiques dépendant de l'écoulement.

## 3.3 Interpolation de la fonction de corrélation

La fonction définie par (1) est une fonction discrète de l'espace, une estimation sous pixels de son maximum est obtenue en interpolant localement cette fonction autour de son maximum.

Westerweel (Westerweel 93) a observé une erreur systématique dans les mesures 2D, ce biais est causé par les fonctions d'interpolations dont le maximum qui est toujours plus proche d'une valeur entière (pixel) qu'il ne devrait l'être (Fincham et Delerce 99). Afin de minimiser ce biais nous privilégions l'utilisation d'un spline lissé à des interpolations par une fonction gaussienne. Ce biais, nommé dans la littérature « peak-locking », est largement diminué en utilisant une méthode itérative; cette méthode sera prochainement implémentée dans cet algorithme.

La réduction de cette erreur systématique est très importante car elle détériore fortement les mesures des petits déplacements ce qui a pour effet de réduire la plage des vitesses mesurables simultanément dans un champ turbulent.
## 3.4 Parallélisation

Le nombres d'opérations pour calculer les corrélations 3D étant important et les systèmes de clusters informatiques se démocratisant nous avons parallélisé ce code afin de minimiser le temps d'attente pour l'obtention d'un champ 3D. Cette approche permet aussi d'utiliser des PC ayant peu de mémoire vive pour le calcul, les volumineux tableaux correspondants aux données et aux résultats étant alloués sur une machine ayant plus de mémoire qui distribue les calculs (sous régions) aux noeuds disponibles.

L'utilisation de caméras rapides et le coût important des nouveaux algorithmes 2D et 3D (Scarano03) font que la distribution des calculs sur un réseau de PC (dédiés ou non) devient incontournable.

Tous les codes développés et utilisés sont distribués sous la licence GNU/GPL (http://www.fsffrance.org/gpl/gpl.fr.html). Chacun peut consulter et améliorer les sources des logiciels, les redistribuer ce qui en assure à terme la pérennité et une large diffusion: http://www.livretdulibre.org

## 3.5 Calibration en milieu stratifié

Afin de pouvoir créer une correspondance précise entre nos volumes de données et l'espace physique lié au laboratoire il est important de suivre une procédure de calibration.

Une première méthode consiste à placer trois tiges orientées à 45° dans le champ de mesure; durant l'acquisition on observe sur les images que la reflection de la tranche laser se déplace et il est aisé de trouver une concordance entre ces déplacements et la position verticale de la tranche. A partir de ces relevés nous définissons un plan par tranche dans l'espace physique et obtenons ainsi une correspondance entre notre volume de mesure et l'espace physique.

Cette première méthode est simple mais est peu adaptée aux fluides stratifiés, en effet la stratification entraîne une courbure de la tranche laser et il est important de la prendre en compte.

### 3.5.2 Deuxième méthode, calibration stéréoscopique

Afin de déterminer précisément la position de chaque point du volume de fluide dans l'espace physique nous utilisons une méthode stéréoscopique.

En plus de la caméra principale située à la verticale du champ, nous plaçons une deuxième caméra faisant un angle de 45° avec la verticale. Pour différents niveaux z, nous plaçons une mire géométrique et capturons une image par caméra et niveau afin de déterminer les caractéristiques de notre système optique. Le modèle optique a été décrit par Tsai (Tsai 86 et Tsai 87).

Nous procédons ensuite à une acquisition continue d'un volume de fluide ensemencé en particule en ayant une synchronisation parfaite des deux caméras.

Pour chaque niveau du volume nous allons maintenant calculer la position suivant z de chaque point :

Pour chaque caméra nous transformons les images (dewrapping) en leur appliquant la transformation Espace_Image vers Espace_Réel. Si la tranche laser est placée dans le même plan que la calibration les images provenant de chaque caméra seront identiques. En utilisant un algorithme de PIV 2D sur cette paire d'image nous déterminons un déplacement apparent qui est du à une différence entre la coordonnée z de la mire de calibration et la tranche laser. Par triangulation nous obtenons la position exacte de la tranche laser en tout point des images.

En répétant cette opération pour chaque plan laser nous créons une relation précise entre les deux espaces en prenant en compte toutes les paramètres physiques : déformation de la tranche par la stratification, passage d'interface par les rayons lumineux, déformations radiales de l'objedif.

# 4- Expériences en laboratoire : Etude des ondes internes sur la plaque tournante de Coriolis (LEGI-Grenoble)

Les ondes de gravité internes dans les océans participent au bilan énergétique en dissipant de l'énergie et mélangent dans les couches profondes de l'océan.
Afin d'étudier ce processus de mélange nous avons modélisé sur la plateforme Coriolis du LEGI à Grenoble la génération d'ondes internes créées par l'arrivée d'une onde barotrope sur un talus continental. L'onde barotrope est créée par l'oscillation sinusoïdale d'un batteur visible à gauche du schéma 1; le talus continental est visible sur la droite de ce même schéma. Le milieu est linéairement stratifié tel que la fréquence de Brunt-Vaisala est $N= 0.72s^{-1}$ ce qui définit l'angle $\omega=N\cos\theta$ de propagation des ondes internes. L'ensemble du système est placé en rotation, la période de rotation est de 90 secondes et la période d'oscillation du batteur est de 25 secondes. L'ensemble du système est placé dans un canal d'une longueur de 14 mètres et de largeur 4 mètres.

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

La zone de mesure est située à proximité talus, la dimension du champ de mesure horizontale est de 180 x 180 cm² et le balayage vertical est de 50cm. La tranche lumineuse est créée par un faisceau laser (argon 8W continu) se réfléchissant sur un miroir oscillant (Cambridge Technology). Ce miroir est placé sur un support de translation motorisé, le faisceau vertical crée est réfléchi sur un miroir à 45° et permet de balayer un volume verticalement par une tranche horizontale.

A partir des mesures tridimensionnelles réalisées, nous pouvons suivre dans le temps l'évolution des ondes créées, en plus de la mesure des angles de propagation des ondes et de leurs harmoniques visible sur des coupes verticales (fig 1) nous pouvons déterminer la direction des déplacements dans les trois dimensions
<figure><img src="image_4.png" /></figure>

La figure 3 est une iso surface de Wx (rotationnel suivant x), les ondes planes sont nettement identifiables, de gauche à droite nous distinguons un groupe d'ondes crées par des jets au niveau du bas du talus, puis l'onde principale se propageant dans la même direction; enfin la première harmonique avec un angle réduit.

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

Les premiers ondes sont crées par des jets espacés de 40 cm, à l'opposé des ondes planes les déplacements ainsi crées ont un rotationnel Wy non nul que l'on peut visualiser sur la figure 4 qui est une iso surface Wy vu par le dessus.
<figure><img src="image_6.png" /></figure>

La tranche laser a une épaisseur de 11 mm, le volume de 500 mm est parcouru en 1.5 secondes et constitué de 93 tranches.

Le rapport (déplacement verticale de la tranche)/(épaisseur de la tranche) est proche de 50%.

Une caméra Dalsa DS41-1M60 (capteur 1024x1024 pixels, fill-factor 100%, résolution 12 bits, fréquence 60Hz) enregistre en continu les séries d'images à la fréquence de 2.5 secondes par volume (1.5 s d'acquisition et 1 s pour le repositionnement de la tranche lumineuse et l'enregistrement des images sur le disque dur). Les calculs ont été réalisés avec des boites de corrélation de 19x19x7 voxels, les déplacements verticaux mesurés sont de 7 pixels. Dans la figure ci-après montrant une coupe verticale du volume mesuré on peut comparer la résolution horizontale et verticale. Malgré une méthode d'acquisition des données privilégiant la résolution horizontale nous obtenons des résultats semblables dans les deux dimensions, on peut par exemple observer la continuité de l'écoulement mesuré.

Ces mesures ont prouvé la faisabilité expérimentale de cette technique, l'écoulement choisi était lent car nous étions limités techniquement par le mécanisme de translation de la tranche lumineuse, par la fréquence maximale d'acquisition des images, ainsi que par la puissance lumineuse du laser. La disponibilité d'une caméra rapide ainsi qu'une technique de balayage différente nous permettait maintenant de nous intéresser à un écoulement turbulent.

**Conclusion:**

Nous avons présenté une méthode d'acquisition tridimensionnelle et son application sur des écoulements complexes de grandes dimensions. L'analyse des simulations et des résultats sur des
expériences réelles ont montré la maturité de la méthode et sa facilité d'utilisation. Le faible coût du système d'acquisition, associé à la souplesse d'utilisation du matériel de balayage, assurera à cette méthode une large diffusion.

Références:

1- Scarano F Iterative image deformation methods in PIV, Meas. Sci. Technol. (2002) 13 R1-19

2- Kähler CJ, Kompenhaus J, Multi-plane stereo PIV – Technical realization and fluidmechanical significance, 3rd International Workshop on PIV, University of California-Santa Barbara, USA, Paper PIV96, Sept. 16-18 1999

3- Trigui et al Characterization of intake generated fluid flow fields in IC engines using 3-D particle tracking velocimetry (3-D PTV) SAE-paper 940279

4- Hinsh K Three-dimensionnal particle velocimetry Meas. Sci. Technol. 6 742-53 (1995)

5- Brücker CH. 3D scanning PIV applied to an air flow in a motored engine using digital high-speed video, Meas. Sci. Technol 8(1997) 1480-1492

6- Lester K.Su and Werner J.A. Dahm, Scalar imaging velocimetry measurements of the velocity gradient tensor field in turbulent flows. Phys. Fluids 8, July 96

Fincham, A., and G. Spedding, 1997, Low cost, high resolution DPIV for measurement of turbulent fluid flow. Exp. Fluids, 23, 449-462 defined the CIV techniques and several acquisition optimizations.

7- Hart DP DPIV error correction Exp. Fluids 29 13-22 (2000)

8- Rignot E.J.M. & Spedding G.R. Performance Analysis of automated image processing and grid interpolation techniques for fluid flows, USCAE 143, March 88

9- Scarano F. A super-resolution particle image velocimetry interrogation approach by means of velocity second derivatives correlation, Meas. Sci. Technol(2004) 15 475-486

10- Mowbray, D.E. and B.S.H.Rarity (1967). A theoretical and experimental investigation of the phase configuration of internal waves of small amplitude in a density stratified liquid. *J.Fluid Mech.*, **28**, 1-16.

11- References instabilités de Görtler : Cousteix, J., 1989, Aérodynamique II: turbulence et couche limite, CEPADUES eds., Toulouse

"An Efficient and Accurate Camera Calibration Technique for 3D Machine Vision", Roger Y. Tsai, Proceedings of IEEE Conference on Computer Vision and Pattern Recognition, Miami Beach, FL, 1986, pages 364-374,

and

"A versatile Camera Calibration Technique for High-Accuracy 3D Machine Vision Metrology Using Off-the-Shelf TV Cameras and Lenses", Roger Y. Tsai, IEEE Journal of Robotics and Automation, Vol. RA-3, No. 4, August 1987, pages 323-344.