FLUIDOS

Desarrollo de un sistema de velocimetría por seguimiento de partículas para la caracterización lagrangiana de flujos

Development of a particle tracking velocimetry system for the lagrangian characterization of fluid flows

 

Lucía Cappelletti^1,2, Nicolás Del Grosso^1,2, and Pablo Cobelli^*2,3

¹ Estudiantes de la Licenciatura en Ciencias Físicas
² Departamento de Física J.J. Giambiagi, Facultad de Ciencias Exactas y Naturales, UBA, Pabellón 1, Ciudad Universitaria, C.A.B.A., Argentina
³ Instituto de Física de Buenos Aires (IFIBA), CONICET
* cobelli@df.uba.ar.

Recibido: 16/06/18;
Aceptado: 07/08/18

---

En este trabajo se presenta el desarrollo y la caracterización de un sistema de velocimetría por seguimiento de partículas para la medición lagrangiana de flujos. El sistema está basado en el uso conjunto de dos cámaras rápidas que observan sincronizadamente una región común del flujo a estudiar. Se describen en detalle el procedimiento de calibración conjunta del sistema de cámaras, la evaluación de los errores de medición y el procesado de las imágenes individuales para el cálculo de trayectorias de partículas en el espacio 3D. En último lugar, se discuten los resultados de la aplicación de este sistema al estudio de la dinámica de clústeres de partículas flotantes sobre la superficie libre de un campo de ondas turbulento.

Palabras Clave: Velocimetría; Seguimiento de partículas; Descripción lagrangiana de flujos.

In this paper we present the development and characterization of a velocimetry system based on individual particle tracking, specially designed for the lagrangian description of fluid flows. This system is based on the joint use of two high-speed cameras which synchronically observe a common region on the flow under study. We describe in detail the joint calibration procedure of the system, as well as the processing of individual images for calculating particle trajectories in three-dimensional space. Lastly, results are shown concerning the application of this measuring system to the dynamics of buoyant particle clustering on the free surface of a turbulent wave field.

Keywords: Velocimetry; Particle tracking; Lagrangian description of flows.

---

 

I. Introducción

La dinámica de partículas en flujos complejos es una temática de fundamental importancia tanto desde la física básica como desde el punto de vista de su implicancia para las ciencias atmosféricas y/u oceanográficas. Esto es particularmente cierto cuando se consideran fenómenos tales como el transporte de cenizas volcánicas, el transporte en meddies¹ y la dispersión de islas de basura. Estos tres ejemplos arquetípicos tienen en común la presencia de partículas cuyo movimiento se da en un flujo complejo (alta variabilidad espacial y temporal). Por ello, resulta esencial disponer de mediciones de la dinámica de estas partículas, tales como su velocidad y aceleración, lo que permitiría profundizar la comprensión de procesos de física fundamental, como por ejemplo, las signaturas en turbulencia ^{[1, 2]}. Por esto, se propuso desarrollar la técnica de PTV (Particle Tracking Velocimetry) 3D multipartícula y ensayarla en una situación que mimique el movimiento de partículas en la superficie oceánica en condiciones controladas de laboratorio. Dicha técnica consiste en registrar el movimiento de partículas en un flujo dado mediante varias cámaras y a partir de esto triangular su posición en el espacio tridimensional. Esta última etapa permite recuperar las trayectorias de las cuales, posteriormente, se derivan la velocidad y aceleración lagrangianas, para así estudiar el flujo en cuestión ^{[3, 4]}.
Para el estudio y caracterización de flujos complejos, existen diversas técnicas experimentales que permiten determinar el campo de velocidades euleriano instantáneo del fluido o la trayectoria de las partículas que lo conforman ^[3]. Una de estas técnicas, conocida como PIV (Particle Image Velocimetry), consiste en inseminar con partículas trazadoras al flujo que se busca estudiar.
Estas últimas, cuyas dimensiones y masa deben ser tales que su movimiento acompañe al de las partículas elementales del fluido, son iluminadas y su movimiento es registrado por una cámara ultra rápida. Dependiendo del tipo de iluminación (se suelen emplearse láseres de alta potencia) y la cantidad de cámaras que registren el movimiento de las partículas (por ejemplo, microesferas huecas de vidrio cubiertas de plata), la técnica PIV permite recuperar el campo euleriano instantáneo de velocidades ya sea en un plano (PIV 2D) o en volumen (PIV 3D). Esta técnica, sin embargo, resulta inadecuada para abordar y comprender problemáticas más complejas como el mezclado y el transporte de materia, cuyo entendimiento requiere de la capacidad experimental de seguir a las partículas en su movimiento ^[4]. Por ello, en muchos casos de interés, tanto en flujos de laboratorio como industriales y/o naturales es imprescindible disponer de una descripción lagrangiana de flujo. Esta última trata de una caracterización del transporte de una sustancia en un flujo, en la cual las partículas de fluido guardan identidad y son seguidas en su movimiento a medida que el tiempo transcurre. Estas razones han llevado, en los últimos años, al surgimiento de nuevas técnicas para la caracterización lagrangiana de flujos. Una de ellas es la PTV 3D multipartícula, que consiste en seguir a determinadas partículas del flujo, para hallar entonces su trayectoria y a partir de ello su velocidad y aceleración. Esta última es una herramienta que en la actualidad es empleada en laboratorios de investigación en dinámica de fluidos del mundo entero, sobre todo para el estudio de escenarios físicos donde la turbulencia tiene un rol preponderante. En particular, en este trabajo se aplicó la técnica PTV 3D-MP al estudio de la dinámica de partículas en la superficie libre de un campo de ondas turbulento.

II. Diseño del sistema de velocimetría por seguimiento de partículas

Montaje experimental
Con el objetivo de desarrollar la técnica PTV 3D MP se propuso ensayarla en el estudio de un sistema que mimicara el estado oceánico medio. Para ello se empleó un tanque de ondas construido en PMMA, de 198x77x13,6 cm³. Como líquido de trabajo se utilizó agua destilada, con una profundidad de trabajo (en reposo) de 5 cm. En cuanto a las partículas a seguir, se diseñaron esferas huecas de PVC realizadas por impresión 3D, como la que se muestra en la Figura 1. Estas presentaban 15 mm de diámetro, 1 mm de espesor y un peso de (0, 84 ± 0, 04) g. Las partículas fueron diseñadas procurando su buena flotabilidad, baja interacción mutua y con el fluido, además de realizarlas en color negro para facilitar su detección. Asimismo, se agregó TiO2 al agua, al fin de teñirla de blanco (sin modificar significativamente sus propiedades reológicas) a fin de aumentar considerablemente el contraste partículas–fluido. El campo de ondas turbulento se generó mediante una agitación controlada por dos wavemakers independientes. Para ello, dos motores lineales tubulares servocontrolados independientemente suministran un movimiento a dos paletas de acrílico parcialmente sumergidas que ofician como generadoras de oleaje. Dichos motores, cuya fuerza pico es de 47 N y cuya precisión en posicionamiento es menor que 0,01 mm, aseguran la generacion del oleaje con las características deseadas y la repetibilidad de las experiencias.

Figura 1: Panel izquierdo: Floaters empleados en el presente estudio. Se diseñaron esferas huecas de PVC, realizadas por impresión 3D, de 15 mm de diámetro, 1 mm de espesor y un peso de (0, 84 ± 0, 04) g. Panel derecho: Dispositivo experimental montado para llevar a cabo las mediciones. Se señala en colores el sistema de coordenadas empleado.

En cuanto al movimiento impuesto a los wavemakers, conviene mencionar que los mismos fueron excitados con una señal de ruido blanco filtrado en frecuencia en el rango 0 - 4 Hz, y de amplitud (i.e., potencia cuadratica media) variable en forma controlada. Se empleó este tipo de forzado de manera tal de reducir los efectos de capilaridad del sistema y generar turbulencia débil de ondas de gravedad, característica del flujo que se buscó reproducir. La señaal empleada para el forzado, a través de la relación de dispersión de ondas de gravedad, determinó una escala temporal de 0,25 s y una espacial de 10 cm. Para registrar la posición de las esferas se filmaron videos con dos cámaras ultra rápidas marca Photron modelos 1024 PCI y SA3. éstas se ubicaron con sus ejes ópticos cruzados a aproximadamente un ángulo de 90º entre sí. La extensión de las filmaciones fue de 2726 cuadros, siendo la memoria interna de las cámaras (de 4 Gb cada una) el principal limitante en dicho aspecto. El montaje experimental dispuesto se muestra en la Figura 2.

Figura 2: Panel izquierdo: Montaje completo con sistema PTV y tanque de ondas, construído en el Laboratorio de Turbulencia Geofísica (DF, FCEN, UBA) para el presente estudio. Panel derecho: Vista en detalle sobre la cuba de ondas, el sistema generador de oleaje y las cámaras que integran el sistema de PTV.

Considerando el objetivo de recuperar las trayectorias 3D de los flotadores, resultó necesario que la captura de imágenes por ambas cámaras se realizara en forma sincrónica. Dadas las escalas temporales involucradas se procuró que las incertezas en la sincronización fuesen, en todos los casos, menores a 0,01 s. Esto se logró mediante el empleo conjunto de una señal de trigger externo provista por un generador de funciones, que permitió utilizar una frecuencia de adquisición de 120 Hz. Por otro lado, dichas escalas temporales forzaron la utilización de iluminación halógena de alta potencia y bajo parpadeo.
En último lugar, la reconstrucción de las trayectorias 3D a partir del arreglo de visión estéreo antes mencionado requiere de un modelo matemático que describa adecuadamente el registro que del mundo hacen las cámaras, permitiendo pasar de observaciones en dos dimensiones (los CCDs de cada cámara) a posiciones en el espacio tridimensional. La descripción de dicho modelo, su caracterización y utilización es objeto de la sección siguiente.

Visión estéreo
Las cámaras empleadas fueron modeladas por una cámara pinhole con distorsión radial. El modelo de la cámara pinhole consiste, básicamente, en aproximar el funcionamiento de una cámara compleja por aquel asociado a una cámara oscura; es decir, una caja en forma de paralelepípedo recto, sin lentes y con una pequeña abertura en una de sus caras. En este modelo, como se muestra en la Figura 3, la luz ingresa por la abertura y proyecta una imagen invertida en el lado opuesto de la caja. En base a esto un punto 3D de una escena, con coordenadas dadas por [X, Y, Z], se relaciona con su proyección 2D en una imagen, de coordenadas [x, y, 1] (en pixels) mediante la relación:

donde w es una constante de escala y P es una matriz de 4x3. La matriz P contiene los denominados parámetros extrínsecos e intrínsecosde la cámara, los cuales son esenciales para describir el funcionamiento de la cámara pinhole. Los parámetros extrínsecos refieren a la ubicación de la cámara en la escena tridimensonal y se encuentran dados por una rotación R y una traslación T respecto del sistema de coordenadas del laboratorio. Por otro lado, los parámetros intrínsecos se reúnen usualmente en la denominada matriz de la cámara, K, y describen la posición del centro óptico (U0, V0) y la distancia focal de la cámara. Tanto K como R son matrices de 3x3, mientras que T es un vector de 1x3. Las mismas se relacionan por medio de la expresión

Figura 3: Diagrama de rayos para mostrar el funcionamiento de una cámara pinhole. Se caracteriza por ser un modelo de cámara simple sin lentes y con una pequeña abertura. En el marco de este modelo simple, y para evitar invertir las imagenes, se considera al plano imagen como ubicado a una distancia f por delante del centro proyectivo, según se obseva en el panel izquierdo. El panel derecho muestra una vista lateral del mismo esquema.

De esta manera, los puntos del mundo se transforman a coordenadas de la cámara a través de los parámetros extrínsecos, mientras que los parámetros intrínsecos permiten asignar las coordenadas de la cámara en el plano de la imagen [5, caps. 2 y 6]. En contraposición con el modelo sencillo que se viene describiendo, una cámara real posee lentes, las cuales indefectiblemente introducen distorsiones u aberraciones ópticas sobre la imagen registrada, a menudo en forma no lineal. Por este motivo, el modelo de cámara empleado en este trabajo contempla la posibilidad de la existencia de dichas aberraciones ópticas para lograr así una reproducción más precisa del funcionamiento de las cámaras empleadas en las experiencias. De esta forma, resulta posible compensar las distorsiones ópticas introducidas involuntariamente por el uso de dichas cámaras, mejorando considerablemente la precisión en la triangulación de las partículas.
La distorsión más importante a considerar es la radial, que afecta a los puntos de la imagen de coordenadas (x, y) de acuerdo a

donde (xd, yd) son los puntos distorsionados, k1 y k2 son los coeficientes de la distorsión radial y r2 = x2 + y2. Los parámetros k1 y k2 se obtienen a partir del análisis de imágenes de calibración, lo que permite luego desdistorsionar las imágenes capturadas por las cámaras que integran el arreglo. Una vez completada dicha etapa, la obtención de la trayectoria tridimensional conociendo los puntos de la imagen requiere invertir la matriz P. Esta inversa no es única, dado que existen infinitos puntos 3D que se proyectan al mismo punto 2D. Para solucionar esto, resulta necesario disponer de varias cámaras para luego poder recuperar la posición en el mundo (3D) mediante triangulación ^[6].
El proceso de triangulación consiste en, conociendo los parámetros intrínsecos y extrínsecos del sistema estéreo y las proyecciones (Pi, Pd) de un punto P del mundo en los planos de imagen de ambas cámaras, determinar la ubicación 3D del punto del mundo del que provienen. El punto P en coordenadas de la cámara izquierda será P = a Pi (con a un número real) mientras que en coordenadas de la cámara derecha será P = b RtPd + T (donde R y T son la rotación y traslación relativa entre las cámaras y b un número real). Sin embargo, en la práctica este sistema no tendrá solución debido a las imperfecciones que introduce la utilización de cámaras reales. Por esta razón es necesario considerar una traslación adicional, perpendicular a ambos vectores, de la forma c (Pi RT Pd), siendo c un número real. Con esta modificación el sistema final a resolver será

de donde se obtienen los factores a, b y c, a partir de los cuales se obtiene el punto triangulado dado por la expresión simple P = a Pi + c2 (Pi RT Pd).
La concresión de esta última etapa plantea el desafío experimental de establecer unívocamente la correspondencia entre un punto observado por una cámara y su correspondiente en la región de observación de la otra cámara. Este problema de correspondencia puede ser resuelto mediante el uso de la geometría epipolar ^{[5, cap. 7]}. En este sentido, la geometría epipolar permite reducir el problema de la correspondencia, que supone la búsqueda 2D de un punto sobre la imagen de una cámara a un problema unidimensional de búsqueda de dicho punto sólo sobre una recta (ver ^{[7] y [8, caps. 9 y 11]}). Para ello se le asigna a un punto de la cámara i una recta (línea epipolar) sobre la cámara d, dicha recta representa todos los posibles candidatos que pueden corresponder al primero, tal como se muestra en la Figura 4.

Figura 4: Representación de la geometría epipolar de dos imágenes. II es el plano imagen donde se hallan las proyecciones de los objetos 3D y O se refiere al centro óptico de la cámara (i ó d).

Matemáticamente la línea epipolar proyectada en la cámara d por el punto de imagen pi se encuentra dada por la expresión

donde F se conoce como matriz fundamental y se encuentra definida por

donde Tx, Ty y Tz son las componentes cartesianas del vector de traslación T; mientras que Kd y Ki son las matrices de parámetros intrínsecos de las cámaras d e i respectivamente. Estas herramientas teóricas desarrolladas en el presente apartado se emplean en la sección siguiente, en donde se realiza la determinación experimental de los parámetros intrínsecos y extrínsecos del sistema de PTV 3D-MP, junto con sus incertezas.

Calibración
Según se expuso en el apartado anterior, a fin de obtener la posición 3D de un objeto a partir de imágenes de 2D de múltiples cámaras resulta necesario conocer los parámetros extrínsecos e intrínsecos del arreglo estéreo. El proceso por el cual se determinan estos se conoce como calibración, y constituye el objeto de la presente sección. Para llevar a cabo este proceso de calibración existe una gran diversidad de algoritmos ^{[9, 10, 11, 12]}. En el presente estudio experimental se optó por emplear el método de Zhang^[13]. El principio de trabajo de este método puede descomponerse en las siguientes cuatro etapas. La primera consiste en tomar varias imágenes de un blanco 2D de características conocidas (por ejemplo, una greilla regular de puntos donde se conocen las distancias entre vecinos) empleando la cámara que se desea calibrar. A continuación se realiza el reconocimiento de dichos puntos en la imagen 2D, se les asigna a cada uno de ellos los puntos correspondientes del mundo (en 3D) y mediante la resolución de un sistema de ecuaciones se obtienen los valores de los parámetros buscados.
Esta técnica permite calibrar de manera fácil una cámara, además de ser más flexible con respecto a otros métodos. Su flexibilidad radica en que sólo requiere que la cámara a calibrar observe un patrón plano mostrado en distintas orientaciones (en teoríia, al menos dos). Por otro lado, tanto la cámara como el blanco 2D pueden ser movidos libremente, sin la necesidad de que este movimiento sea conocido por el usuario. Una ventaja adicional de este método, es que el blanco de calibración no requiere un diseño de características particulares, emple´andose com´unmente un arreglo regular de cuadros alternativamente blancos y negros en ambas direcciones del plano (en adelante llamado damero).
Matem´aticamente, el algoritmo de Zhang puede expresarse en la siguiente forma: dados los puntos de la imagen, ˜m = [u, v, 1]T , y los puntos del mundo ˜M , se busca establecer la homografía H que los vincula mediante la ecuación

siendo s un factor de escala y

donde r1 y r2 denotan la primera y segunda columna de una matriz de rotación R, respectivamente. Dicha homografía H se obtiene mediante un proceso de ajuste no lineal de cuadrados mínimos. Para ello, se propone un guess inicial calculado a partir de la siguiente expresión

De esta igualdad se nota que se puede aplicar una descomposición en valores singulares (SVD) a la matriz LTL y obtener la homografía buscada, dada por el vector asociado valor singular de menor magnitud. No obstante, dado que L está, por lo general, mal condicionada, resulta imprescindible normalizar los puntos previamente (restando la media y dividiendo por la desviación estándar). Seguidamente, se aplica el proceso de minimización de cuadrados no lineales al funcional

donde mˆi = 1 h3 ˜Mi [h1 ˜Mi; h2 ˜Mi]. En el caso de la implementación llevada a cabo, se empleó el algoritmo ‘Trust Region’ ^[14] como minimizador, dado que provee una rápida convergencia y puede manejar tanto problemas grandes y ralos como pequeños y densos. Asimismo, es posible establecer con mayor precisión los parámetros intrínsecos, considerando homografías de fotos del blanco en distintas orientaciones. A partir de los parámetros intrínsecos, se pueden deducir los correspondientes extrínsecos mediante la expresión

donde b = [B11,B12,B22,B12,B23,B33], siendo B = K-TK-1 y vij = [hi1hj1, hi1hj2 + hi2hj1, hi2hj2, hi3hj1+hi1hj3, hi3hj2+hi2hj3, hi3hj3] y hij la i-ésima columna de la j-ésima homografía. Si dicha relación se escribe como V b = 0 se puede resolver haciendo SVD de V y hallar el vector correspondiente al menor valor singular. De esta manera se obtiene un guess inicial para los parámetros intrínsecos y extrínsecos, los cuales sirven de entrada a un proceso de ajuste no lineal de cuadrados mínimos, dado por la minimización del funcional definido como

donde ˆm(A,Ri, Ti,Mj) representa la proyección del punto Mj en la imagen i, de acuerdo a la ecuación 7. De esta forma se obtienen los parámetros extrínsecos e intrínsecos de cada una de las cámaras del arreglo.

Implementación de los algoritmos de PTV
Para lograr tener un control meticuloso sobre cada paso del proceso de calibración, se decidió programar códigos propios. El procedimiento implementado para calibrar cada una de las cámaras empleadas consistió en las siguientes etapas:

Diseñar un blanco 2D y fotografiarlo:
Se diseñó un damero de 10x7 casilleros de 25,4 mm de lado para emplear como blanco 2D. Se optó por este tipo de patrón ya que asegura una fácil identificación de los cruces de los casilleros que lo conforman, aún frente a grandes rotaciones. Se tomaron aproximadamente 100 fotografías del damero, variando su posición y orientación, procurando cubrir todo el campo de visión de las cámaras y que las imágenes se tomaran de tal forma que el patrón abarcase la mayor superficie de la imagen (Figura 5).

Figura 5: Pequeña muestra de las imágenes utilizadas para la calibración. Se observa la amplia variación angular y traslacional empleada para la disposición del damero a la hora de tomar las fotografías.

Identificar los cruces de los casilleros:
Se implementó un código basado en el método de Harris ^[15] para realizar la detección de los cruces de los casilleros, por ser robusto y de fácil implementación respecto a otros algoritmos de detección de esquinas^{[16, 17]}.

Numerar los cruces identificados:
Como resulta natural, los cruces de los casilleros identificados en la etapa anterior se obtienen sin ningún orden predeterminado. A fin de poder emplear estos puntos como sistema de referencia del mundo, resulta necesario ordenarlos y asignarle a cada uno su vector del mundo 3D correspondiente (ver Figura 6).

Figura 6: Detección de los cruces de los casilleros del damero (círculos rojos) junto con algunos de los puntos del mundo asignados (todos se encuentran sobre un plano siendo su coordenada Z igual a 0).

Obtener los parámetros intrínsecos y extrínsecos:
El código para la calibración de las cámaras fue estructurado considerando el método de Zhang (op. cit.). éste último se sirve de las coordenadas de los puntos situados en el blanco 2D y devuelve los parámetros intrínsecos y extrínsecos de una de las imágenes para cada cámara.

Determinar la imagen óptima como sistema de referencia del mundo:
Para esta etapa se generó un código que selecciona el par de imágenes de ambas cámaras que presenta el menor error de reproyección (en valor absoluto). A partir de esta imagen se determinaron experimentalmente los parámetros extrínsecos del arreglo.

Una vez completado este procedimiento, se obtuvieron los parámetros intrínsecos de las cámaras utilizadas. Los resultados obtenidos fueron los siguientes: U0 = (530 ± 150) px, V0 = (627 ± 150) px y focal (55, 35 ± 0, 5) mm para la cámara i y U0 = (529 ± 150) px, V0 = (477±150) px y focal (29, 01±0, 5) mm para la cámara d. Resulta conveniente destacar que estos valores son consistentes con los esperados para la posición del centro ´óptico de las lentes empleadas y para la distancia focal de las mismas. Tanto las incertezas que surgen del proceso de calibración, así como también los del método se presentan en la Sección III.

III. Caracterización del sistema PTV 3D-MP

En esta sección se describe el proceso de caracterización experimental del sistema de velocimetría por seguimiento de partículas. En primera instancia se llevó a cabo una evaluación cualitativa de los parámetros adquiridos mediante la calibración. Para ello se renderizó de manera virtual la posición del damero y las cámaras junto con sus CCD en el espacio tridimensional, sobre los cuales se proyectaron las imágenes adquiridas y se trazaron los rayos de proyección (ver Figura 7). De este modo, se comprobó que la intersección de los rayos con el CCD coincidían con los cruces de las esquinas, presentando discrepancias de magnitud inferior al pixel. Además, se realizó un estudio cuantitativo de las incertezas experimentales. Para éste, se tomaron los puntos del mundo y se reproyectaron para cada imagen empleando sus parámetros correspondientes. De esta manera se obtuvieron, de la diferencia entre estos últimos y los medidos, los errores de reproyección (Er) para cada cámara.

Figura 7: Renderización virtual del damero junto con las dos cámaras y sus CCDs. Se proyectó una fotografía sobre el CCD donde se observan los rayos de proyección (líneas azules) pasando correctamente por los cruces de las esquinas (círculos rojos).

Una vez adquiridos los Er, se realizaron PDFs de los mismos a partir de 11340 puntos para ambas cámaras, donde las componentes X e Y fueron consideradas como dos puntos distintos entre sí. Las PDFs se ajustaron de manera adecuada con una distribución normal, lo cual se puede notar en el panel izquirdo de la Figura 8. De esto se deduce que los errores son independientes de las mediciones y de la dirección.

Figura 8: Panel izquierdo: PDFs de los errores de reproyección de ambas cámaras, cada una ajustada por una distribución normal. Panel derecho: PDF de los errores acumulados del método de detección.

Asimismo, se obtuvo para los errores de la cámara i una desviación estándar de 0, 33 px y una media de 2, 30 x 10^-5 px; mientras que para la cámara d, la desviación estándar fue 0, 24 px y la media 6, 66 x 10^-6 px. Estos resultados nos permiten asegurar que el método de calibración presenta un nivel de precisión sub-pixel. Por otro lado, para analizar los errores acumulados de la técnica, se trianguló la posición 3D de 34 imágenes de un damero conocido. Evaluando la diferencia entre las distancias del triangulado y del real, se obtuvieron los errores totales de una posición en 3D (E3D). Se realizó una PDF de estos errores (Figura 8, panel derecho), la cual no ajusta a una distribución normal sino que decae más lentamente. El comportamiento lineal en semi-log parece indicar, asimismo, un decaimiento exponencial. Para esta PDF se obtuvo una desviación estándar de 0, 33 mm y una media de 0, 008 mm, valores dos órdenes debajo de la escala más pequeña a estudiar.
Cada uno de estos tableros triangulados se ajustaron por un plano y se calculó la distancia de cada punto al mismo, obteniendo así, los errores Ep. La desviación estándar de Ep resultó ser de 0, 56 mm, en acuerdo con el error acumulado del método e indicando una buena coplanaridad de los puntos. A partir de lo expuesto en esta sección, se puede afirmar que las evaluaciones de las incertezas experimentales mostraron que el método de detección es submilimétrico. De esto, y teniendo en cuenta que la longitud característica del problema a estudiar es del orden de 10 cm, se puede concluir que el sistema posee una precisión satisfactoria para los objetivos de este trabajo.

IV. Obtención de las trayectorias individuales de partículas en 3D

Haciendo uso del montaje descripto en la Sección II se obtuvieron videos del movimiento de los flotadores, a partir de los cuales se buscó recuperar las trayectorias 3D de las partículas. Para ello, empleando lo expuesto en la Sección II, se desarrolló la estrategia de procesamiento de datos que se detalla en los siguientes párrafos.
Compensación de la distorsión: Fue necesario compensar la distorsión radial de las imágenes asociadas a cada cámara para tratar el problema de la correspondencia (detallado previamente).
Detección de círculos: Se detectaron círculos oscuros en cada fotograma, hallando en forma subsiguiente su centro y radio. Seguidamente, se recortó la imagen de cada círculo y se realizó una nueva detección de círculos en la imagen reducida con el fin de mejorar los resultados de la detección original y aumentar la confiabilidad del algoritmo.
Obtención de trayectorias 2D en cada CCD: Una vez detectadas las partículas, se obtuvieron las trayectorias 2D sobre cada imagen de ambas cámaras. Para ello, dado un centro de un círculo en un fotograma, se consideró que se correspondía con otro centro en el fotograma siguiente siempre que la distancia entre ambos fuese inferior a una tolerancia ejustable a criterio del usuario (en nuestro caso, dicho umbral se estableció en un valor de 10 px, que representa un 1% de la dimensión lineal de las imágenes registradas). De esta forma, se vincularon entre sí las posiciones de las partículas en todas las imágenes siguientes, dando como resultado las trayectorias 2D.
Correspondencia entre imágenes: Dadas dos colecciones de trayectorias 2D en cada cámara, la siguiente etapa consistió en establecer cuáles se corresponden a una misma partícula. Para ello se tomaron curvas de a dos de cada colección, que se solaparan temporalmente. Es decir, que el tiempo inicial de la primera fuese menor que el final de la segunda. Entonces, para cada par de puntos de las trayectorias se trazaron las líneas epipolares y se pudo establecer o no la correspondencia, dependiendo de si los puntos estaban a una distancia de la línea epipolar generada por el otro menor al valor umbral antes mencionado. De esta forma, para cada trayectoria se obtuvo un porcentaje de correspondencia. Se tomó como criterio que dos trayectorias 2D correspondían a una misma partícula si este porcentaje superaba el 90%.
Triangulación: Dadas dos trayectorias 2D en cada cámara de la misma partícula, es posible recuperar la trayectoria 3D triangulando instante a instante la posición 3D, para lo que se requiere conocer los parámetros intrínsecos y extrínsecos de ambas cámaras determinados durante el proceso de calibración.
Determinación de la superficie libre: En el sistema a estudiar existe un eje de coordenadas cartesianas cuya orientación resulta natural y que corresponde al formado por el plano de la superficie libre del fluido y la dirección perpendicular a esta. Por ello, fue necesario determinar la posición de la superficie libre. Para llevarlo a cabo se dispusieron alrededor de diez flotadores en el campo de visión de las cámaras, con el fluido en reposo, y se triangularon sus posiciones. Seguidamente, estas posiciones se ajustaron por un plano y se determinó su centro. A partir de éstos, se determinó la rotación entre el sistema de referencia usado y el plano hallado.
Rectificación: Una vez obtenidas las trayectorias 3D con el sistema de referencia obtenido de la calibración, se las rectificó. Este procedimiento consiste en transformar sus componentes a un nuevo sistema de referencia. Dicho sistema se eligió tomando como origen el centro encontrado en el inciso anterior, como ejes a X e Y en el plano de la superficie libre, y a Z perpendicular este.
El procesamiento de datos detallado se expone de manera ilustrativa en la Figura 9.

Figura 9: Cuatro etapas principales de la determinación experimental de trayectorias individuales de partículas, a partir de imágenes, mediante PTV.

V. Obtención de magnitudes derivadas: velocidad y aceleración 3D

Una vez obtenidas las trayectorias 3D de las partículas, resulta deseable calcular magnitudes derivadas de éstas, como la velocidad y la aceleración. Sin embargo, al realizar la operación de derivación de forma numérica empleando diferencias finitas, el error propio de la trayectoria se ve amplificado con cada derivación. Esto introduce ruido de alta frecuencia a las mediciones y hace necesario un post-procesamiento de los datos que permita disminuir el efecto. Dicho post-procesamiento consistió en filtrar las trayectorias, haciendo uso de un spline cúbico con una ventana adaptada a las condiciones físicas de la dinámica de la partícula. Se desarrollaron códigos para implementar el filtrado, basados en el trabajo de Luthi^[18]. Estos consistieron en, para cada instante de tiempo t y alrededor de cada punto de la trayectoria medida, ajustar un polinomio de tercer orden entre t - m"t y t + m"t para cada componente x1, x2, x3, donde 2m es el tamaño de la ventana utilizada. Las constantes ci para el polinomio de tipo

están determinados por

donde

y

Luego del filtrado, la posición, velocidad y aceleración xi, ui y ai son obtenidas como

Para la implementación del código se generaron diferentes trayectorias a las que se les añadió ruido de alta frecuencia (Figura 10, panel superior). Seguidamente, se comparó el resultado de derivarlas numéricamente con aquel proveniente del filtrado propuesto (Figura 10, panel central). Se empleó una trayectoria sinusoidal con una frecuencia de 3 Hz, a la que se le agregó ruido de una amplitud del 25% de la amplitud original con un espectro plano entre 0 y 50 Hz. De esta forma, el error en la velocidad calculada a partir de diferencias finitas resultó ser del 205%, mientras que utilizando una ventana de un largo de 6% de la longitud de la curva y filtrándola, se obtuvo un error en la velocidad de sólo el 14%.

Figura 10: Panel superior: Se observa una trayectoria sintética, a la que se le añade ruido sintético y se procede a filtrarla empleando la estrategia propuesta en el presente trabajo. Panel central: Velocidad correspondiente a la trayectoria expuesta en la figura anterior, obtenida numéricamente a partir de diferencias finitas y empleando el filtrado propuesto. Panel inferior: Espectro de la velocidad original, con ruido sintético agregado y luego de ser filtrada por el método propuesto.

Además, se observaron mejoras sustanciales en el espectro de la energía cinética (ver Figura 10, panel inferior). Se encontró que esta ventana de filtrado disminuye especialmente el ruido de alta frecuencia (superior a 6 Hz), haciéndolo particularmente apto para el sistema de forzado empleado, el cual presenta un espectro filtrado en frecuencias en el rango de 0 a 4 Hz. Adicionalmente, en una segunda instancia, se utilizaron dos métodos adicionales de medición, denominados de double y triple shutting. Estos consisten en capturar conjuntos de dos o tres fotografías con una rápida cadencia temporal (en particular se utilizó una frecuencia de 250 Hz) a intervalos regulares (en el presente trabajo se empleó una frecuencia de 80 y 60 Hz, respectivamente). Asimismo, se compararon las velocidades obtenidas con los métodos de double y triple shutting, encontrándose una buena correspondencia entre ambas, lo que convalida el método de filtrado implementado.

VI. Aplicación del sistema PTV 3D MP al estudio de la dinámica de partículas flotadoras

Una vez diseñado, construido y caracterizado el sistema PTV 3D-MP se lo empleó para estudiar la dinámica de flotadores en la superficie libre de un campo de ondas turbulento. Para ello, se hizo uso del montaje descripto en la Sección II y del protocolo de medición que se describe brevemente a continuación. En primer lugar, se realiza la calibración individual de cada cámara del arreglo, en condiciones de tanque de ondas vacío. Luego se llena el tanque con el líquido de trabajo y se introducen las partículas. A continuación se hace uso del sistema de PTV para, a partir de la posición de las partículas flotando en reposo sobre la interfaz, determinar la ubicación de dicho plano. Luego se pone en marcha el sistema de generación de oleaje, y se espera el tiempo necesario para establecer un campo de ondas turbulento completamente desarrollado (dicha transición usualmente implica un forzado de transición de 10 minutos). En estas condiciones, se daba comienzo a la adquisición de imágenes para la medición de trayectorias de partículas. Haciendo uso de este protocolo de medición se llevaron a cabo 56 corridas.
Cabe destacar que cada medición genera aproximadamente 8 Gb de datos, que se corresponden a 2726 pares de imágenes. Dichos datos fueron procesados empleando lo desarrollado en la Sección IV, lo cual dio como resultado las trayectorias 3D de las partículas, cuyo análisis se detalla en la sección subsiguiente.

Trayectorias
La campaña de mediciones que constituye la base de los resultados expuestos en este trabajo experimental consta de un volumen de datos que excede ligeramente los 450 GB de datos, de los cuales se obtuvieron un total de 11391 trayectorias de más de 50 puntos. Se encontró que éstas presentan una alta variabilidad tanto espacial como temporal.
Sus longitudes resultaron muy variadas; observándose trayectorias de hasta 513 puntos para las mediciones realizadas con una amplitud de 20 mm de forzado. La distribución de longitudes de trayectorias medidas se muestra en la Figura 11 en forma de histograma. En dicha figura se observa que la cantidad de trayectorias decrece con la longitud. Esto era esperable, dado que si en un instante de tiempo la partícula deja de ser detectada por el sistema, la trayectoria se da por concluida. De esta forma la distribución de longitudes resultó compatible con una distribución exponencial de longitud de decaimiento 40 y con una media de 98. De esto se interpreta que en cada instante hay una probabilidad de dejar de detectar una partícula de 1/98. Esto, puede deberse a que la misma desaparece del campo de visión de alguna de las cámaras o a que es ocultada por una ola o por otra partícula. Por otro lado, para optimizar la frecuencia de muestreo, se estudio la distancia recorrida por las partículas sobre las imágenes de cada cámara, medida en pixeles. Esto se debe a que si las partículas se desplazaran entre cada snapshot mucho menos de 1 px, su desplazamiento sería comparable con el error de detección y se estarían desperdiciando snapshots. Se analizó la PDF de estos desplazamientos 2D (Figura 12) y se encontró que, en el peor escenario (que corresponde a una amplitud de forzado de 5 mm con una frecuencia de muestreo de 120 Hz) las partículas se desplazaban en promedio 0,94 px. Este resultado muestra que la frecuencia de adquisición resulta adecuada a los propósitos de este estudio.

Figura 11: Histograma de las longitudes de las trayectorias obtenidas y ajuste por medio de una exponencial.

Figura 12: PDF de la distancia recorrida por los flotadores sobre las imágenes de ambas cámaras, para una amplitud de forzado de 5 mm y frecuencia de adquisición de 120 Hz.

Dinámica de clusters de partículas en la superficie libre
Finalmente, se aplicó el sistema de PTV 3D-MP al estudio de la formación de clusters de partículas en el plano XY . Para ello se proyectaron las trayectorias sobre el plano y se realizó, cuadro a cuadro para cada instante de tiempo, un teselado de Delaunay² con la posición de las partículas como vértices (Figura 13, panel izquierdo). Se calcularon luego las áreas de los triángulos, las cuales se normalizaron por la media del cuadro, y se analizó la distribución de las mismas (Figura 14). Se encontró que la función de distribución de probabilidad resulta compatible con una distribución gamma de parámetro a, indicativo de la formación de clusters. Por otro lado, para identificar los clusters, se comparó la PDF medida de los flotadores con una obtenida distribuyendo puntos aleatoriamente en un cuadrado. Se encontró que ambas diferían significativamente, y utilizando el criterio propuesto por ^[19] se pudo determinar un área adimensional crítica de valor Ac = 0, 23. Se pudieron identificar los diferentes clusters en los que se distribuyen las partículas tomando como pertenecientes a un cluster los flotadores que formaran un tringulo de área adimensional menor a un valor crítico, denominado Ac (Figura 13, panel derecho).

Figura 13: Panel izquierdo: Teselado de Delaunay típico para un cuadro, con los flotadores en negro en los vértices. Panel derecho: Identificación de 3 clusters de partículas, cada cluster se indica con un color distinto.

Figura 14: PDF típica de las áreas de los teselados. Se muestran los resultados de las mediciones, su promedio y su ajuste por una distribución gamma de parámetro (adimensional) a.

VII. Conclusiones y Perspectivas

Con el objetivo de obtener una caracterización lagrangiana de flujos de laboratorio se desarrolló el sistema de PTV 3D multipartícula y se lo aplicó al estudio de un sistema que mimicase la dinámica de partículas flotantes en el océano.
Para ello, se diseñó un blanco de calibración y se implementó el método de Zhang para la calibración de las cámaras. También, se desarrollaron e implementaron algoritmos para la corrección de la distorsión radial, la detección de las partículas en las imágenes, la obtención de sus trayectorias 2D, la resolución del problema de la correspondencia mediante geometría epipolar y la triangulación de su posición 3D. Se caracterizó este sistema de medición, por medio de los errores de reproyección y los errores de posicionamiento 3D. Estos resultaron ser de (0, 3 ± 0, 1) px y (0, 3±0, 1) mm respectivamente, confirmando que el método de medición posee resolución submilimétrica, lo que resulta satisfactorio para el seguimiento de las partículas flotadoras. Haciendo uso del montaje de generación de oleaje y del sistema PTV 3D-MP, fue posible recuperar las trayectorias 3D de las partículas. éstas resultaron ser complejas y presentaron alta variabilidad espacial y temporal.
Finalmente, se aplicó la técnica desarrollada a la formación de clusters de partículas en un flujo turbulento. Para ello, se realizaron teselados de Delaunay de las posiciones de las partículas para cada cuadro y se analizó la distribución de las áreas de los triángulos. Se encontró que dichas PDFs eran compatibles con una distribución gamma de parámetro a = (1, 7 ± 0, 3). Esto resultó consistente con la formación de clusters de flotadores, los cuales se pudieron identificar cuadro a cuadro.
Los resultados experimentales presentados en este trabajo muestran la factibilidad de realizar estudios detallados de la dinámica lagrangiana de partículas en condiciones controladas de laboratorio, y ofrecen evidencia clara de la potencialidad que el sistema de PTV 3D-MP desarrollado e implementado en este trabajo presenta para dicha clase de estudio. En particular, a futuro se propone realizar un análisis estadístico de la velocidad y aceleración de las partículas flotadoras en el sistema considerado, analizando el tiempo de vida de los clusters y la dinámica del intercambio de partículas que los integran a lo largo de su vida. Este tipo de estudios permitirán avanzar considerablemente en el entendimiento de fenómenos de transporte de interés tanto en física fundamental como en ciencias atmosféricas y/o oceanográficas.

Los autores dejan constancia de que los resultados mostrados en este trabajo fueron presentados en formato de póster y en presentación oral en la Reunión de División de Fluidos y Plasmas de la Reunión AFA 2017. P.J.C. agradece el financiamiento provisto por la Agencia Nacional de Promoción Científica y Tecnológica a través del PICT 2015-3530. En último lugar, se hace explícito que todos los autores del presente estudio trabajaron por igual en cada una de las etapas del proyecto.

Notas

¹ Remolinos con dimensión de submesoescala identificados como agua mediterránea que se separan de los flujos principales. Estos se caracterizan por la alta salinidad y temperatura que presentan.

² Red de triángulos conexa y convexa que cumple la condición de Delaunay. Esta condición dice que la circunferencia circunscrita de cada triángulo de la red no debe contener ningún vértice de otro triángulo.

Referencias

1. PL Richardson, AS Bower, and Walter Zenk. A census of meddies tracked by floats. Progress in Oceanography, 45(2):209–250, 2000.         [ Links ]

2. Haitao Xu, Alain Pumir, Gregory Falkovich, Eberhard Bodenschatz, Michael Shats, Hua Xia, Nicolas Francois, and Guido Boffetta. Flightcrash events in turbulence. Proc. N. Acad. Sci., 111(21):7558–7563, 2014.         [ Links ]

3. HG Maas, A Gruen, and D Papantoniou. Particle tracking velocimetry in three-dimensional flows. Exp. Fluids, 15(2):133–146, 1993.         [ Links ]

4. NA Malik, Th Dracos, and DA Papantoniou. Particle tracking velocimetry in three-dimensional flows. Exp. Fluids, 15(4-5):279–294, 1993.         [ Links ]

5. Emanuele Trucco and Alessandro Verri. Introductory techniques for 3-D computer vision, volume 201. Prentice Hall Englewood Cliffs, 1998.         [ Links ]

6. Richard Hartley and Peter Sturm. Triangulation. Computer vision and image understanding, 68(2):146–157, 1997.         [ Links ]

7. Zhengyou Zhang. Determining the epipolar geometry and its uncertainty: A review. Int. J. Computer Vision, 27(2):161–195, 1998.         [ Links ]

8. Richard Hartley and Andrew Zisserman. Multiple view geometry in computer vision. Cambridge University Press, 2003.         [ Links ]

9. B. Caprile and V. Torre. Using vanishing points for camera calibration. International Journal of Computer Vision, 4(2):127–139, 1990.         [ Links ]

10. Stephen Maybank and Olivier Faugeras. A theory of self-calibration of a moving camera. International Journal of Computer Vision, 8(2):123–151, 1992.         [ Links ]

11. Richard Hartley. Self-calibration from multiple views with a rotating camera. In Proceedings of the Third European Conference on Computer Vision (Vol. 1), pages 471–478, 1994.         [ Links ]

12. J.Weng, P. Cohen, and M. Herniou. Camera calibration with distortion models and accuracy evaluation. IEEE Transactions on Pattern Analysis and Machine Intelligence, 14(10):965–980, 1992.         [ Links ]

13. Zhengyou Zhang. Flexible camera calibration by viewing a plane from unknown orientations. In Proceedings of the Seventh IEEE International Conference on Computer Vision, volume 1, pages 666–673. IEEE, 1999.         [ Links ]

14. A.R. Conn, N.I.M. Gould, and P.L. Toint. Trust Region Methods. MPS-SIAM Series on Optimization. Society for Industrial and Applied Mathematics, 2000.         [ Links ]

15. Junying Xia, Jiu-long Xiong, Xiaoquan Xu, and Haiyan Qin. A multiscale sub-pixel detector for corners in camera calibration targets. In Int. Conf. on Intelligent Computation Technology and Automation, pages 196–199. IEEE, 2010.         [ Links ]

16. H. Freeman and L. Davis. A corner-finding algorithm for chain-coded curves. IEEE Trans. Comput., 26(3):297–303, 1977.         [ Links ]

17. Jianping Hua and Qingmin Liao. Wavelet-based multiscale corner detection. Signal Processing, 1:341 – 344, 2000.         [ Links ]

18. Beat Luthi. Some aspects of strain, vorticity and material element dynamics as measured with 3D particle tracking velocimetry in a turbulent flow. PhD thesis, Twente University, Netherlands, 2002.         [ Links ]

19. Pablo Gutiérrez and Sébastien Aumaître. Clustering of floaters on the free surface of a turbulent flow: An experimental study. European Journal of Mechanics - B/Fluids, 60(Supplement C):24         [ Links ]