Repositorio Institucional de Documentos

Análisis de modelos sencillos para los términos difusivos del sistema conjunto gradiente de velocidad - gradiente de un escalar en mezcla escalar turbulenta mediante simulaciones Montecarlo Trabajo fin de Máster presentado en la Universidad de Zaragoza para la obtención del grado de Máster en Mecánica Aplicada por Luis Hernando Cifuentes Rubio Director: Dr. Jesús Joaquín Martín Yagüe POP en Ingeniería Mecánica y de Materiales Curso Académico 2010 - 2011 Centro Politécnico Superior Universidad de Zaragoza Noviembre, 2010 Análisis de modelos sencillos para los términos difusivos del sistema conjunto gradiente de velocidad - gradiente de un escalar en mezcla escalar turbulenta mediante simulaciones Montecarlo Resumen El presente trabajo comienza con un capítulo introductorio donde se expone brevemente el contexto que motiva la investigación, los antecedentes en turbulencia, los objetivos y estructura del trabajo. Posteriormente se describe la dinámica del flujo turbulento no reactivo, formulándose las ecuaciones de transporte del gradiente de velocidad en turbulencia homogénea e incompresible. Se detallan las condiciones de homogeneidad e isotropía que se suponen en el estudio y se analiza el significado físico del tensor gradiente de velocidad y el papel que desempeña en la dinámica de líneas, superficies, volúmenes fluidos y evolución en un campo escalar. Se plantea el problema de la evolución del sistema conjunto gradiente de velocidad - gradiente del escalar en la mezcla escalar turbulenta y se afronta éste por medio del método de la función densidad de probabilidad. El cierre del problema planteado se lleva a cabo con un modelo isótropo para la presión, una estimación lineal (LMSE) modificada para el término difusivo del gradiente de velocidad, la cual contiene una parte aleatoria para darle variabilidad al término y otra parte no lineal en función del gradiente de velocidad inicial e instantáneo que asegura la convergencia del sistema, y una estimación lineal para el término difusivo del gradiente del escalar. Se analiza la solución numérica de los métodos para la evolución del sistema conjunto gradiente de velocidad - gradiente del escalar, teniendo en cuenta los errores tanto globales como locales, la exactitud y estabilidad. Cada sistema se resuelve mediante un método numérico de Runge - Kutta de segundo orden. Se implementa una modelización Montecarlo para la evolución de una distribución de sistemas de gradientes, la cual reproduce numerosas características de la estadística de las pequeñas escalas, como la intermitencia en las funciones densidad de probabilidad de los tensores, así como para las magnitudes dependientes de las mismas. La estadística del campo final se analiza y se valida con resultados conocidos de LMSE, DNS y experimentales. Finalmente se presentan las conclusiones del trabajo de forma compacta en un capítulo final. Índice general 1. Introducción 1 1.1. Motivación . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 1.2. Estado del arte . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2 1.2.1. Cascada de energía . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2 1.2.2. Intermitencia . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2 1.2.3. Simulación Numérica Directa . . . . . . . . . . . . . . . . . . . . . . . . 3 1.3. Objetivos del trabajo . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4 1.4. Estructura del Trabajo . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5 2. Dinámica del flujo turbulento no reactivo 7 2.1. Descripción del comportamiento en un flujo turbulento . . . . . . . . . . . . . . 7 2.2. Tensor gradiente de velocidad, vorticidad y deformación. Invariantes . . . . . . 9 2.2.1. Tensores gradiente de velocidad, velocidad de deformación, velocidad de rotación y vector vorticidad . . . . . . . . . . . . . . . . . . . . . . . . . 9 2.2.2. Autovalores, autovectores e invariantes del los tensores gradiente de velocidad, velocidad de deformación y velocidad de rotación . . . . . . . . 10 2.2.3. Influencia del tensor gradiente de velocidad en la evolución de un campo escalar..................................... 13 3. Planteamiento del problema, modelización y resultados 15 3.1. Evolución de las PDFs ............................... 15 3.2. Modelos para los términos difusivos . . . . . . . . . . . . . . . . . . . . . . . . . 17 3.3. Implementación del modelo . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20 3.4. Campo resultante . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23 4. Conclusiones 39 i i i Capítulo 0 Bibliografía 41 Capítulo 1 Introducción La predicción de la mezcla escalar turbulenta constituye una parte esencial en numerosos problemas de la ingeniería, la industria y en proyectos medioambientales. La comprensión de la estructura de los campos de velocidades y escalares en turbulencia es todavia un problema no resuelto satisfactoriamente. No se ha podido establecer la teoría que englobe todos los fenómenos que aparecen en los flujos turbulentos y permita predecir su comportamiento en un caso general. Este trabajo de investigación toma como herramienta el método de la función densidad de probabilidad [1,2] para analizar el comportamiento estadístico que presenta la mezcla escalar turbulenta. En este capítulo se expone brevemente el contexto que motiva la investigación, los antecedentes en turbulencia y finalmente, los objetivos y estructura del trabajo. 1.1. Motivación En general, las técnicas y modelos para la predicción de flujos turbulentos constituyen una parte esencial en numerosos problemas de la ingeniería, la industria y en proyectos medioambientales. Los flujos turbulentos con escalares inertes y reactivos se presentan en numerosos procesos naturales y en muy variadas tecnologías. Sin embargo, por el momento el problema está muy lejos de ser resuelto. La enorme complejidad de estos flujos no solo excluye soluciones analíticas, sino que impide un tratamiento matemático general único y exacto, por lo que las técnicas principales para su estudio son la modelización teórica, la simulación numérica computacional y la experimentación. En un flujo turbulento la dinámica del transporte del escalar es muy compleja dado que las ecuaciones que describen el comportamiento del campo de velocidades, que son las de Navier- Stokes [3,4,5,6], son un conjunto de ecuaciones en derivadas parciales no lineales. En ellas aparecen términos que especifican fenómenos de difusión, estiramientos, cortes y rotaciones. 1 2Capítulo 1 Conocer los efectos y la importancia relativa de estos términos es esencial en la modelización de mezcla escalar turbulenta. 1.2. Estado del arte 1.2.1. Cascada de energía Kolmogorov [7,8,9] introduce en 1941 con sus dos hipótesis de semejanza una teoría para las pequeñas escalas de los flujos turbulentos a altos números de Reynolds que establece una distribución universal de la estadística de diferencias de velocidades entre puntos separados por una distancia pequeña comparada con la escala integral de la turbulencia. La primera hipótesis supone que los momentos estadísticos de estas diferencias de velocidades, llamados funciones de estructura, son funciones isótropas universales dependientes únicamente de la interdistancia, de la viscosidad cinemática y de la disipación media de energía cinética por unidad de masa. La segunda propone que sí además la separación es grande comparada con las escalas disipativas de longitud, la estadística de estas distribuciones es independiente de la viscosidad. Ambas hipótesis conducen inmediatamente, por análisis dimensional a una relación funcional explícita del enésimo momento de las diferencias de velocidades con la disipación y la interdistancia. En el espacio de Fourier lo anterior implica la famosa ley de la pendiente −5/3 [10] en escala logarítmica para el espectro de energía en el subrango inercial, siempre que se cumpla la segunda hipótesis, es decir que la separación entre puntos sea grande frente a las escalas disipativas. La imagen que se deduce de las ideas de Kolmogorov para la dinámica de la turbulencia es la de transferencia de energía de las escalas grandes a las pequeñas, donde actúa la disipación, es lo que se ha llamado Cascada de energía de Onsager - Kolmogorov. 1.2.2. Intermitencia Cuando se considera estadística de ordenes más altos, los momentos normalizados (adimensionalizados dividiendo por la potencia adecuada de la varianza) de las diferencias de velocidades deberían ser, de acuerdo con las hipótesis de semejanza de Kolmogorov, constantes universales independientes de la interdistancia r, supuesta ésta en el subrango inercial. Por el contrario se observa una dependencia de esos momentos adimensionales con el inverso de la separación. Además esta dependencia aumenta con el número de Reynolds, es decir, con la viscosidad, y con el orden del momento considerado, lo que implica fuerte intermitencia en las pequeñas escalas [11,12,13]. La imagen física asociada es la de zonas de intensos gradientes de velocidad concentradas en una fracción limitada del espacio en láminas, bandas o tubos, de dimensiones tanto más pequeñas cuanto mayor es el numero de Reynolds. Esto es, las fluctuaciones de la velocidad de mayor amplitud aumentan su probabilidad, y su función de densidad de probabilidad (en adelante PDF, Probability Density Function) es mucho más abierta que una gaussiana, Introducción 3 con momentos normalizados muy diferentes a los de aquella. (Batchelor y Townsend, 1949; Corrsin, 1972; Saffman, 1968) [14,15,16]. Cuando la distancia de separación tiende a cero, la estadística de diferencias de velocidades equivale a la de gradientes de velocidad. Si la distancia aumenta, por el contrario, se obtiene comportamiento próximo al gaussiano (Vincent y Meneguzzi, 1991; Chen et al., 1993) [17,18]. Tsinober et al., 1992 [19,20]; Zocchi et al., 1994 [21]; Maurer et al., 1994 [22]; Tabeling et al., 1996 [23] y Belin et al, 1997 [24] presentan medidas experimentales en turbulencia generada por una rejilla, en la región externa de una capa límite y en un flujo entre dos discos a contra - rotación. Confirman las relaciones de isotropía entre las componentes del gradiente de velocidad predichas por Taylor (1938) [25], el momento de tercer orden (asimetría) negativo para la derivada longitudinal y la intermitencia de la disipación. Belin et al, obtienen distribuciones del gradiente de velocidad adimensionalizada con su desviación estándar adecuada, demuestran el aumento de la asimetría de la PDF con el aumento del número de Reynolds, Tsinober et al, obtienen distribuciones de producción de enstrofía (el cuadrado del vector vorticidad) de media positiva y alineación de la vorticidad con la dirección principal asociada al autovalor intermedio del tensor velocidad de deformación (parte simétrica del tensor gradiente de velocidad) y demuestran la alineación en dirección y sentido de la vorticidad con el vector "vortex - stretching"(estiramiento de la vorticidad). Como se indicó, la noción de intermitencia implica que existen unas zonas espacialmente muy localizadas que presentan valores anormalmente altos de las variables consideradas. Si estas regiones ocupan una fracción mínima del volumen total, la contribución de estos valores altos a los momentos de orden bajo será inapreciable, pero aumentarán considerablemente en los momentos de orden más alto. Esta explicación, comúnmente aceptada, refleja por qué la teoría de Kolmogorov reproduce bien el espectro de energía, relacionado con momentos de segundo orden, y no la forma de las funciones de densidad de probabilidad, determinada por los momentos de orden mayor para valores altos (alejados más de una desviación típica) de la variable. 1.2.3. Simulación Numérica Directa La simulación numérica directa (en adelante DNS, Direct Numerical Simulation) de las ecuaciones de Navier - Stokes es una poderosa herramienta para el estudio de la turbulencia y , en concreto, para la investigación en la estructura de las pequeñas escalas. Las tres componentes de la velocidad y las nueve de su tensor gradiente pueden conocerse simultáneamente en todos los puntos del campo para cualquier instante de tiempo, lo que proporciona información completa de la estadística y de la dinámica local. Como todas las escalas son resueltas sin aproximación puede estudiarse tanto el rango inercial del espectro como el rango disipativo donde los gradientes son más altos y se producen los fenómenos de intermitencia. Sin embargo la técnica de DNS está limitada por el momento a números de Reynolds moderados, ya que la 10 Capítulo 2 Tensor velocidad de rotación: Wij(x, t) = 1 2(aij −aji) (2.8) Respectivamente Wij tiene asociado un vector axial ω(también denominado eje del tensor) y se puede expresar en función de éste como: Wij =−1 2εijkωk(2.9) donde εijk es el tensor alternador o pseudotensor de Levi - Civitá. Aquí ωse conoce como vector vorticidad y, como se puede ver fácilmente, viene definido por ω=∇×u(2.10) o en notación de subíndices ωi=εijk ∂uk ∂xj =εijkakj (2.11) La relación inversa de 2.9 permite escribir ωi=εijkWkj (2.12) 2.2.2. Autovalores, autovectores e invariantes del los tensores gradiente de velocidad, velocidad de deformación y velocidad de rotación Para estudiar el significado físico de los tensores gradiente de velocidad, velocidad de deformación y velocidad de rotación, se hace un estudio detallado de la clasificación de todas las posibles geometrías espaciales asociadas al sistema dinámico definido por el tensor correspondiente en función de los invariantes, Chong, Perry y Cantwell (1990) [30]. Para el tensor aij los autovalores satisfacen la ecuación característica: λ3+Pλ2+Qλ +R= 0 (2.13) donde λes un valor propio o autovalor del tensor ayP, Q, R sus invariantes. Se denominan invariantes del tensor aij dado que no cambian ante una transformación de coordenadas. Identificando coeficientes en la ecuación 2.13 se obtiene: P=−(a11 +a22 +a33) = −traza(a) Q=1 2P2−traza(a2)=1 2P2−aikaki(2.14) R=−det(a) = 1 3−P3+ 3PQ −traza(a3) Dinámica del flujo turbulento no reactivo 11 primero, segundo y tercer invariante respectivamente. Si se descompone aen sus partes simétricas S, y antisimétrica W, estas expresiones de 2.14 resultan: P=−(S11 +S22 +S33) = −traza(S) Q=1 2P2−SijSji −WijWji=1 2P2−SijSij +Wij Wij  =1 2P2−|S|2+|W|2(2.15) R=1 3−P3+ 3PQ −SijSjkSki −3WijWjkSki =1 3−P3+ 3PQ −traza(S3)−3traza(W2S) Introduciendo la condición de incompresibilidad, traza (a) = traza (S) = 0, P= 0 Q=1 2|W|2−|S|2(2.16) R=−traza(a3) = −1 3aikakjaji En efecto los valores propios de avienen completamente determinados por los tres invariantes. En el espacio P, Q, R, la superficie que separa regiones con valores propios reales de aquellas con valores propios complejos viene dada por el discriminante D= 0. Para el caso en estudio (flujo incompresible - figura 2.1): D=R2 4+Q3 27 (2.17) En el caso real (D <0, velocidad de deformación >vorticidad), las trayectorias de las partículas fluidas en un entorno del punto crítico pueden presentar topologías de nodos inestables, nodos estables, puntos de sillas con trayectorias en expansión, contracción y fijos. En el caso complejo (D >0, velocidad de deformación <vorticidad), se presentaran topologías de focos, si el eje del foco esta a compresión se presentara la configuración de foco inestable y si el eje del foco esta a tracción la topología es de foco estable. 12 Capítulo 2 Fig. 2.1: Topología tridimensional en el plano Q - R (P=0). Dato de Soria et al. (1992) [34] Fig. 2.2: Esbozo para la deformación de un cubo de fluido en función de λ2. Dinámica del flujo turbulento no reactivo 13 El tercer invariante Raporta información, en el caso real, sobre el carácter del valor propio intermedio λ2del tensor aij (figura 2.2). Si ′R→0, λ2→0′,si ′R > 0, λ2>0′,y si ′R < 0, λ2<0′(supuesto λ1> λ2> λ3). Entonces Rpositivo implica dos direcciones expansivas y una compresiva que deformarán un cubo de fluido hacia una lámina. Rnegativo produce, por el contrario, dos direcciones compresivas y una expansiva y el cubo en este caso se transformará en un prisma alargado. Para el tensor Sij las expresiones del los invariantes son: QS=1 2(PS−SijSji) = −1 2|S|2=−ε 4ν(2.18) RS=−1 3−P3 S+ 3PSQS−SijSjkSki=−1 3traza(S3) PSes nulo por incompresibilidad. Para el tensor Wij: QW=1 2(−WijWji) = 1 2WijWij =1 2|W|2=ω2 4(2.19) El único invariante no nulo de Wij es QWque resulta ser la mitad de la enstrofía, (el cuadrado de la vorticidad dividido por dos, Ver Apéndice ??), y es siempre positivo. Con respecto a los invariantes de aij, se cumplen las relaciones: Q=QS+QW R=RS−1 4Sijωiωj(2.20) donde Sijωiωjes el llamado término de producción de enstrofía que aparece en la ecuación de variación temporal de la enstrofía ω2 2, (Ver Apéndice ??). 2.2.3. Influencia del tensor gradiente de velocidad en la evolución de un campo escalar Los gradientes de c, ∇c, influyen decisivamente en la evolución del campo escalar hacia un estado de homogenización espacial. Si derivamos la ecuación de conservación de un escalar 2.3 con respecto a xiy denotando ∂c ∂xi por c,i(de aquí en adelante), se obtiene: Dc,i Dt =−c,iuj,i+D∇2c,i(2.21) El primer término del lado derecho de la ecuación 2.21 sin el signo menos se puede expresar de la siguiente manera: −c,iuj,i=Sijc,j+1 2εijkωjc,k(2.22) 14 Capítulo 2 El primer término representa la contribución del tensor velocidad de deformación a la variación con el tiempo de ∇c; el segundo término expresa la variación de ∇cdebida a la velocidad angular ω 2; esto nos hace referencia a que en la dinámica del transporte del escalar existen deformaciones de estiramiento, cortes y rotaciones debidos a los tensores velocidad de deformación y rotación. Si se multiplica la ecuación 2.21 por c,iy se contraen índices, resulta la ecuación para el módulo al cuadrado del gradiente del escalar, esta variable es la llamada disipación del escalar y se denota como εc. La ecuación obtenida es: ∂ ∂t(c,ic,i) + uj ∂ ∂xj (c,ic,i) = −2uj,ic,ic,j+D∇2(c,ic,i)−2D∇c,i·∇c,i(2.23) donde se ha utilizado la identidad c,i∇2c,i=∇2(c,ic,i)−2∇c,i∇c,i promediando la ecuación 2.23 se obtiene la ecuación de variación de la disipación (ó módulo del gradiente al cuadrado) del escalar, a saber: d dthc,ic,ii=−2huj,ic,ic,ji−2Dh∇c,i·∇c,ii(2.24) donde por homogeneidad se han anulado las derivadas espaciales de promedios. El segundo término de 2.24 es la disipación de la disipación del escalar, siempre negativo ya que el promedio es sobre un cuadrado. El primer término, que procede del estiramiento y rotación del gradiente del escalar, es el responsable de la amplificación o reducción del módulo del gradiente de c,i. Actúa similarmente al término de producción Sijωiωjen la ecuación de la vorticidad, ya que ui,jc,ic,j=Sij c,ic,j. Capítulo 3 Planteamiento del problema, modelización y resultados En este capítulo se plantea el problema en la evolución del flujo turbulento, los términos no lineales, locales, cerrados y no cerrados que se presentan en las ecuaciones del gradiente de velocidad y el gradiente del escalar, al igual que la modelización y el método para sus evoluciones en el tiempo. 3.1. Evolución de las PDFs El acoplamiento explicito entre las ecuaciones del escalar y su gradiente, 2.3 y2.21, es a través del término difusivo, para campos escalares inertes, en el caso en estudio la evolución de ∇cdada por la ecuación 2.21 no depende explícitamente de cy ésta se puede tratar independientemente de la ecuación 2.3. Combinando en este caso las ecuaciones de evolución para c,iy para aij, los términos que representan los efectos de las velocidades de deformación y de rotación sobre los gradientes escalares, probablemente los definitorios y caracterizantes del efecto decisivo de la turbulencia sobre el campo escalar, se reducen a términos no lineales, locales y cerrados, similares a los términos fuente en reacciones químicas autocatalíticas. Parece apropiado atacar este problema en el marco de los métodos de PDF o, equivalentemente, de los métodos estocásticos. Se desea ahora obtener la ecuación de transporte de la PDF del tensor gradiente de velocidad y gradiente del escalar. Sea aij(x, t) el valor de ∂ui ∂xj en el punto xen tiempo tyAij su correspondiente variable asociada en el espacio de probabilidad. Siguiendo a Pope (1985) [2] y Dopazo (1994) [1], la evolución de la PDF conjunta del gradiente del escalar y del gradiente de velocidad Pa,∇cen los espacios de probabilidades Aij para el 15 16 Capítulo 3 gradiente de velocidad y C,ipara el gradiente del escalar, se tiene: DPa,∇c Dt +∂ ∂Aij Daij Dt a=A, ∇c=∇CPa,∇c +∂ ∂C,iDc,i Dt a=A, ∇c=∇CPa,∇c= 0 (3.1) Para cerrar esta ecuación de transporte los promedios condicionados deben ser expresados en términos conocidos. El primer término condicionado en 3.1 es debido a la evolución del tensor gradiente de velocidad en el espacio de probabilidad: ∂ ∂Aij Daij Dt a=A, ∇c=∇CPa,∇c(3.2) Esta ecuación 3.2 es de magnitud instantánea y el gradiente del escalar no influye en ella, se presenta el desarrollo de la ecuación desde una referencia lagrangiana siguiendo la trayectoria de una partícula, la evolución de la característica en el espacio de probabilidad Acorrobora la independencia condicionada al gradiente del escalar: Daij Dt a=A=h−aikakj|a=Ai +−1 ρ ∂2p ∂xi∂xja=A +hν∇2aij|a=Ai(3.3) El primer término del lado derecho en la ecuación 3.3 es simplemente −AikAkj, que es cerrado y no necesita modelización. Los términos segundo y tercero son, respectivamente, los correspondientes a la influencia de la presión y al efecto de la difusión viscosa de la ecuación del tensor gradiente de velocidad, ambos no cerrados y contienen promedios condicionados a que el gradiente de velocidad sea Aij. El sistema de Lagrange asociado a la ecuación 3.2 es: dt 1=dxi ui(x, t)=dAij Daij Dt aij =Aij (3.4) De la primera igualdad en 3.4 resultan las trayectorias de las partículas fluidas en el sistema fluido y de la segunda pueden obtenerse las trayectorias de Aij en el espacio de probabilidad. Planteamiento del problema, modelización y resultados 17 Estas ecuaciones son de magnitud instantáneas, a saber: dAij ="Daij Dt aij =Aij #dt dAij dt =Daij Dt aij =Aij =Daij Dt a=A y de la ecuación 3.3 tenemos: dAij dt =−AikAkj +−1 ρ ∂2p ∂xi∂xja=A +hν∇2aij|a=Ai(3.5) La ecuación 3.5 produce una evolución determinista de Aij en el espacio de probabilidad o de fases, A, mientras su variable instantánea asociada aij(x, t) es estocástica. El segundo término condicionado en la ecuación 3.1 es debido a la evolución del tensor gradiente del escalar en el espacio de probabilidad y tiene la forma: ∂ ∂C,iDc,i Dt a=A, ∇c=∇CPa,∇c(3.6) Para esta ecuación del gradiente del escalar 3.6 se realiza el mismo análisis que se planteo para el gradiente de velocidad en la ecuación 3.2 y se llega a que en una referencia lagrangiana que sigue la partícula, la evolución de las características en el espacio de probabilidad, C, asociado al gradiente del escalar es: dC,i dt =−hc,jaji|a=A, ∇c=∇Ci+hD∇2c,i|a=A, ∇c=∇Ci(3.7) donde Aji es la variable asociada a aji(x, t) en el espacio de probabilidad A. En la ecuación 3.7 el primer término del lado derecho resulta −C,jAji y es cerrado, mientras que el segundo es no cerrado, este es correspondiente al efecto de la difusión viscosa del tensor gradiente del escalar. dC,i dt ≡ −C,jAji +hD∇2c,i|a=A, ∇c=∇Ci(3.8) 3.2. Modelos para los términos difusivos El término viscoso en las ecuaciones 3.5 y3.8 es formalmente idéntico al de difusión másica que aparece en problemas de mezcla de escalares. Autores como Dopazo [35,36], Girimaji [37], Pope [2] y otros, han propuestos cierres para este 18 Capítulo 3 término difusivo y han contribuido formidablemente en la interpretación de las topologías de flujos turbulentos. Dopazo [35,36] sugiere un modelo determinista para la difusión condicionada al valor del escalar C(x, t) = φexpresable como: hD∇2C(x, t)|C(x, t) = φi=−k(φ−hCi) (3.9) donde φes la variable asociada a C(x, t) en el espacio de probabilidad. En este modelo la variable φse desplaza hacia la media, hCi, con una velocidad proporcional a su distancia a ella. La constante kpuede considerarse el inverso de un tiempo característico de decaimiento de las fluctuaciones del campo escalar, que se elige de modo que ajuste la caída de la varianza del escalar. Este modelo, que se denominará [35]LMSE (Linear Mean Square Estimation), no altera ni la media ni la forma de la PDF. De lo anterior, para la ecuación 3.5 el cierre de los dos últimos términos de promedios condicionados para este trabajo fueron modelados usando una aproximación isotrópica del tensor Hessiano para la presión (J. Martín, C. Dopazo, and L. Valiño, 1998 [38])y una aproximación estimada para la difusión viscosa inspirada en el modelo LMSE (Dopazo, 1974 [35]). El resultado de la ecuación 3.5 queda de la siguiente manera: dAij dt ≡ −AikAkj +1 3δijAlmAml −ω0Aij (3.10) donde la constante ω0, que tiene dimensiones de inverso de tiempo, puede asociarse con un tiempo característico de decaimiento τa=ω−1 0, de los gradientes del campo de velocidad. Asimismo el término difusivo de la ecuación 3.8 se cerró con un modelo LMSE y su resultado es: dC,i dt ≡ −C,jAji −ω1C,i(3.11) ω1se considera como la frecuencia asociada al decaimiento viscoso del gradiente del escalar. En el caso del tensor aij, la convergencia para el sistema dependerá, sin embargo, del valor de ω0; pero, para todo valor inicial Aij (0),(i= 1,...3, j = 1,...3),existe un valor de ω0tal que Aij(t→ ∞)→0, este valor no necesita ser mayor de dos ó tres veces el módulo kAk,(se usa la definición de módulo para matrices: kAk= (Aij Aij ) 1 2), (J. Martín, 1995 [39]). En este trabajo el modelo propuesto de ω0se compone de una parte constante/aleatoria y otra parte no lineal en función del módulo inicial kA(n) ij (0)ky el módulo de la partícula montecarlo en cada instante kA(n) ij (t)k. La primera parte le da variabilidad a la relajación del sistema y la segunda le da convergencia debido a que los términos cuadráticos que existen en la ecuación 3.10 pueden dominar sobre la difusión lineal empleada en la evolución del modelo; cabe resaltar que usando el modelo LMSE lineal y constante el gradiente de velocidad diverge en tiempo Planteamiento del problema, modelización y resultados 19 finito. ω0=ωa+ωb(3.12) De la ecuación 3.10 tenemos: dAij dt =−AikAkj +1 3δijAlmAml −ω0Aij dAij dt =−AikAkj +1 3δijAlmAml −(ωa+ωb)Aij (3.13) La parte constante/aleatoria en 3.12 es llamada ωa: ωa=µa(1 + rand) (3.14) rand es una variable gaussiana de media cero y varianza uno, por tanto es la parte aleatoria de ωa, y µaes la parte constante. Las dos partes (rand y µa)cumplen con la condición de tener un valor en cada paso de tiempo menor a tres veces el módulo kAk(J. Martín, 1995 [39]). Bajo los parámetros empleados, el valor de µa= 2,55s−1. La parte que garantiza la convergencia de los gradientes Aij es ωb: ωb=µbkAij(t)k kAij(0)k2 (3.15) La frecuencia ωbesta regulada con el módulo al cuadrado de la partícula montecarlo en cada instante. De acuerdo a la ecuación 3.13, con ωbyAij en el último término, la difusión contiene un término cubico que dominará sobre el término cuadrático de estiramiento si los gradientes crecen. Por tanto esta parte garantiza la no divergencia a ∞de los gradientes Aij. Si el módulo kAij(t)krespecto a kAij(0)ken la ecuación 3.15 es muy grande, su resultado elevado al cuadrado lleva a ωba valores más altos para relajar el gradiente en el instante (ecuación 3.13)y viceversa. µbes la parte que hace que ωbsea pequeña dado que el aporte de ωbaω0 debe cumplir con la condición de tener un valor en cada paso de tiempo menor a tres veces el módulo kAk. Bajo los parámetros empleados se ha tomado µb= 0.05s−1. El último término en la ecuación 3.13 queda: ω0Aij = (ωa+ωb)Aij = (ωa11 +ωb11 )A11 (ωa21 +ωb21 )A21 (ωa31 +ωb31 )A31 (ωa12 +ωb12 )A12 (ωa22 +ωb22 )A22 (ωa32 +ωb32 )A32 (ωa13 +ωb13 )A13 (ωa23 +ωb23 )A23 (ωa33 +ωb33 )A33  (3.16) Por lo descrito, con el modelo propuesto existirán frecuencias ω0distintas en cada paso de tiempo para cada componente del tensor aij de las partículas montecarlo. Cabe señalar que para el tensor aij, en cada paso de tiempo, se ajustan las componentes de su diagonal de modo que a11 +a22 +a33 = 0, respetando de esta manera la condición de incompresibilidad, ∇·u= 0. 26 Capítulo 3 0.0001 0.001 0.01 0.1 1 -10 -5 0 5 10 PDF (a11) a11 Modelo Belin et. al., 1997 Gaussiana Fig. 3.4: PDF final para a11 la gaussianidad, y por lo visto en la figura 3.4 el modelo representa de manera cualitativa y cuantitativamente aproximada este comportamiento. 0.0001 0.001 0.01 0.1 1 -10 -5 0 5 10 PDF (a12) a12 DNS Modelo Fig. 3.5: PDF final para a12 La P(a12) de a12 es comparada en la figura 3.5 con resultados de DNS de tamaño 1283 (J. Martín, C. Dopazo, and L. Valiño) [45] con un Reλ≃47. El modelo presenta caídas exponenciales de forma simétrica que siguen el comportamiento conocido de DNS. La PDF de una componente de la vorticidad se muestra en la figura 3.6, esta resulta casi gaussiana, Planteamiento del problema, modelización y resultados 27 1e-05 0.0001 0.001 0.01 0.1 1 -15 -10 -5 0 5 10 15 PDF (ω(1)) ω(1) Modelo Gaussiana Fig. 3.6: PDF final para ω(1) contrariamente al hecho establecido de que las componentes de la vorticidad son variables intermitentes con distribuciones similares a las de las componentes no diagonales de aij, figura 3.5. EL modelo falla en la estadística para momentos altos de la vorticidad. 0.0001 0.001 0.01 0.1 1 -15 -10 -5 0 5 10 15 PDF α β γ Fig. 3.7: PDF final de DNS para los autovalores α, β yγde Sij El tensor Sij, por ser simétrico, tiene siempre tres autovalores propios reales que además, por la incompresibilidad, suman cero. Se denotarán estos autovalores por α,βyγ, con α > 0 y γ < 0. βes el intermedio, que puede tener cualquier signo, figura 2.2.Betchov (1956) [46] 28 Capítulo 3 0 0.5 1 1.5 2 2.5 3 3.5 4 0 0.2 0.4 0.6 0.8 1 PDF cos (ω | αSij) cos (ω | βSij) cos (ω | γSij) Fig. 3.8: PDF final de DNS para los cosenos de los ángulos de la vorticidad ωy los autovectores de Sij 1e-06 1e-05 0.0001 0.001 0.01 0.1 1 -15 -10 -5 0 5 10 15 PDF (β) β inicial final Fig. 3.9: PDF del modelo para el autovalor βde Sij demostró que el promedio hαβγies proporcional al skewness de a11, típicamente negativo (Ver Apéndice ??), por lo que la configuración más probable sería con β > 0. La figura 3.7 confirma este extremo. Además, presentan fuerte correlación Sij y la vorticidad ω, teniendo esta última tendencia a la alineación con la dirección principal asociada al autovalor intermedio β, figura 3.8. Los resultados del modelo para los autovalores de Sij se representan en la figura 3.10. Planteamiento del problema, modelización y resultados 29 0.001 0.01 0.1 1 -15 -10 -5 0 5 10 15 PDF α β γ Fig. 3.10: PDF final del Modelo para los autovalores α, β yγde Sij P(α), P(β),yP(γ) son inicialmente gaussianas y el modelo las lleva hacia estados finales intermitentes con aumento de la asimetría y del aplastamiento. En la figura 3.9 se observa el desplazamiento de P(β) hacia la derecha de modo que su media es positiva, hecho conocido de DNS, figura 3.7. 0 0.5 1 1.5 2 2.5 3 0 0.2 0.4 0.6 0.8 1 PDF cos (ω | αSij) cos (ω | βSij) cos (ω | γSij) Fig. 3.11: PDF final del Modelo para los cosenos de los ángulos de la vorticidad ωy los autovectores de Sij La vorticidad en el modelo resulta fuertemente alineada con el autovector iβasociado al autovalor intermedio βde Sij, según puede verse en la figura 3.11.Pcos(ω|αSij )yPcos(ω|γSij ) 30 Capítulo 3 presentan máximos entorno al cero indicando tendencia a la perpendicularidad. Estos resultados coinciden con los estudios de DNS, figura 3.8, y con los resultados experimentales de Tsinober et al. (1992) [19,20]. 1e-05 0.0001 0.001 0.01 0.1 1 10 -5 0 5 10 15 PDF (Sij ωi ωj) Sij ωi ωj Fig. 3.12: PDF final del modelo para la producción de enstrofia Sij ωiωj La producción de enstrofía (Sijωiωj- Ver Apéndice ??), responsable del crecimiento del cuadrado de la vorticidad presenta una PDF asimétrica extremadamente abierta, con caída más lenta para valores positivos, figura 3.12. Aquí cabe señalar que el modelo, a partir de un campo inicial donde únicamente se fijan las correlaciones hasta orden dos, produce hSij ωiωji(una correlación de orden tres) positiva en concordancia con lo expuesto en la ecuación 2.6, ya que el modelo genera también promedio del cubo de a11 negativo. Por lo tanto el término Sijωiωj, efectivamente produce aumento del módulo de la vorticidad, según el resultado. Esta característica del campo, anticipada por Betchov (1956) [46] aparece en todos los resultados de DNS. La distribución final del segundo invariante Qy tercer invariante Rde aij se representan en las figuras 3.13 y3.14. La P(Q) desplazada hacia valores negativos nos confirma el comportamiento dominante del estiramiento en la estadística del modelo, (ecuación 2.16). P(R) aparece desplazada hacia valores positivos, el resultado es coherente con el promedio positivo del autovalor intermedio β, a la vista de la ecuación 2.16. La ecuación 3.11 para las características del gradiente del escalar se implementó mediante el método Montecarlo descrito en la sección 3.3.C,(n) i(0) se tomó de componentes gaussianas independientes, con matriz de covarianzas unidad y sin correlación con A(n) ij (0). La PDF final de una componente del gradiente del escalar se muestra en la figura 3.15, el resultado es simétrico con caídas exponenciales, alejada de la gaussianidad, similar en su forma a la de las componentes no diagonales de Aij pero menos abierta. Su media es nula, al igual que el Planteamiento del problema, modelización y resultados 31 1e-05 0.0001 0.001 0.01 0.1 1 10 -14 -12 -10 -8 -6 -4 -2 0 2 PDF (Q) Q Fig. 3.13: PDF final del modelo para el segundo invariante de aij 1e-05 0.0001 0.001 0.01 0.1 1 10 0 2 4 6 8 10 PDF (R) R Fig. 3.14: PDF final del modelo para el tercer invariante de aij coeficiente de asimetría y su coeficiente de aplastamiento alcanza valores mayores a 3 indicando la no gaussianidad. El comportamiento del gradiente del escalar, con la modelización utilizada, resulta intermitente, como sucede con la del gradiente de velocidad. En la figura 3.16 se representa la PDF final del producto aij c,ic,jprocedente de la ecuación 2.24. Como se mencionó, este término actúa similarmente al término de producción Sijωiωj en la ecuación de la vorticidad, ya que aijc,ic,j=Sij c,ic,j. La media se hace rápidamente 32 Capítulo 3 0.0001 0.001 0.01 0.1 1 10 -15 -10 -5 0 5 10 15 PDF (c,i) c,i DNS Modelo Fig. 3.15: PDF final para c,i negativa al principio de la evolución, alcanzando un valor absoluto relativamente alto. Se aprecian en la figura 3.16 la fuerte asimetría y el alto aplastamiento de la PDF, resultando así el comportamiento intermitente. Además, el hecho de que esta variable tenga promedio negativo durante toda la evolución produce un término de producción de disipación del escalar positivo (ecuación 2.24), por lo que amplifica de manera continua el módulo del gradiente del escalar. 1e-06 1e-05 0.0001 0.001 0.01 0.1 1 10 -6 -5 -4 -3 -2 -1 0 PDF (aij c,i c,j) aij c,i c,j Fig. 3.16: PDF final del modelo para la producción del modulo del gradiente del escalar, (aij c,ic,j) Planteamiento del problema, modelización y resultados 33 0.0001 0.001 0.01 0.1 1 10 0 2 4 6 8 10 12 14 PDF (εc) εc DNS Modelo Fig. 3.17: PDF final para la disipación del escalar εc La PDF de la disipación del escalar se aprecia en la figura 3.17. Observamos que el modelo en comparación con el resultado de DNS concuerda cualitativamente, ya que P(εc) presenta caída más lenta para valores altos de la variable. 0 0.5 1 1.5 2 2.5 3 3.5 0 0.2 0.4 0.6 0.8 1 PDF cos ( c,i | αSij) cos (c,i | βSij) cos (c,i | γSij) Fig. 3.18: PDF final de DNS para los cosenos de los ángulos del gradiente del escalar c,iy los autovectores de Sij Las figuras 3.18 por DNS y 3.19 del modelo muestran las PDFs de los cosenos de los ángulos que forma el gradiente del escalar con las tres direcciones principales del tensor velocidad de deformación, correspondientes a los autovalores α, β yγ. Se observa una gran probabilidad de 34 Capítulo 3 0 1 2 3 4 5 6 0 0.2 0.4 0.6 0.8 1 PDF cos ( c,i | αSij) cos (c,i | βSij) cos (c,i | γSij) Fig. 3.19: PDF final del Modelo para los cosenos de los ángulos del gradiente del escalar c,iy los autovectores de Sij alineación de c,icon la dirección correspondiente al autovalor negativo γ. Anteriormente se vió que β, el intermedio, presenta promedio positivo, lo que producía mayor probabilidad de dos direcciones expansivas y una compresiva para la deformación en torno al punto fluido, lo que tendería a deformar un cubo a una estructura en forma de lámina (sección 2.2.2. El resultado para la alineación de c,icon la dirección compresiva indica que, en promedio, el gradiente del escalar sería perpendicular a dicha lámina. Este tipo de topología refleja el hecho de que ∇ces perpendicular a la superficie isoescalar c(x, t) = constante, favoreciendo el transporte molecular del escalar. Para el estudio de correlación entre magnitudes de interés en el campo resultante del modelo se han empleado diagramas de dispersión (’scatter - plots’) de 5000 partículas, que dan una idea, al menos cualitativa, de la correlación entre las variables. Todas las variables presentadas en las dispersiones están normalizadas (adimensionalizadas) con el segundo invariante del tensor velocidad de rotación QW. Jiménez et al., 1993 [28], definieron una magnitud para la dinámica de la vorticidad llamada Streching rate: σω≡Sijωiωj ωkωk (3.21) yMartín et al., 2001 [45], definieron otra magnitud en la dinámica del gradiente del escalar llamada Scalar gradient streching rate: σ∇c≡Sijc,ic,j c,kc,k (3.22) Planteamiento del problema, modelización y resultados 35 Estas variables, σωyσ∇cbrindan información de las magnitudes relativas de la enstrofía ω2 y la disipación del escalar εcfrente a los respectivos términos de producción. -15 -10 -5 0 5 10 15 -20 -15 -10 -5 0 5 10 15 20 Q R (a) DNS -15 -10 -5 0 5 10 15 -20 -15 -10 -5 0 5 10 15 20 Q R (b) Modelo Fig. 3.20: Dispersión Rvs Q 0 2 4 6 8 10 12 14 0 10 20 30 40 50 60 QW - QS (a) DNS 0 2 4 6 8 10 12 14 0 10 20 30 40 50 60 QW - QS (b) Modelo Fig. 3.21: Dispersión −QSvs QW El diagrama de dispersión del segundo y tercer invariantes de aij (R, Q) se muestra en la figura 3.20; vemos que el modelo llega a valores de Qpositivos entrando en la zona compleja donde la velocidad de deformación es mucho menor que la vorticidad, muestra también un llenado en la región de discriminante nulo y por debajo de ella misma donde los autovalores son reales, en la región de discriminante nulo se extiende hacia la zona de R > 0, lo que explica el promedio positivo del autovalor intermedio de Sij . Para la correlación de −QSyQW, figura 3.21, el modelo parece reproducir la geometría conocida de DNS, sin embargo se aprecia que el estiramiento es dominante ante la vorticidad debido a que no hay muchos valores altos para QWen la dispersión. En las figuras 3.22 y3.23 se muestra la correlación que −QScon las magnitudes σ∇cyσω. Los valores altos del estiramiento hacen que se amplifique la disipación del escalar y la vorticidad, 42 Capítulo 4 [14] G. K. Batchelor and A. A. Townsend. The nature of turbulent motion at large wave - numbers. Proc. R. Soc. Lond. A., 199:238–255, 1949. 3 [15] A. Kuo and S. Courrsin. Experiment on the geometry of the fine - structure regions in fully turbulent fluid. J. Fluid Mech., 56:447–479, 1972. 3 [16] P. G. Saffman. In Vortex dynamics, California institute of technology. Cambridge University, 1992. 3 [17] A. Vincent and M. Meneguzzi. The spatial structure and statistical properties of homogeneous turbulence. J. Fluid Mech., 225:1–20, 1991. 3,4 [18] R. H. Kraichnan S. Chen, G. D. Doolen and S. She. On statistical correlations between velocity increments and locally averaged dissipation in homogeneous turbulence. Phys. Fluids A., 5:458–463, 1993. 3 [19] E. Kit A. Bershadskii and A. Tsinober. On universality of geometrical invariants in turbulence - experimental results. Phys. Fluids A., 5:1523–1525, 1993. 3,30 [20] E. Kit A. Tsinober and T. Dracos. Experimental investigation of the field velocity gradients in turbulence flows. J. Fluids Mech., 242:162–192, 1992. 3,30 [21] J. Maurer G. Zocchi, P.Tabeling and H. Willaime. Measurement of the scaling of the dissipation at high reynolds numbers. Phys. Rev., 50:3693–3700, 1994. 3 [22] P. Tabeling J. Maurer and G. Zocchi. Statistics of turbulence between two counterrotating disks in low-temperature helium gas. Europhys. Lett., 26:31–36, 1994. 3 [23] F. Belin J. Maurer P. Tabeling, G. Zocchi and H. Willaime. Probability density functions, skewness, and flatness in large reynolds number turbulence. Phys. Rev. E, 53:1613–1621, 1996. 3 [24] P. Tabeling F. Belin, J. Maurer and H. Willaime. Velocity gradient distributions in fully developed turbulence: An experimental study. Phys. Fluids, 9:3843–3850, 1997. 3,25 [25] G. I. Taylor. Production and dissipation of vorticity in a turbulent fluid. Proc. R. Soc. Lond. A., 164:15–23, 1938. 3 [26] E. D. Siggia. Point-vortex simulation of the inverse energy cascade in two-dimensional turbulence. Phys. Fluids, 24:171–173, 1981. 4 [27] R. M. Kerr. Higher - order derivative correlations and the alignment of small - scale structures in isotropic numerical turbulence. J. Fluid Mech., 153:31–58, 1985. 4 [28] P. Saffman J. Jiménez, A. Wray and R. Rogallo. The structure of intense vorticity in homogeneous turbulence. J. Fluid Mech., 255:65–90, 1993. 4,34 Bibliografía 43 [29] J. Jiménez. Kinematic alignment effects in turbulent flows. Phys. Fluids A., 4:652–654, 1992. 4 [30] A. E. Perry M. S. Chong and B. J. Cantwell. A general classification of three - dimensional flow fields. Phys. Fluids A., 2:765–777, 1990. 4,10 [31] B. J. Cantwell. Exact solution of a restricted euler equation for the velocity gradient tensor. Phys. Fluids, 4:782–793, 1992. 4 [32] B. J. Cantwell. On the berhavior of velocity gradient tensor invariants in direct numerical simulations of turbulence. Phys. Fluids, 5:2008–2013, 1993. 4 [33] J. O. Hinze. In Turbulence. Mc. Graw Hill, 1975. 9,20 [34] R. Sondergaard A. E. Perry J. Soria, M. S. Chong and B. J. Cantwell. A study of the fine scale motions of incompressible time - developing mixing layers. Center for Turbulence Research Proceedings of the Summer Program, pages 101–121, 1992. 12 [35] C. Dopazo and E. E. O’Brien. An approach to the autoignition of a turbulent mixture. Acta Astronaut., 1:1239–1266, 1974. 17,18 [36] C. Dopazo. Probability density function approach for a turbulent axisymmetric heated jet. centerline evolution. Phys. Fluids, 18:397–404, 1975. 17,18 [37] S. S. Girimaji and S. B. Pope. A diffusion model for velocity gradients in turbulence. Phys. Fluids A nº2, 2:242–256, 1990. 17 [38] C. Dopazo J. Martín and L. Valiño. Dynamics of velocity gradient invariants in turbulence: Restricted euler and linear diffusion models. Phys. Fluids, 10:2012–2025, 1998. 18 [39] J. J. Martín. Modelización de la intermitencia interna del gradiente de velocidad y del gradiente de un escalar en turbulencia homogenea. PhD thesis, Universidad de Zaragoza, España, 1995. 18,19,20,24 [40] S. B. Pope. A monte carlo method for the pdf equations of turbulent reactive flow. Combustion science and technology, 25:159–174, 1981. 20 [41] S. B. Pope. Monte carlo calculations of premixed turbulent flames. Eighteenth symposium (international) on combustion, pages 1001–1010, 1981. 20 [42] L. Valiño. Computación de flujos turbulentos homogeneos con reacción química. Simulación Numérica Montecarlo de velocidades, concentraciones y gradientes de concentraciones. PhD thesis, Universidad de Zaragoza, España, 1989. 20 [43] J. P. Hierro. Métodos estadísticos para simular la evolución de campos escalares y de gradientes en turbulencia homogénea e isótropa. PhD thesis, Universidad de Zaragoza, España, 2002. 20 44 Capítulo 4 [44] S. C. Chapra and R. P. Canale. Métodos de runge - kutta. In McGraw Hill Interamericana, editor, Métodos numéricos para ingenieros, pages 713 – 759. McGraw - Hill, México, 2003. 22 [45] C. Dopazo J. Martín and L. Valiño. Averaged dynamics in the phase planes of a scalar field using dns data. Journal of Turbulence, 2001. 25,26,34,37 [46] R. Betchov. An inequality concerning the production of vorticity in isotropic turbulence. J. Fluid Mech., 1:497–504, 1956. 27,30