Regístrese
Buscar en
RIBAGUA - Revista Iberoamericana del Agua
Toda la web
Inicio RIBAGUA - Revista Iberoamericana del Agua Simulación directa de turbulencia en corrientes de gravedad con efecto Coriolis
Información de la revista
Vol. 1. Núm. 1.
Páginas 26-37 (Enero 2014)
Compartir
Compartir
Descargar PDF
Más opciones de artículo
Visitas
2371
Vol. 1. Núm. 1.
Páginas 26-37 (Enero 2014)
DOI: 10.1016/S2386-3781(15)30005-0
Open Access
Simulación directa de turbulencia en corrientes de gravedad con efecto Coriolis
Direct numerical simulations of gravity currents with Coriolis effect
Visitas
2371
J.S. Salinas, M.I. Cantero
Autor para correspondencia
mcantero@cab.cnea.gov.ar

Autor para correspondencia. Correo electrónico:
, E.A. Dari
Instituto Balseiro, San Carlos de Bariloche, Río Negro, Argentina
Consejo Nacional de Investigaciones Científicas y Técnicas, San Carlos de Bariloche, Río Negro, Argentina
Centro Atómico Bariloche, Comisión Nacional de Energía Atómica, San Carlos de Bariloche, Río Negro, Argentina
Este artículo ha recibido
2371
Visitas

Under a Creative Commons license
Información del artículo
Resumen
Texto completo
Bibliografía
Descargar PDF
Estadísticas
Figuras (8)
Mostrar másMostrar menos
Resumen

Las corrientes de gravedad son flujos generados por gradientes de presión horizontales resultantes del efecto de la gravedad sobre fluidos de diferente densidad. Cuando ocurren en la naturaleza, estas corrientes tienen un fuerte comportamiento no lineal y presentan un amplio rango de escalas temporales y espaciales. Adicionalmente, la rotación de la Tierra eleva el nivel de complejidad para el estudio de las corrientes de gravedad debido al efecto de la fuerza de Coriolis. En este trabajo se abordó el estudio de las corrientes de gravedad en geometría plana en rotación mediante simulación directa de turbulencia. Las simulaciones realizadas permitieron un análisis detallado de la evolución del flujo, tanto en los parámetros macroscópicos como en la estructura de la turbulencia. Las simulaciones se realizaron mediante un código pseudoespectral que emplea expansiones de Fourier en las dos direcciones horizontales y expansiones de Chebyshev en la dirección vertical. En este trabajo se documentó en detalle el código de cálculo desarrollado y se validó comparando una simulación con observaciones experimentales y predicciones teóricas. Se realizaron dos simulaciones con distintas velocidades de rotación. Se comprobó que dicha rotación restringe el desarrollo de la corriente en la dirección de propagación principal y se observaron oscilaciones en la posición del frente cuya frecuencia es proporcional a la velocidad de rotación. Además, se observaron estructuras turbulentas de tipo Kelvin-Helmholtz verticales en el frente de la corriente producidas por la rotación del sistema.

Palabras clave:
Coriolis
Corrientes de gravedad
Rotación
Simulación directa de turbulencia
Abstract

Gravity currents are flows generated by horizontal pressure gradients resulting from the effect of gravity on fluids of different density. When they occur in nature, gravity currents have a strong nonlinear behavior and have a wide range of temporal and spatial scales. In addition, the rotation of the earth raises the level of complexity of gravity currents due to the effect of the Coriolis force. This work addresses rotating gravity currents in planar geometry by direct numerical simulation (DNS). The simulations allow for a detailed analysis of the flow development, macroscopic parameters of the flow and turbulence structure. Simulations were performed using a pseudospectral code that uses Fourier expansions in the two horizontal directions and Chebyshev expansions in the vertical direction. This work documents in detail the code developed and the validation performed by comparing a simulation with experimental observations and theoretical predictions. This work also reports on two simulations with different rotation speeds. It was found that the flow rotation restricts the development of the current in the propagating direction, and induces oscillations in the front position. The frequency of these oscillations varies linearly with the rotation speed. Finally, this work also reports on Kelvin-Helmholtz-like turbulent structures at the front of the current produced by the rotation of the system.

Keywords:
Gravity currents
Direct numerical simulations
DNS
Coriolis
Rotational flows
Texto completo
1Introducción

Las corrientes de gravedad (también llamadas corrientes de densidad) son flujos generados por gradientes de presión horizontales resultantes del efecto de la gravedad sobre fluidos de diferente densidad 1. Estos flujos se manifiestan como una corriente horizontal de fluido liviano por encima de un fluido pesado, o como una corriente de fluido pesado debajo de un fluido liviano. Estos fenómenos naturales son importantes en el transporte de masa, cantidad de movimiento y energía; a grandes escalas, que derivan en altos números de Reynolds, son altamente turbulentos y presentan un amplio rango de escalas temporales y espaciales.

Las corrientes de gravedad están presentes en la naturaleza en muchas situaciones y, en muchos de los casos de interés, son producidas por pequeñas diferencias de densidad. En la atmósfera, por ejemplo, la mayoría de las ráfagas intensas asociadas con las tormentas son causadas por la llegada de una enorme corriente de gravedad de aire más frío. Muchas veces estas corrientes se pueden identificar por la presencia de partículas de arena y polvo que han sido resuspendidas de la superficie terrestre por los fuertes vientos. Las avalanchas de nieve en las montañas son corrientes de gravedad en las cuales la diferencia de densidad se debe a la suspensión de partículas de nieve en el aire. Debido a los inmensos daños que pueden provocar estos fenómenos, en la actualidad existen centros de investigación dedicados especialmente al estudio de este tipo particular de corrientes de gravedad. Otro ejemplo de corrientes de gravedad es la marea negra. Un derrame de petróleo de un barco produce una corriente de gravedad de este fluido sobre la superficie marítima.

En las últimas décadas se han realizado gran cantidad de trabajos de simulación. Se han llevado a cabo estudios basados en simulaciones numéricas bidimensionales y tridimensionales con el objetivo de explorar la dinámica de las corrientes de gravedad2–11. Se han realizado simulaciones tanto en configuraciones planas 7,12-15 como cilíndricas 12, 13. Simulaciones tridimensionales de alta resolución fueron realizadas por Cantero et al 16 para corrientes de gravedad cilíndricas, y por Cantero et al 17 para corrientes de gravedad planas.

Trabajos como los de Cenedese y Adduce 18 y Marino et al 19 estudiaron las corrientes de gravedad por medio de experimentos de laboratorio.

Por su parte, los trabajos de Fannelop y Waldman 20 y Hoult 21 estudiaron el comportamiento de corrientes de gravedad mediante modelos teóricos simples. Estos modelos predicen la velocidad del frente durante las fases de hundimiento, inercial y viscosa. En la fase de hundimiento, el frente se mueve a una velocidad prácticamente constante (uF ≈ cte). La teoría hidráulica para corrientes de gravedad sin efectos de Coriolis propone uF ≈ 0,5 22, 23. Basándose en experimentos de laboratorio, Huppert y Simpson 24 reportaron uF ≈ 0,45, y Cantero et al 12, basado en simulaciones directas de turbulencias y resultados experimentales, ha reportado uF = 0,45. Tras la fase de hundimiento, se produce la transición a la fase inercial. El comportamiento asintótico de la corriente en la fase inercial se ha establecido como 21,24–26:

Aquí, uF, x0yh0 son la velocidad del frente en la dirección de propagación de la corriente, la altura inicial de la corriente y el largo inicial, respectivamente. La diferencia entre las teorías se encuentra en la constante ξp. En nuestro estudio utilizamos el valor ξp = 1,47 propuesto por Cantero et al 12 en base al análisis de una gran cantidad de datos experimentales.

Realizando el balance entre la fuerza boyante y las fuerzas viscosas a lo largo de la interfaz entre el fluido de la corriente, el fluido ambiente y una superficie rígida horizontal, Huppert 27 obtuvo la siguiente solución para la fase viscosa:

En este trabajo utilizamos el valor ξpHp= 3,2 propuesto por Cantero et al 12.

La rotación de la Tierra eleva el nivel de complejidad para el estudio de las corrientes de gravedad debido a que la fuerza de Coriolis modifica drásticamente el desarrollo de la corriente. A diferencia de lo que ocurre con las corrientes de gravedad sin efectos de rotación, en las corrientes rotantes la literatura especializada es mucho más escasa. El trabajo más completo sobre los efectos de rotación en corrientes de gravedad es probablemente el de Hallworth et al 28, quienes reportaron resultados experimentales y numéricos (en dos dimensiones) de corrientes de gravedad cilindricas en sistemas en rotación, focalizándose principalmente en el comportamiento general de las corrientes y en las características cuantitativas macroscópicas de lasmismas. No se reportaron trabajos similares para corrientes de gravedad en geometría plana.

En este trabajo se abordó el estudio de las corrientes de gravedad en geometría plana en rotación mediante simulación directa de turbulencia. Las simulaciones se realizaron con un código pseudoespectral que utiliza expansiones de Fourier en las dos direcciones horizontales, y expansiones de Chebyshev en la dirección vertical. El trabajo se focalizó principalmente en documentar de manera detallada el modelo numérico utilizado, y en validarlo mediante la comparación de una simulación tridimensional de una corriente de gravedad sin rotación con datos experimentales. Para el análisis de las corrientes de gravedad bajo los efectos de la rotación se realizaron dos simulaciones tridimensionales con una resolución de 22 millones de puntos de grilla, número de Reynolds Re = 4.000, número de Schmidt Sc = 1 y parámetros de Coriolis C = 0,15 y C = 0,25. Este último parámetro se define como la relación entre la fuerza de Coriolis y la inercial. El detalle alcanzado en las simulaciones permitió un análisis detallado de la dinámica del flujo. Además, los valores elegidos para el parámetro de Coriolis fueron los adecuados para visualizar de forma evidente el efecto de la fuerza de Coriolis sobre la corriente. Se presentaron resultados sobre la máxima distancia de propagación del frente, la frecuencia de las oscilaciones del frente, el estado cuasi-estacionario de la corriente y las estructuras turbulentas presentes en el flujo.

2Modelo matemático

El problema que consideramos se presenta de forma esquemática en la figura 1. Consiste en un dominio rectangular rotando a una velocidad angular constante ΩZ alrededor del eje vertical z que pasa por el centro del dominio. Un fluido pesado de densidad pv inicialmente en una región rectangular de dimensiones 2x0×Ly×H (región gris sombreada en la figura 1), se encuentra separado del fluido ambiente de densidad ρ01 > ρ0). Los dos fluidos se encuentran inicialmente en co-rotación con el dominio. Al comienzo de la simulación, el fluido más denso se libera produciendo un flujo que se propaga principalmente en la dirección x. Para nuestras simulaciones utilizamos x0= H.

Figura 1.

Representación esquemática del fenómeno físico. El dominio de dimensiones 2Lx × Ly × H se encuentra rotando a una velocidad angular Ωz alrededor del eje vertical z que pasa por el centro del dominio. El fluido de densidad r1, en la región rectangular de dimensiones 2x0 × Ly × H (región gris sombreada), se encuentra separado del fluido ambiente de densidad ρ010). Los dos fluidos se encuentran inicialmente en co-rotación con el dominio. Al comienzo de la simulación el fluido más denso se libera, produciendo un flujo que se propaga principalmente en la dirección x.

(0,07MB).

Se consideran flujos en los que la diferencia de densidad es lo suficientemente pequeña como para que la aproximación de Boussinesq sea válida. Bajo estas circunstancias, las ecuaciones de conservación adimensionales para un fluido newtoniano en régimen incompresible en un sistema de referencia rotante son

En estas ecuaciones ũ={ũx,ũy,ũz}={ũ,v˜,w˜} es el campo de velocidades adimensional, p˜ la presión adimensional, y ρ˜ la densidad adimensional. El vector unitario gˆ apunta en la dirección de la fuerza de gravedad. Los parámetros adimensionales en las ecuaciones (3)–(5) son el número de Reynolds Re, de Schmidt Sc, y el parámetro de Coriolis C˜ definidos como:

donde v es la viscosidad cinemática y κ la difusividad molecular.

Para la adimensionalización de las ecuaciones (3)–(5) se emplea como escala de longitud la altura del dominio H, escala de velocidad U=β|g|H con β=(ρ10)/ρ0, y escala de tiempo T = H/U. La densidad adimensional ρ˜ se define como:

donde ρ=ρ0(1+βρ˜). La densidad adimensional se encuentra en el intervalo 0≤ρ˜≤1, siendo ρ=ρ0 cuando ρ˜=0 y ρ=ρ1 cuando ρ˜=1.

Las ecuaciones de conservación (3)–(5) se resuelven en un dominio rectangular de dimensiones 2L˜x×L˜y×H˜. En las direcciones horizontales x˜ e y˜ se aplican condiciones de contorno periódicas para todas las variables. Esto significa que las variables del flujo son iguales en las paredes verticales x˜=−L˜x y x˜=L˜x, así como en y˜=0 y y˜=−L˜y.

En la dirección vertical z˜ se aplican condiciones de contorno de Neumann paraρ˜ y de “libre deslizamiento”/“no deslizamiento” para ũ.

3Modelo numérico

Por simplicidad de notación, en esta sección se omitirá el uso del ~ para indicar si una variable es adimensional.

3.1Geometría

El dominio está compuesto por dos paredes paralelas separadas una distancia H en la dirección vertical z (Fig. 2). El modelo asume que el flujo es periódico en las direcciones x e y. En la dirección vertical, el flujo es no homogéneo por la presencia de las paredes. El dominio de cálculo tiene dimensiones 2Lx×Ly×H en las direcciones x, y y z, respectivamente.

Figura 2.

Sección del dominio periódico de dimensiones 2Lx × Ly × H.

(0,06MB).
3.2Discretización temporal

Para resolver las ecuaciones de conservación (3)-(5) se utiliza un método de pasos fraccionados 29. El campo de velocidades del flujo es avanzado desde el instante de tiempo t(n) al t(n+1) en dos pasos sucesivos. Primero se utiliza una ecuación de advección-difusión para avanzar desde el instante t(n) al paso intermedio t(n*). Esta ecuación consiste en la ecuación (3) sin el término del gradiente de presión. Las condiciones de contorno para este paso intermedio serán discutidas en la siguiente sección. Luego de este paso intermedio (*) se utiliza un paso de corrección de presión para avanzar al instante de tiempo t(n+1).

La discretización temporal del paso de advección-difusión y de la ecuación de transporte del escalar se logra utilizando un esquema mixto Runge-Kutta de tercer orden y Crank-Nicolson 30. El esquema se lleva a cabo en tres etapas. El paso de tiempo de t(n) a t(n+1)t) es fraccionado en tres pasos más pequeños, con corrección de presión al final de cada subpaso de Runge-Kutta. En cada subpaso de Runge-Kutta, la corrección de la presión se lleva a cabo en dos etapas. En la primera se resuelve una ecuación de Poisson para la presión y, en la segunda, se actualiza el campo de velocidades.

Para simplificar la notación se definen los operadores de advección y difusión. El operador de difusión para la ecuación de conservación de momento es:

Para el campo de velocidades se utilizan dos operadores de advección:

Esto corresponde a las formas convectivas y divergentes del operador de advección. Alternando entre estos dos operadores de advección para sucesivos pasos de tiempos, se obtienen resultados similares a los que se ven con formulaciones antisimétricas (skew-symmetric, en inglés) de las ecuaciones, pero con la ventaja de que el número de derivadas que se necesita evaluar por paso de tiempo se reduce a la mitad. De aquí en adelante, el operador de advección será referido como A.

Además, se definen los operadores de advección y difusión para la ecuación de transporte del escalar. El operador de advección es:

y el de difusión:

El vector de fuerzas se define como:

Para mostrar un paso de tiempo completo se define un número de etapa m de Runge-Kutta, un coeficiente para el primer término no lineal cnl1, un coeficiente para el segundo término no lineal cnl2 y un coeficiente para el término difusivo cd:

Los campos en la etapa de nivel 0 son iguales a los del tiempo t(n):

Las ecuaciones empleadas en un subpaso de Runge-Kutta son:

Primero se calcula el campo de velocidades auxiliar (m*) con la ecuación (24), y el campo escalar de la etapa (m) con la ecuación (25). Luego, se calcula la presión correspondiente resolviendo la ecuación de Poisson (26) que hace uso del campo de velocidad auxiliar. Paso seguido, se corrige el campo de velocidades con la ecuación (27), obteniéndose el campo de velocidades de la etapa (m) con divergencia nula. Una vez que se completa la última etapa, los campos son actualizados:

y el proceso puede comenzar nuevamente.

3.3Discretización espacial

En las direcciones homogéneas horizontales x e y se utilizan expansiones de Fourier, mientras que en la dirección inhomogénea vertical z se utilizan expansiones de Chevyshev. La expansión utilizada para discretizar una variable u es:

donde Tkll se refiere al kl-ésimo polinomio de Chevyshev, kx=πkjLx y ky=2πkkLy 2son los números de onda, y uˆ (kj,kk,kl,t) es la amplitud del modo correspondiente.

En la dirección vertical z se utilizan puntos de cuadratura de Gauss-Lobatto, los cuales proveen una mayor resolución cerca de las paredes. La posición de los puntos de grilla en las direcciones horizontales x e y es equiespaciada. Las posiciones de los puntos de grilla son

donde nx, ny y n son la cantidad de puntos de grilla en las direcciones x, y y z, respectivamente.

Derivadas espectralmente precisas se obtienen utilizando técnicas de transformadas rápidas de Fourier en las direcciones horizontales y técnicas de multiplicación matricial en la dirección vertical.

Un paso de tiempo se completa utilizando las ecuaciones (20) a (29). Las ecuaciones (22) a (27) deben ser resueltas para cada una de las tres etapas del esquema. Las ecuaciones (24) y (25) se pueden escribir como ecuaciones de Helmholtz para las tres componentes de la velocidad y el campo escalar. La componente j-ésima de la ecuación de conservación de momento puede reescribirse como:

El lado derecho de la ecuación (34) puede ser calculado con información conocida y, para simplificar el análisis, se escribirá como rhˆsj(m*). La transformada de Fourier en las direcciones x e y de la ecuación (34) es:

donde “^” indica coeficiente de Fourier. De manera similar se puede proceder para la ecuación (25) para el campo escalar:

La discretización en la dirección z se expresa en términos de multiplicación matricial 31 combinando una transformación discreta de Chevyshev, diferenciación recursiva y una transformación discreta de Chevyshev inversa. De manera simbólica, la derivada primera con respecto a z se puede escribir como:

La matriz D══1() es antisimétrica de n filas y n columnas. El operador para realizar la segunda derivada parcial con respecto a z se formula simplemente multiplicando dos operadores de derivada primera 31.

El desarrollo que sigue se realiza para una variable cˆ (en el espacio de Fourier) que puede representar el campo escalar ρ˜ o una de las tres componentes del campo vectorial de velocidades (uˆ,vˆ o wˆ ).

3.4Tratamiento de las condiciones de contorno

Las ecuaciones (35) y (36) son solo válidas para el interior del dominio. El sistema de ecuaciones se completa con condiciones de contorno. En las direcciones horizontales, las condiciones de contorno son periódicas y son capturadas directamente por las expansiones de Fourier. En esta sección se detalla la implementación de las condiciones de contorno en la dirección vertical, para la cual se especifican condiciones de contorno mixtas. Condiciones tipo Dirichlet o Neumann son casos particulares de las condiciones mixtas. El tratamiento presentado aquí se aplica a cada subpaso de Runge-Kutta.

Una vez transformado al espacio de Fourier, el sistema de ecuaciones completo que se ha de resolver es:

donde el subíndice “b” se refiere al borde inferior del dominio y “t” al borde superior. El coeficiente α de la ecuación (39) puede adquirir distintos valores según qué variable represente cˆ:

Las constantes αb, αt, βb y βt pueden tomar cualquier valor numérico. En la forma más general, los valores de las condiciones de contorno pueden variar punto a punto en las paredes, es decir γˆb=γˆbkx,ky y γˆt=γˆtkx,ky. En el caso que cˆ=uˆ o cˆ=vˆ:

donde at = 0, βt = 1 y yˆt=0 representa la condición de contorno de “libre deslizamiento”, mientras que at = 1, βt = 0 y yˆt=0 la condición de “no deslizamiento”. Para el caso cˆ=wˆ:

la cual representa la condición de contorno de “no penetración”.

El sistema de ecuaciones (39)–(41) puede presentarse en su forma discreta:

que en forma matricial es:

donde:

Para obtener la solución a este sistema de ecuaciones, se descompone en tres subsistemas: la ecuación de la primer fila (k = 1), el subsistema de ecuaciones k = 2, …, n – 1 y la ecuación de la última fila (k = n). Tras el trabajo algebraico, se obtiene una ecuación para el valor del campo cˆ en el borde inferior del dominio (cˆn) y otra para el valor en el borde superior del dominio (cˆ1):

Empleando (52) y (53), el sistema (51) se puede reescribir en forma reducida como:

donde la matriz B══ (de dimensión n – 2 × n – 2) y los vectores cˆ* y Nˆ (de dimensión n – 2) son:

con i = 2, ..., n – 1 y k = 2, …, n – 1.

Utilizando las ecuaciones (52) y (53) junto con el sistema de ecuaciones (54) se obtienen los valores de cˆ para nuestro problema con condiciones de contorno mixtas. El procedimiento de resolución se explicará en la sección 3.5.

3.4.1Condiciones de contorno para el campo de velocidades intermedio

El esquema de paso fraccionado demanda la especificación de una condición de contorno para el campo de velocidades intermedio. Esta condición de contorno debe ser elegida de modo tal que la velocidad al final del paso temporal satisfaga las condiciones de contorno necesarias. El esquema de paso fraccionado también demanda la especificación de una condición de contorno adecuada para el paso de presión. Aunque no existe una condición de contorno natural para la presión, es de práctica usual el uso de una condición de Neumann pura. Detalles sobre su implementación se pueden encontrar en Deville et al 29.

Siendo t˘ y n˘ los versores en las direcciones tangencial y normal a la pared, respectivamente, las componentes tangencial y normal del paso de corrección de presión (27) son:

Debido a que ▽pˆ(m) no ha sido calculado cuando se necesitan las condiciones de contorno de velocidad intermedias, se debe utilizar una aproximación de ▽pˆ(m). Para esta aproximación utilizamos una serie de Taylor alrededor de t = t(m – 1), donde la derivada temporal de ▽pˆ(m) se puede aproximar utilizando diferencias finitas hacia atrás, obteniéndose una aproximación de segundo orden para ▽pˆ(m).

Reordenando términos en la ecuación (56), las condiciones de contorno para las velocidades intermedias en las direcciones horizontales (direcciones tangenciales a la pared) son:

En el caso particular que el segundo término del lado derecho de las ecuaciones (58) y (59) sea cero:

se está imponiendo una condición de “no deslizamiento” en las paredes.

En el caso de utilizar condiciones de contorno del tipo Neumann para el campo de velocidades, la condición de contorno para el paso intermedio (m*) en las direcciones horizontales (direcciones tangenciales a la pared) es:

En el caso particular donde el segundo término del lado derecho de las ecuaciones (61) y (62) sea cero, estamos imponiendo una condición de “libre deslizamiento” en las paredes.

Para la componente vertical (componente normal a la pared) de la velocidad intermedia se debe tener en cuenta que el paso de corrección de presión se resuelve con una condición de contorno de Neumann pura:

En este caso, la condición de contorno para la velocidad vertical (componente normal a la pared) intermedia es simplemente cero:

3.5Procedimiento de resolución

Para cada paso de Runge-Kutta y combinación de números de onda horizontales se debe resolver el sistema de ecuaciones (54) junto con las ecuaciones (52) y (53) para las tres componentes de velocidad intermedia y el campo escalar. Además, se debe resolver la ecuación de Poisson (26) para la presión. Esto quiere decir que, para cada paso de tiempo completo, se deben resolver 3 × 5 × 2 × nx × ny / 2 veces las ecuaciones de Helmholtz o Poisson (tres pasos de Runge-Kutta por cinco campos escalares complejos en nx × ny / 2 combinaciones de números de onda). Por ejemplo, para una grilla de resolución 256 × 256 × 257 se deben resolver 983.040 veces por paso de tiempo sistemas de ecuaciones de 255 × 255 para las ecuaciones de Helmholtz o Poisson. Claramente, se hace necesaria una técnica para resolver estas ecuaciones de manera veloz. En su forma actual, las ecuaciones requieren una inversión matricial para cada combinación de números de onda. Diagonalizando la matriz B══ (ver ecuación [54]), es posible eliminar la dependencia de la inversión con respecto a los números de onda y, por ende, será posible invertir la matriz una sola vez en todo el cálculo. Definimos la matriz modal T══ de la matriz B══, T══−1 la inversa de la matriz modal, λ el vector cuyos elementos son los autovalores de B══ y Λ══ la matriz diagonal de dimensiones n – 2 × n – 2 cuyos elementos diagonales son los autovalores λ2, .... λn-1 de B══. La matriz B══ puede ser descompuesta como:

De esta forma, la ecuación (54) puede ser reescrita como:

El coeficiente α (véanse las ecuaciones [(42)] y [43]) puede ser combinado con Λ══. Sea I══ una matriz identidad de dimensiones n – 2 × n – 2, la ecuación (66) puede invertirse para dar una expresión para cˆ*(m*):

donde:

Una vez calculado cˆ*(m*), se pueden calcular cˆn(m*) y cˆ1(m*) con las ecuaciones (52) y (53), respectivamente.

En la figura 3 se presenta el diagrama de flujo de los pasos requeridos para el avance de un paso temporal Δt, para un par de números de onda (kx, ky). Cabe destacar que la matriz B══ se obtiene antes de comenzar el bucle temporal presentado en la figura 3.

Figura 3.

Diagrama de flujo de los pasos requeridos para el avance de un paso temporal Δt, para un par de números de onda (kx, ky). Las siglas “C.C.” se refieren a las condiciones de contorno.

(0,3MB).

Más detalles sobre la implementación del código y su paralelización híbrida con OpenMP-MPI pueden encontrarse en Salinas 32.

4Resultados

Con el código de cálculo presentado en la sección anterior, se realizaron tres simulaciones de corrientes de gravedad planas en un canal periódico. La primera simulación, con un parámetro de Coriolis C = 0, se utilizó para validar el código. La segunda y la tercera simulaciones, con C = 0,15 y C = 0,25, respectivamente, se analizaron para determinar el efecto de la rotación en las corrientes de gravedad planas. El número de Reynolds es Re = 4.000 y el de Schmidt es Sc = 1 para todas las simulaciones. Para la simulación (1) se utilizó un dominio de dimensiones 2Lx × Ly × H = 25 × 1 × 1 con una resolución de nx × ny × n = 1.024 × 80 × 128. Para las simulaciones (2) y (3) se utilizó un dominio de dimensiones 2Lx × Ly × H = 18 × 1 × 1 con una resolución de nx × ny × n = 1.536 × 80 × 180. Se utilizaron condiciones de contorno de “no deslizamiento” en las paredes superior e inferior del dominio para el campo de velocidades en la simulación (1), mientras que para las simulaciones (2) y (3) se utilizaron condiciones de contorno de “libre deslizamiento”. Para el campo escalar r se utilizó la condición de contorno de gradiente nulo en las paredes superior e inferior, para todas las simulaciones. El paso de tiempo Δt utilizado en cada iteración se eligió de forma que se cumpla la condición CFL< 0,7, donde CFL es el número de Courant-Friedrich-Levy.

Para rastrear la posición media del frente de la corriente x¯F se adoptó la definición presentada por Shin et al 23 y Marino et al 19 de la altura equivalente local, la cual se promedió a lo largo de la dirección transversal y, definiéndose de esta forma la altura equivalente local media:

Las variables con una barra deben entenderse como cantidades adimensionales promediadas en la dirección y. La posición media del frente x¯F se define como la ubicación donde la altura equivalente local media h¯ se hace más pequeña que un umbral δ. La velocidad del frente media se determinó como:

4.1Validación del código de cálculo

Para validar el código descrito en la sección anterior se realizó la simulación (1) con un parámetro de Coriolis C = 0. La figura 4 muestra la velocidad del frente media u¯F en función del tiempo para la simulación (1), junto con datos experimentales reportados por Marino et al 19 para Re = 6.360 y predicciones teóricas de las fases de hundimiento, inercial 24 y viscosa para Re = 4.000 27. La velocidad del frente media en la fase de hundimiento se encuentra de acuerdo con los experimentos de laboratorio presentados por Marino et al 19 y con valores reportados por Huppert y Simpson 24 y Cantero et al 12. Luego de esta fase, la velocidad del frente disminuye rápidamente, siguiendo la predicción teórica realizada por Huppert 27 para la fase viscosa.

Figura 4.

Velocidad del frente media u¯F en función del tiempo para la simulación (1). La figura incluye los datos experimentales de Marino et al19 para Re = 6.360 y las predicciones teóricas de las fases de hundimiento, inercial y viscosa para Re = 4.000.

(0,11MB).
4.2Efecto Coriolis en las corrientes de gravedad

Para visualizar el efecto de la rotación en las corrientes de gravedad planas se realizaron dos simulaciones: la simulación (2) con un parámetro de Coriolis C = 0,15 y la simulación (3) con C = 0,25. El número de Reynolds es Re = 4.000 y el de Schmidt es Sc = 1 para ambas simulaciones.

Definimos la función de densidad media promediada en la dirección transversal a la corriente como:

En la figura 5 se muestra la función densidad media ρ¯ para distintos instantes de tiempo para la simulación (2). Como puede verse en esta figura, la corriente desarrolló inicialmente la estructura de cabeza y cuerpo presentando vórtices tipo Kelvin-Helmholtz en la interfaz superior de la misma (véase t = 5 en la figura 5). Esta estructura es idéntica a la que se visualizó en corrientes de gravedad sin rotación17. El efecto de la rotación se hizo presente para tiempos mayores, modificando drásticamente el desarrollo de la corriente. Se modificó la estructura de cabeza y cuerpo, y la estructura de vórtices en la interfaz superior (véase t = 10 en la figura 5). Más aún, la rotación del sistema produjo que la corriente detenga su avance, llegando a un estado cuasi-estacionario (véase fig. 5 para t = 100).

Figura 5.

Densidad media ρ¯ para diferentes instantes de tiempo, para la simulación (2).

(0,12MB).

La figura 6 muestra la posición media del frente con respecto a la posición inicial en función del tiempo para diferentes umbrales δ de h¯, para las simulaciones (2) y (3) con C = 0,15 y C = 0,25, respectivamente. Se pudo observar la restricción en la propagación del frente y las oscilaciones de su posición producidas por la rotación del sistema. Después de un cierto tiempo, la propagación de la corriente en rotación se desaceleró hasta que la velocidad del frente u¯F cambió de signo y el frente comenzó a moverse hacia el interior del dominio. Después de un período de tiempo, la velocidad del frente cambió otra vez de signo y la corriente comenzó a propagarse nuevamente hacia afuera. Este movimiento oscilatorio se observó varias veces después de esta primera oscilación, siendo la máxima distancia de propagación del frente en la primera oscilación x¯max = 5,6 y x¯max = 3,8, para las simulaciones (2) y (3), respectivamente. Como podemos ver en la figura 6, x¯max es influenciado por la velocidad de rotación del sistema y, por tanto, por el parámetro de Coriolis. Con el incremento de este último se produce la disminución de la propagación del frente. Se llegó al estado cuasi-estacionario antes en la simulación (3) (con un mayor parámetro de Coriolis) que en la simulación (2). La máxima distancia de propagación de las oscilaciones varió solamente en un 3% en todo nuestro dominio temporal, por lo que se la consideró constante. Para la determinación de la distancia máxima de propagación del frente x¯max se utilizó el umbral δ = 0,01.

Figura 6.

Posición del frente relativa a la posición inicial x¯F−x¯0 en función del tiempo t para diferentes umbrales δ de h¯, desde δ = 0,01 a δ = 0,1, para: a) Simulación 2; b) Simulación 3.

(0,2MB).

El carácter oscilatorio del flujo puede explicarse del siguiente modo. Para t = 5, el fluido más denso se desplazó en la dirección longitudinal positiva x, lo que produjo una fuerza de Coriolis en la dirección transversal negativa y y, por tanto, un aumento en la velocidad en este sentido. Esta velocidad uy también produjo una fuerza de Coriolis, pero en la dirección longitudinal negativa. Cuando esta fuerza superó la fuerza boyante en la dirección opuesta, provocó que el fluido más denso comience a moverse en la dirección longitudinal negativa x. Este proceso se repite generando el carácter oscilatorio del flujo que muestra la figura 6.

El período medio de las oscilaciones T¯p se define como el intervalo de tiempo promedio en el que las sucesivas oscilaciones hacia afuera alcanzan una determinada distancia x¯. La frecuencia de las oscilaciones ω¯p se define como:

La frecuencia de las oscilaciones del frente es ω¯p = 0,316 y ω¯p = 0,529, para C = 0,15 y C = 0,25, respectivamente. Se logró encontrar una relación lineal entre el parámetro de Coriolis y la frecuencia de las oscilaciones, siendo esta relación ω¯p = 2,1C. Resultados similares fueron reportados por Hallworth et al 28 con experimentos de laboratorio de corrientes de gravedad en rotación.

En la figura 7 se presenta una visualización compuesta de una iso-superficie de densidad r = 0,05 junto a dos cortes del campo de densidad: uno para el plano y = 0,3 y otro para el plano z = 0,05, para distintos instantes de tiempo, para la simulación (2). En las primeras etapas de desarrollo del flujo (t = 2,5), el frente se caracterizó por ser principalmente bidimensional y la interfaz entre el fluido ambiente y la corriente se vio modificada solamente por los núcleos de vórtices transversales tipo Kelvin-Helmholtz en la región del cuerpo de la corriente. En t = 5,75, se pudieron ver los vórtices transversales tipo Kelvin-Helmholtz en el cuerpo de la corriente, cuya existencia es independiente de la presencia de la fuerza de Coriolis. Sin embargo, para este mismo instante de tiempo, se pudo encontrar un nuevo grupo de vórtices tipo Kelvin-Helmholtz verticales que aparece solamente en corrientes de gravedad en rotación. Estos vórtices Kelvin-Helmholtz verticales se produjeron por el esfuerzo de corte en la dirección transversal y entre el fluido ambiente y el de la corriente, mientras que los vórtices Kelvin-Helmholtz transversales en el cuerpo de la corriente son producidos por el esfuerzo de corte en la dirección longitudinal x. En la figura 8 presentamos una vista superior del corte del campo de densidades en z = 0,05 para t = 5,75, donde se puede ver el efecto de los vórtices verticales Kelvin-Helmholtz en el frente. Estructuras similares fueron observadas por Ungarish y Huppert 33 en experimentos de laboratorio, mientras que Hallworth et al 28 mencionaron la existencia de algún tipo de inestabilidad en el frente, aunque no las describieron en detalle. Después de t = 5,75, el flujo se caracterizó por ser altamente tridimensional en el frente, el cuerpo y la cola de la corriente (Fig. 7).

Figura 7.

Estructura tridimensional de la interfaz entre la corriente y el fluido ambiente visualizada mediante una iso-superficie de densidad ρ = 0,05 para distintos instantes de tiempo. En esta figura también se muestran dos cortes del campo de densidad: uno para el plano y = 0,3 y otro para el plano z = 0,05.

(0,15MB).
Figura 8.

Vista en detalle del efecto de las estructuras turbulentas presentes en el frente visualizada mediante un corte del campo de densidades en z = 0,05, para t = 5,75.

(0,05MB).
5Conclusiones

Se analizó el efecto de la fuerza de Coriolis en corrientes de gravedad planas mediante simulación directa de turbulencia. Se explicó en detalle la implementación del código pseudoespectral utilizado y se presentó la validación del mismo mediante la comparación de una simulación de una corriente de gravedad plana sin efecto de Coriolis con observaciones experimentales y predicciones teóricas.

En las dos simulaciones con efecto de Coriolis se pudo observar un comportamiento oscilatorio del flujo. Esto fue determinado mediante el seguimiento de la posición media del frente x¯F. Producto de los efectos de la fuerza de Coriolis, la corriente llegó a un estado cuasi-estacionario (en nuestro dominio temporal) en donde las oscilaciones del frente se mantuvieron aunque disminuyendo levemente su amplitud. La máxima distancia alcanzada por el frente fue x¯max = 5,6 y x¯max = 3,8 para C = 0,15 y C = 0,25, respectivamente. La frecuencia de las oscilaciones en la posición del frente fue ω¯p = 0,316 y ω¯p = 0,529, para C = 0,15 y C = 0,25, respectivamente. Se logró determinar una relación lineal entre la frecuencia de estas oscilaciones y el parámetro de Coriolis, siendo esta ω¯p = 2,1C. Resultados similares fueron reportados por Hallworth et al28 con experimentos de laboratorio de corrientes de gravedad en rotación.

Se identificaron las estructuras de la turbulencia presentes en las corrientes de gravedad planas en rotación. Los vórtice transversales Kelvin-Helmholtz fueron el resultado del esfuerzo de corte en la dirección longitudinal x entre el fluido de la corriente y del ambiente. Por otro lado, los vórtices verticales Kelvin-Helmholtz en el frente de la corriente fueron producidos por el esfuerzo de corte en la dirección transversal y entre los dos fluidos, producto del efecto de la rotación del sistema. Estructuras similares fueron observadas por Ungarish y Huppert33 en experimentos de laboratorio.

Bibliografía
[1]
M. García.
Turbidity currents.
Encyclopedia of Earth System Science. 4, pp. 399-408
[2]
K. Droegemeier, R. Wilhelmson.
Kelvin-Helmholtz instability in a numerically simulated thunderstorm outflow.
Bulletin of the American Meteorological Society., 67 (1986), pp. 416-417
[3]
D. Terez, O. Knio.
Numerical study of the collapse of an axisymmetric mixed region in a pycnoclyne.
Physics of Fluids., 10 (1998), pp. 1438-1448
[4]
D. Terez, O. Knio.
Numerical simulation of large-amplitude internal solitary waves.
J Fluid Mechanics., 362 (1998), pp. 53-82
[5]
C. Härtel, L. Kleiser, M. Michaud, C. Stein.
A direct numerical simulation approach to the study of intrusion fronts.
J Engineering Mathematics., 32 (1997), pp. 103-120
[6]
C. Härtel, E. Meiburg, F. Necker.
Vorticity dynamics during the start-up phase of gravity currents.
Il Nuovo Cimento., 22 (1999), pp. 823-833
[7]
C. Härtel, E. Meiburg, F. Necker.
Analysis and direct numerical simulation of the flow at a gravity-current head. Part 1. Flow topology and front speed for slip and no-slip boundaries.
J Fluid Mechanics., 418 (2000), pp. 189-212
[8]
C. Härtel, F. Carlsson, M. Thunblom.
Analysis and direct numerical simulation of the flow at a gravity-current head. Part 2. The lobe-and-cleft instability.
J Fluid Mechanics., 418 (2000), pp. 213-229
[9]
F. Necker, C. Härtel, L. Kleiser, E. Meiburg.
High-resolution simulations of particle-driven gravity currents.
International J Multiphase Flow., 28 (2002), pp. 279-300
[10]
T. Tokyay, G. Constantinescu, E. Meiburg.
Lock-exchange gravity currents with a high volume of release propagating over a periodic array of obstacles.
J Fluid Mechanics., 672 (2011), pp. 570-605
[11]
L. Espatha, L. Pintoa, S. Laizetb, J. Silvestrini.
Two- and three-dimensional direct numerical simulation of particle-laden gravity currents.
Computers & Geosciences., 63 (2014), pp. 9-16
[12]
M. Cantero, J.R. Lee, S. Balachandar, M. García.
On the front velocity of gravity currents.
J Fluid Mechanics., 586 (2007), pp. 1-39
[13]
M. Cantero, S. Balachandar, M. García, J. Ferry.
Direct numerical simulations of planar and cylindrical density currents.
J Applied Mechanics., 73 (2006), pp. 923-930
[14]
M. Cantero, S. Balachandar, A. Cantelli, C. Pirmez, G. Parker.
Turbidity current with a roof: direct numerical simulation of self-stratified turbulent channel flow driven by suspended sediment.
J Geoph Res-Oceans., 114 (2009), pp. C03008
[15]
M. Cantero, S. Balachandar, G. Parker.
Direct numerical simulation of stratification effects in a sediment-laden turbulent channel flow.
J Turbulence., 10 (2009), pp. 1-28
[16]
M. Cantero, S. Balachandar, M. García.
Highly resolved simulations of cylindrical density currents.
J Fluid Mechanics., 590 (2007), pp. 437-469
[17]
M. Cantero, S. Balachandar, M. García, D. Bock.
Turbulent structures in planar gravity currents and their influence of the flow dynamics.
J Geoph Res Oceans., 113 (2008), pp. C08018
[18]
C. Cenedese, C. Adduce.
Mixing in a density-driven current flowing down a slope in a rotating fluid.
J Fluid Mechanics., 604 (2008), pp. 369-388
[19]
B. Marino, L. Thomas, P. Linden.
The front condition for gravity currents.
J Fluid Mechanics., 536 (2005), pp. 49-78
[20]
T. Fannelop, G. Waldman.
The dynamics of oil slicks: or ‘creeping crude’.
AIAAJ., 41 (1971), pp. 1-10
[21]
D. Hoult.
Oil spreading in the sea.
Ann Rev Fluid Mechanics., 4 (1972), pp. 341-368
[22]
T. Benjamin.
Gravity currents and related phenomena.
J Fluid Mechanics., 31 (1968), pp. 209-248
[23]
J. Shin, S. Dalziel, P. Linden.
Gravity currents produced by lock exchange.
J Fluid Mechanics., 521 (2004), pp. 1-34
[24]
H. Huppert, J. Simpson.
The slumping of gravity currents.
J Fluid Mechanics., 99 (1980), pp. 785-799
[25]
J. Fay.
The spreads of oil slicks on a calm sea.
Oils in the sea., pp. 53-63
[26]
J. Rottman, J. Simpson.
Gravity currents produced by instantaneous releases of a heavy fluid in a rectangular channel.
J Fluid Mechanics., 135 (1983), pp. 95-110
[27]
H. Huppert.
The propagation of two-dimensional and axisymmetric viscous gravity currents over a rigid horizontal surface.
J Fliud Mechanics., 121 (1982), pp. 43-58
[28]
M. Hallworth, H. Huppert, M. Ungarish.
Axisymmetric gravity currents in a rotating system: experimental and numerical investigations.
J Fluid Mechanics., 447 (2001), pp. 1-29
[29]
M.O. Deville, P.F. Fischer, E.H. Mund.
High-order methods for incompressible fluid flow.
Cambridge University Press, (2002),
[30]
C. Canuto, M. Hussaini, A. Quarteroni, T. Zang.
Spectral methods: fundamentals in single domains. 563.
Springer-Verlag, (2006),
[31]
C. Canuto, M. Hussaini, A. Quarteroni, T. Zang.
Spectral methods in fluid dynamics. 557.
Springer-Verlag, (1988),
[32]
Salinas J. Modelado y simulación de corrientes de gravedad con efectos de rotación. Tesis de Maestría, Instituto Balseiro, CNEA-UNC, Bariloche, Río Negro, Argentina; 2014.
[33]
M. Ungarish, H. Huppert.
The effects of rotation on axisymmetric gravity currents.
J Fluid Mechanics., 362 (1998), pp. 17-51
Copyright © 2014. IAHR y WCCE
Opciones de artículo
Herramientas
es en pt
Política de cookies Cookies policy Política de cookies
Utilizamos cookies propias y de terceros para mejorar nuestros servicios y mostrarle publicidad relacionada con sus preferencias mediante el análisis de sus hábitos de navegación. Si continua navegando, consideramos que acepta su uso. Puede cambiar la configuración u obtener más información aquí. To improve our services and products, we use "cookies" (own or third parties authorized) to show advertising related to client preferences through the analyses of navigation customer behavior. Continuing navigation will be considered as acceptance of this use. You can change the settings or obtain more information by clicking here. Utilizamos cookies próprios e de terceiros para melhorar nossos serviços e mostrar publicidade relacionada às suas preferências, analisando seus hábitos de navegação. Se continuar a navegar, consideramos que aceita o seu uso. Você pode alterar a configuração ou obter mais informações aqui.