3.1 FUNDAMENTO MATEMÁTICO

Las ecuaciones gobernantes de la dinámica de fluidos para flujo compresible, son la ecuación (1) de la conservación de la masa, ecuación (2) de la cantidad de movimiento, ecuación (3) de la conservación de la energía y la ecuación (4) de estado. En su forma compacta, las ecuaciones (1), (2) y (3) en condiciones estacionarias, se expresan como:

Siendo: 𝜌 la densidad; 𝑢 la velocidad; 𝑝 la presión; 𝜏 el tensor de tensiones, 𝐸 la energía total; 𝑇 la temperatura; 𝑅 la constante del gas. Además, 𝑘 𝑒𝑓𝑓 la conductividad térmica efectiva y 𝜏 𝑒𝑓𝑓 el tensor de tensiones efectivo. De la ecuación (4), la densidad se puede expresar como 𝜌=𝑝/𝑅𝑇. Cabe señalar, la ecuación (2) de la cantidad de movimiento está cerrada, ya que incluye el término de tensiones de Reynolds, la cual se expresa como −𝜌 𝑢 𝑖 ′ 𝑢 𝑗 ′ , donde 𝑢 ′ son las fluctuaciones. Las ecuaciones (1) y (2), en su conjunto, son conocidos como la ecuación de Navier-Stokes de número de Reynolds promedio (RANS, por sus siglas en inglés).

El parámetro dominante es el número de Mach (M). Las consideraciones del número de Mach son las siguientes: si 𝑀<0.3 el flujo se considera incompresible, si es 𝑀>0.3 el flujo es compresible. Y se clasifica como: flujo subsónico 0.3<𝑀<0.8; flujo transónico 0.8<𝑀<1.2; flujo sónico 𝑀=1; flujo supersónico 1,2<𝑀<5 ⁽⁵⁾, ⁽²⁴⁾; acotando que, White ⁽²⁴⁾ considera hasta Mach 3 como flujo supersónico. Para un valor mayor de 𝑀>5, el flujo es hipersónico, el cual es analizado a altas temperaturas ⁽²⁵⁾.

La ecuación (5) de relaciones de presiones y la ecuación (6) de relaciones de temperaturas, están en función del número de Mach, así como, de la relación de calor específico 𝛾; siendo 𝑝 0 la presión total y 𝑇 0 la temperatura total

Para el flujo de aire, para la viscosidad en función de la temperatura, se considera la ley de Sutherland (6), la cual se expresa como

Donde, 𝜇 0 =1.716 kg/(m.s) es la viscosidad de referencia, 𝑇 0 =273.11 K es la temperatura de referencia, y 𝑆=110.56 K la temperatura efectiva. La ley de Sutherland, es el resultado de una aproximación de la teoría cinética de los gases, de la idealización del potencial de la fuerza intermolecular. Donde, la temperatura efectiva 𝑆 es conocida como la constante de Sutherland el cual es una característica del gas, y están tabulados para diferentes rangos de temperatura para diferentes tipos de gases ⁽⁶⁾, ⁽²⁶⁾.

En la literatura se reportan diferentes modelos de turbulencia ⁽¹⁰⁾, y para la simulación de la turbulencia del flujo compresible se emplea el modelo de turbulencia SST 𝑘−𝜔 de Menter ⁽¹²⁾, el cual es ensamblado en RANS. El modelo de turbulencia SST 𝑘−𝜔 ha sido evaluado para diferentes condiciones del régimen del flujo compresible ⁽²⁷⁾, ⁽²⁸⁾, ⁽²⁹⁾, y es un modelo adecuado para simular la turbulencia del flujo en el presente trabajo.

El modelo de turbulencia SST 𝑘−𝜔 contiene dos ecuaciones, una para la energía cinética específica 𝑘, y la otra para la tasa de disipación específica 𝜔, por lo cual, logra mejorar las respuestas en presencia de gradientes de presión adversos, y separación del flujo, permitiendo el transporte de esfuerzos cortante en las regiones adyacentes a las paredes. Las primeras investigaciones de la turbulencia fueron desarrolladas por Kolmogorov (1941) en base a los resultados de Reynolds en 1883.

3.2 DOMINIO COMPUTACIONAL, MALLADO Y CONDICIONES DE BORDE

El dominio computacional 2D con simetría axial en el eje x para simular el flujo de aire, se muestra en la Fig. 2. No obstante, para resultados con mayor precisión numérica, se debe considerar un dominio 3D. Sin embargo, la consideración para simular el flujo en un dominio 2D es por el ahorro de tiempo de iteración y costo computacional, además, permite obtener las formas de las ondas de choque para ser analizados en la proyección del plano 2D.

El dominio computacional 2D con simetría axial comprende de una campana y la brida de admisión. Cabe señalar, las bridas de admisión poseen una campana incorporada para direccionar el flujo, y en la entrada de la misma se encuentra un filtro para evitar el ingreso de impurezas al cilindro del motor durante la etapa de admisión del flujo de aire que lo toma de la atmósfera, y dicha sección de la campana se considera en la simulación.

El mallado del dominio 2D con simetría axial se muestra en la Fig. 3, siendo una malla no estructurada con 22755 celdas cuadriláteras. La malla se refinó hacia las paredes de la campana y de la brida de admisión debido a la presencia del esfuerzo cortante en esas regiones. Antes de obtener la malla final, se refinó cuatro veces, hasta obtener una convergencia numérica, para una referencia de control en la sección de la garganta. La calidad de la malla se determinó mediante el sesgo equiángulo, el cual está definido en el rango de 0≤𝑄 𝐸𝐴𝑆 ≤1; se obtuvo un 𝑄 𝐸𝐴𝑆 =0.36, siendo la malla de buena calidad. El mallado se realizó en Workbench, en la plataforma ANSYS-Meshing.

En la misma Fig. 2, en el dominio, se señalan donde se aplican las condiciones de borde:

Para el flujo de aire a la entrada (presión y temperatura de la atmósfera), se establecen como: presión 101.325 𝑘𝑃𝑎 y temperatura 300 𝐾.

Para el flujo de aire a la salida de la brida (presión y temperatura), se establecen para seis caídas de presión: 80 𝑘𝑃??; 79 𝑘𝑃𝑎; 70 𝑘𝑃𝑎; 60 𝑘𝑃𝑎; 50 𝑘𝑃𝑎 y 40 𝑘𝑃𝑎. Siendo las relaciones de presión: 𝑝/ 𝑝 0 =𝑟𝑝 =0.7895; 0.7796; 0.6908; 0.5921; 0.4934 y 0.3947. Y la temperatura a 300 𝐾.

En las paredes de la campana y de la brida de admisión las velocidades del fluido son nulas, por la condición de no deslizamiento. Las paredes se consideran adiabáticas.

En la simetría axial, en el eje x, la velocidad del flujo en la dirección radial es nula.

3.4 MÉTODO DE SOLUCIÓN COMPUTACIONAL

Se empleó el código ANSYS-Fluent versión 12.1 que aplica el método de volumen finito (MVF), para obtener los resultados de las simulaciones numéricas. Para el método de solución se optó el análisis basado en densidad; y dominio 2D con simetría axial. Para la turbulencia se seleccionó el modelo de Menter, para la viscosidad la ecuación de la ley Sutherland. Para el flujo, la turbulencia de la energía cinética y para la disipación específica, se seleccionó la opción: Second Order Upwin. Para el monitor residual se determinó un valor fijo de 0.0001, tanto para continuidad, velocidad, y energía. Se realizó 20000-23000 iteraciones, para obtener los resultados finales de los campos de flujo del número de Mach, densidad, velocidad, presión y temperatura.

Para el procesamiento de datos se utilizó un equipo con las siguientes características: Laptop marca Síragon, modelo M54R, Intel Core 2 Duo, dos procesadores de 1.8 GHz y memoria RAM de 3 GB.

4.1 CAMPO DE FLUJO PARA EL NÚMERO DE MACH, DENSIDAD, VELOCIDAD, PRESIÓN ESTÁTICA Y TEMPERATURA ESTÁTICA

Los resultados numéricos de las variaciones de las magnitudes termodinámicas en el campo de flujo bidimensional, así como la evolución de la forma de la onda de choque a medida que se desplaza hacia la salida de la brida, se presentan en esta sección.

En un campo de flujo compresible, cuando se presenta una onda de choque por el cambio abrupto de la presión, los parámetros termodinámicos densidad, temperatura, así como la velocidad, varían drásticamente sus magnitudes, por lo cual, se presenta un patrón de flujo, y un gradiente en todo el campo de flujo, antes y después del choque.

Para cada relación de presión de la simulación del flujo, desde 𝑟𝑝=0.7895 hasta 𝑟𝑝=0.3947, se obtienen los campos del flujo para el número de Mach, densidad, velocidad, presión estática y temperatura estática, para un flujo sobreexpandido.

Al analizar el campo de flujo para el número de Mach, la cual se muestra en la Fig. 4, se observa cómo evoluciona la forma de la onda de choque y el inicio de la separación del flujo. Para 𝑟𝑝=0.7895, el flujo no ha alcanzado la velocidad sónica en la sección de la garganta, aunque no existe una onda de choque, se tiene una separación de flujo, y esto indica que el perfil aerodinámico que tiene una curvatura en la sección divergente de la pared de la brida, contribuye a que el flujo se desprenda de la pared, ocasionado por el desprendimiento de la capa límite, por lo cual, el flujo no puede continuar bordeando la pared.

Para 𝑟𝑝=0.7796, el flujo ha alcanzado la velocidad del sonido en la sección de la garganta, por lo cual el flujo está estrangulado, y se acelera a velocidad transónica. Se observa un choque normal definido al inicio de la sección divergente, el cual ocasiona la desaceleración de flujo aguas abajo. Al igual que en el caso anterior existe una separación del flujo. Como el flujo en la garganta está estrangulado, ya no permite mayor ingreso de flujo másico, por más que se disminuya la presión en la salida de la brida.

Para 𝑟𝑝=0.6908 y 𝑟𝑝=0.5921, el flujo se acelera a velocidad supersónica y se presenta la onda de choque, y la onda tiende a curvarse, y con mayor ángulo hacia la pared, y en el medio del flujo con tendencia a ser normal por la condición de simetría.

Para 𝑟𝑝=0.4934 y 𝑟𝑝=0.3947, a medida que la presión disminuye en la salida de la brida, el flujo se sigue acelerando a velocidad supersónica, la onda de choque presenta una configuración, se presenta el choque oblicuo, el reflejado y el disco de Mach. Por lo cual, por más que se disminuya la presión en la salida de la brida, el flujo se desacelera.

En la sección de la campana, en el extremo superior derecho, en la esquina de la pared, el flujo no sufre cambios significativo en su magnitud de número de Mach. A la entrada de la brida, la cual está bordeada por la campana, la velocidad del flujo se incrementa y se presenta un gradiente de velocidad subsónica, y cuya magnitud de velocidad se incrementa a medida que el flujo se aproxima a la sección de la garganta de la brida, alcanzando el régimen transónico.

El campo de flujo para la densidad se muestra en la Fig. 5, desde 𝑟𝑝=0.7895 hasta 𝑟𝑝=0.3947, las geometrías de las ondas son similares comparando con el campo de flujo para el número de Mach. Las formas de las ondas de choque para el campo de flujo para la densidad muestran las variaciones de sus magnitudes en diferentes regiones del dominio, así como donde se produce la separación de flujo, siendo la densidad de menor magnitud antes del choque.

Para el caso de la velocidad, las regiones en el campo de flujo muestran en la Fig. 6, desde 𝑟𝑝=0.7895 hasta 𝑟𝑝=0.3947, se muestra de cómo varían sus magnitudes antes y después del choque, así como, en las regiones donde se presenta la separación de flujo y en la esquina de la pared de la campana, siendo de mayor magnitud antes del choque.

En la Fig. 7 se muestra el campo de flujo para la presión estática, desde 𝑟𝑝=0.7895 hasta 𝑟𝑝=0.3947. Las distribuciones de presiones se muestran antes y después que se presenta el choque. Antes del choque la magnitud de la presión es menor, por lo cual, producto de la onda de presión abrupta, el flujo se desacelera aguas abajo después de la onda. Para 𝑟𝑝=0.7895, 𝑟𝑝=0.7796 y 𝑟𝑝=0.6908, el comportamiento del gradiente de la presión estática después de la onda tiende a ser uniforme en el centro y hacia las paredes.

Así mismo, se tiene para el caso del campo de flujo para la temperatura estática, la cual se muestra en la Fig. 8, desde 𝑟𝑝=0.7895 hasta 𝑟𝑝=0.3947. Las magnitudes de la temperatura estática sufren cambios en diferentes regiones del dominio, y la menor magnitud se presenta antes del choque, y esta caída de la temperatura está relacionada con la caída de presión.

En el campo de flujo sobreexpandido, la separación de flujo causa cargas laterales, donde la magnitud de las cargas laterales está en función de la relación de presión. Para un caso real, se presenta una distribución asimétrica del flujo debido a su naturaleza inestable, y esto se presentan para dispositivos de flujo frío, en esta caso la brida de admisión que toma el aire a la temperatura del ambiente; así como, para flujos calientes, para toberas supersónicas de cohetes, que recibe el flujo que está a alta temperatura en la cámara de combustión. Sin embargo, para una simulación 2D del campo de flujo de número de Mach, se presenta un patrón de choque, donde las formas de las ondas de choque son simétricas en el plano de simetría con respecto al eje x. No obstante, la simplificación para el análisis de los resultados en un dominio computacional 2D para el flujo sin rotación contribuye a entender la naturaleza del comportamiento del flujo bajo diferentes cargas de presión, por lo cual, tiene un aporte significativo que permite hacer estimaciones de predicciones de las formas y posiciones de las ondas de choque.

4.2 PERFILES DE NÚMERO DE MACH, DENSIDAD, VELOCIDAD, PRESIÓN ESTÁTICA Y TEMPERATURA ESTÁTICA

En esta sección se presentan los perfiles de número de Mach en la Fig. 9, de la densidad en la Fig. 10, de la velocidad en la Fig. 11, de la presión estática en la Fig. 12 y de la temperatura estática en la Fig. 13, todos obtenidos en la simetría axial del eje x, las cuales corresponde para cada campo de flujo presentado en las Fig. 4, fig. 5, fig. 6, fig. 7, fig. 8, respectivamente.

Los perfiles aportan mayor información del comportamiento del patrón del flujo, cuando disminuyen e incrementan sus magnitudes termodinámicas consecuencia de la onda de choque. Donde se presenta el choque, se aprecia una línea recta casi en vertical.

En la Fig. 9, para 𝑟𝑝= 0.7895 el flujo no ha alcanzado la velocidad sónica, por lo cual no está estrangulado. Sin embargo, está en el rango de velocidad transónica. Para las otras relaciones de presión 𝑟𝑝 = 0.7796, 0.7698, 0.6908, 0.5921 y 0.4934, el flujo está estrangulado. Para el flujo con 𝑟𝑝=0.7698, la onda de choque se encuentra en el régimen transónico, en el rango estimado de 0.9−1.1 𝑀𝑎𝑐ℎ. Para 𝑟𝑝=0.6908, el choque se encuentra en el rango de 0.8−1.3 𝑀𝑎𝑐ℎ, donde el flujo alcanzó una velocidad supersónica antes del choque, y se desaceleró a velocidad transónica. Para 𝑟𝑝=0.5921, la onda está en el rango de 0.6−1.4 𝑀𝑎𝑐ℎ, a pesar que el flujo alcanzó un valor supersónico, la caída de número de Mach corresponde para flujo subsónico. Para 𝑟𝑝=0.4934, el choque se encuentra en el rango de 0.6−1.6 𝑀𝑎𝑐ℎ, también el flujo alcanzó una velocidad supersónica, y se desaceleró a velocidad subsónica. Para ??𝑝=0.3947, el choque se produce en el rango de 0.7−1.8 𝑀𝑎𝑐ℎ, donde se evidencia que el flujo se desacelera a velocidad subsónica. Aunque el flujo se desacelera drásticamente producto del choque, los perfiles muestran que el flujo tiende a acelerarse con debilidad tratando de alcanzar la velocidad a 𝑀𝑎𝑐ℎ=1, para 𝑟𝑝=0.5921, 0.4934 y 𝑟𝑝=0.3947; excepto para las otras relaciones de presión, donde el flujo se desacelera. Sin embargo, ya no es posible que recupere su velocidad a lo que tenía antes de choque, y se debe a la influencia que tiene la geometría con curvatura del perfil aerodinámico de la sección divergente de la brida de admisión con forma de trompeta, lo cual hace que haya un desprendimiento de la capa límite.

Los perfiles de densidad se muestran en la Fig. 10, y para cada perfil se puede correlacionar con los perfiles de número de Mach mostrados en la Fig. 9, mostrando así las disminuciones e incremento de las magnitudes donde se presenta el choque.

Para los perfiles de la velocidad del flujo que se muestran en la Fig. 11, el número de Mach está relacionado, por lo cual, se muestran a que magnitud de la velocidad se alcanza antes y después del choque. Siendo el comportamiento de los perfiles similares a los del número de Mach. Donde, en la sección de la garganta, la velocidad del flujo está en el rango estimado de 280−330 𝑚/𝑠. Para la menor caída de presión, 𝑟𝑝 = 0.3947, donde se presenta el choque el flujo alcanza una velocidad máxima y mínima en el rango estimado de 480−235 𝑚/𝑠. Mientras que, para 𝑟𝑝 = 0.7895 y 𝑟𝑝 = 0.7796 el flujo alcanza una velocidad ligeramente mayor de 200 𝑚/𝑠.

En la Fig. 12, correspondiente a los perfiles de presión estática, se muestran las trayectorias de las caídas de presión hasta donde se produce el choque, y cómo se incrementa abruptamente al finalizar el choque. Después del choque, sigue disminuyendo la presión tratando de acelerar el flujo, excepto para el perfil de 𝑟𝑝 = 0.6908 que tiene un ligero incremento de la presión, lo cual indica que el flujo se sigue desacelerando, así como, para 𝑟𝑝 = 0.7895 y 𝑟𝑝 = 0.7796, respectivamente.

En la sección de la garganta de la brida, las temperaturas estática del flujo están en el rango estimado de 245 𝐾 (−28° 𝐶 ) a 255 𝐾 (−13° 𝐶), tal como se muestran en la Fig. 13. Para 𝑟𝑝 = 0.7895 y 𝑟𝑝 = 0.7796, la temperatura a la salida de la brida tiene un valor estimado de 280 𝐾 (7° 𝐶). Para 𝑟𝑝=0.6908, antes y después del choque, la temperatura se encuentra en el rango de 225−265 𝐾; para 𝑟𝑝=0.5921, en el rango de 215−280 𝐾; para 𝑟𝑝=0.4934, en el rango de 200−280 𝐾; para 𝑟𝑝=0.3947, en el rango de 180−275 𝐾, respectivamente. Donde, la menor caída de temperatura estimada de 180𝐾 (−93° 𝐶) se presenta para 𝑟𝑝=0.3947 antes que se presente el choque, y cuyo valor de la presión en la salida de la brida es de 40 𝑘𝑃𝑎. Lo cual equivale a 2.533 veces menor de la presión del flujo (101.325 𝑘𝑃𝑎) que ingresa a la brida.

Los perfiles de velocidad evaluados en la sección transversal de la divergente para diferentes distancias locales, para cada relación de presión estudiada 𝑟𝑝=0.7895, 0.7796, 0.7698, 0.6908, 0.5921 y 0.4934, se muestran en las Fig. 14, fig. 15, fig.16, fig. 17, fig. 18, fig. 19. Cuando la presión disminuye en la salida de la brida de admisión, la onda de choque se mueve hacia la salida, y en los extremos de la onda de choque se produce la separación del flujo y se presentan pequeñas regiones de recirculación del flujo.

En las Fig. 14 y fig. 15, para la distancia local 𝑥=95 𝑚𝑚, se muestra que la velocidad en la dirección radial tiene una magnitud negativa, lo cual indica que una pequeña región del flujo tiene recirculación, por lo tanto, el flujo que viene de aguas arriba se separa de la pared en la sección divergente. Donde, en esa región adyacente a la pared donde el flujo tiene velocidad negativa se presenta un gradiente de presión adverso.

A medida que la presión de salida disminuye, ver Fig. 16-18, para las distancias locales 𝑥=85 𝑦 95 𝑚𝑚, la velocidad negativa del flujo se incrementa, por lo tanto, al incrementarse la recirculación del flujo, el flujo se separa más de la pared en la dirección radial, y el gradiente de presión adverso incrementa su radio de acción.

Para la distancia local 𝑥=95 𝑚𝑚 mostrada en la Fig. 19, se cuantifica la magnitud negativa de la velocidad del flujo; mientras que, para las otras distancias locales, no hay velocidades negativas. Donde, la intensidad del gradiente de presión adverso es mayor con respecto a los casos anteriores.

A partir de la distancia local ??=95 𝑚𝑚, se puede considerar una región crítica en el desarrollo del régimen del flujo para un caso del régimen del flujo bidimensional, ya que se ha considerado que el flujo sale de la brida sin rotación. Sin embargo, se debe verificar estos resultados no concluyentes realizando una simulación en tres dimensiones en un trabajo a futuro, e incluyendo la rotación del flujo en la salida de la brida inducido por las revoluciones del turbocompresor, por lo cual se intuye que presentaría una solución más ajustada a la realidad.

De acuerdo a la consideración para la simulación numérica bidimensional del flujo compresible, se ha logrado determinar las variaciones de las magnitudes de los parámetros termodinámicos, hasta la caída de presión de 40 𝑘𝑃𝑎 en la salida de la brida. Para los casos del régimen del flujo abordado, los resultados del campo de flujo muestran que el flujo esta sobreexpandido.

CONCLUSIONES

Los resultados obtenidos de la simulación del flujo bidimensional, muestran que el flujo está sobreexpandido, y para estas condiciones se concluye que:

La curvatura de la pared de la sección divergente contribuye en el desprendimiento de la capa límite, causando transición en el desarrollo del flujo adyacente a la pared. Por lo cual, la sección divergente de la brida de admisión actúa como un difusor.

Durante las caídas de presión en la salida de la brida, en la sección divergente, la onda de choque normal fue cambiando su configuración hasta tener una configuración de forma, compuesto por un choque oblicuo, reflejado y choque normal.

El flujo en la sección divergente alcanzó una velocidad supersónica Mach 1.8 antes del choque para 𝑟𝑝=0.3947, el cual es el caso para una caída de presión de 40 𝑘𝑃𝑎 en la salida de la brida, sin embargo, se determinó que la velocidad del flujo a la salida de la brida está por debajo de Mach 1.

Los campos de flujo para la densidad, velocidad y presión, mostraron cómo varían sus magnitudes en diferentes regiones del régimen del flujo, antes y después de la onda de choque, así como, donde se produce la separación de flujo, y se determinó que existe una región crítica en la distancia de referencia 𝑥=95 𝑚𝑚, donde se produce una recirculación del flujo, la cual está cerca a la salida de la brida.

En la sección de la garganta, se determinó que la temperatura está por debajo de 255 𝐾 (−13° 𝐶), mientras que, la caída mínima de la temperatura de 180𝐾 (−93° 𝐶) se presentó para 𝑟𝑝=0.3947, esto es, para la presión de 40 𝑘𝑃𝑎 a la salida de la brida.