Método de volúmenes finitos en ecuaciones diferenciales parciales
En el estudio de ecuaciones diferenciales parciales , el esquema MUSCL es un método de volumen finito que puede proporcionar soluciones numéricas de alta precisión para un sistema dado, incluso en casos en los que las soluciones presentan choques, discontinuidades o grandes gradientes. MUSCL significa Esquema Monótono Centrado en la Corriente Arriba para Leyes de Conservación (van Leer, 1979), y el término fue introducido en un artículo seminal de Bram van Leer (van Leer, 1979). En este artículo construyó el primer esquema de variación total decreciente (TVD) de alto orden donde obtuvo precisión espacial de segundo orden.
La idea es reemplazar la aproximación constante por partes del esquema de Godunov por estados reconstruidos, derivados de estados promediados de celdas obtenidos del paso de tiempo anterior. Para cada celda, se obtienen estados izquierdo y derecho reconstruidos y limitados por pendiente, y se utilizan para calcular flujos en los límites de la celda (bordes). Estos flujos pueden, a su vez, utilizarse como entrada para un solucionador de Riemann , después de lo cual se promedian las soluciones y se utilizan para avanzar la solución en el tiempo. Alternativamente, los flujos pueden utilizarse en esquemas sin solucionador de Riemann , que son básicamente esquemas similares a Rusanov.
Reconstrucción lineal
Consideraremos los fundamentos del esquema MUSCL considerando el siguiente sistema simple, escalar, 1D, de primer orden, que se supone que tiene una onda que se propaga en la dirección positiva,
Donde representa una variable de estado y representa una variable de flujo .
El esquema básico de Godunov utiliza aproximaciones constantes por partes para cada celda y da como resultado una discretización de primer orden en contra del viento del problema anterior con centros de celdas indexados como . Un esquema semidiscreto se puede definir de la siguiente manera:
Este esquema básico no es capaz de manejar choques o discontinuidades agudas ya que tienden a difuminarse. Un ejemplo de este efecto se muestra en el diagrama opuesto, que ilustra una ecuación advectiva 1D con una onda escalonada que se propaga hacia la derecha. La simulación se llevó a cabo con una malla de 200 celdas y se utilizó un integrador temporal de Runge-Kutta de cuarto orden (RK4).
Para proporcionar una mayor resolución de las discontinuidades, el esquema de Godunov se puede ampliar para utilizar aproximaciones lineales por partes de cada celda, lo que da como resultado un esquema de diferencia central que tiene una precisión de segundo orden en el espacio. Las aproximaciones lineales por partes se obtienen a partir de
Así, evaluando los flujos en los bordes de la celda obtenemos el siguiente esquema semidiscreto
donde y son los valores aproximados por partes de las variables del borde de la celda, es decir ,
Aunque el esquema de segundo orden anterior proporciona una mayor precisión para soluciones suaves, no es un esquema de disminución de variación total (TVD) e introduce oscilaciones espurias en la solución donde hay discontinuidades o choques. Un ejemplo de este efecto se muestra en el diagrama opuesto, que ilustra una ecuación advectiva 1D , con una onda escalonada que se propaga hacia la derecha. Esta pérdida de precisión es de esperar debido al teorema de Godunov . La simulación se llevó a cabo con una malla de 200 celdas y se utilizó RK4 para la integración temporal.
Los esquemas numéricos basados en MUSCL extienden la idea de utilizar una aproximación lineal por partes para cada celda mediante el uso de estados extrapolados izquierdo y derecho con pendiente limitada . Esto da como resultado el siguiente esquema de discretización TVD de alta resolución:
Lo cual, alternativamente, puede escribirse de forma más sucinta,
Los flujos numéricos corresponden a una combinación no lineal de aproximaciones de primer y segundo orden a la función de flujo continuo.
Los símbolos y representan funciones dependientes del esquema (de las variables de borde de celda extrapoladas limitadas), es decir ,
donde, utilizando pendientes a favor del viento:
y
La función es una función limitadora que limita la pendiente de las aproximaciones por partes para garantizar que la solución sea TVD, evitando así las oscilaciones espurias que de otro modo ocurrirían alrededor de discontinuidades o choques (consulte la sección Limitador de flujo) . El limitador es igual a cero cuando y es igual a la unidad cuando . Por lo tanto, la precisión de una discretización TVD se degrada a primer orden en los extremos locales, pero tiende a segundo orden en partes suaves del dominio.
El algoritmo es muy sencillo de implementar. Una vez que se ha elegido un esquema adecuado , como el esquema de Kurganov y Tadmor (ver más abajo), se puede proceder a la solución utilizando técnicas de integración numérica estándar.
Esquema central de Kurganov y Tadmor
Un precursor del esquema central de Kurganov y Tadmor (KT) (Kurganov y Tadmor, 2000) es el esquema central escalonado de Nessyahu y Tadmor (NT) (Nessyahu y Tadmor, 1990). Es un esquema de alta resolución , de segundo orden y sin solucionador de Riemann que utiliza la reconstrucción MUSCL. Es un método completamente discreto que es fácil de implementar y se puede utilizar en problemas escalares y vectoriales , y se puede ver como un flujo de Rusanov (también llamado flujo local de Lax-Friedrichs) complementado con reconstrucciones de alto orden. El algoritmo se basa en diferencias centrales con un rendimiento comparable a los solucionadores de tipo Riemann cuando se utiliza para obtener soluciones para EDP que describen sistemas que exhiben fenómenos de alto gradiente.
El esquema KT extiende el esquema NT y tiene una cantidad menor de viscosidad numérica que el esquema NT original. También tiene la ventaja adicional de que se puede implementar como un esquema completamente discreto o semidiscreto . Aquí consideramos el esquema semidiscreto.
El cálculo se muestra a continuación:
Donde la velocidad de propagación local , , es el valor absoluto máximo del valor propio del jacobiano de las sobre celdas dado por
Más allá de estas velocidades relacionadas con CFL , no se requiere ninguna información característica.
El cálculo de flujo anterior se denomina con mayor frecuencia flujo de Lax-Friedrichs (aunque vale la pena mencionar que dicha expresión de flujo no aparece en Lax, 1954, sino en Rusanov, 1961).
Un ejemplo de la efectividad de usar un esquema de alta resolución se muestra en el diagrama opuesto, que ilustra la ecuación advectiva 1D , con una onda escalonada propagándose hacia la derecha. La simulación se llevó a cabo en una malla de 200 celdas, utilizando el esquema central de Kurganov y Tadmor con limitador Superbee y se utilizó RK-4 para la integración temporal. Este resultado de simulación contrasta extremadamente bien con los resultados de diferencia central de segundo orden y de primer orden en contra del viento que se muestran arriba. Este esquema también proporciona buenos resultados cuando se aplica a conjuntos de ecuaciones; consulte los resultados a continuación para este esquema aplicado a las ecuaciones de Euler. Sin embargo, se debe tener cuidado al elegir un limitador apropiado porque, por ejemplo, el limitador Superbee puede causar una nitidez poco realista para algunas ondas suaves.
El esquema puede incluir fácilmente términos de difusión, si están presentes. Por ejemplo, si el problema escalar unidimensional anterior se amplía para incluir un término de difusión, obtenemos
para lo cual Kurganov y Tadmor proponen la siguiente aproximación de diferencia central,
Dónde,
Los detalles completos del algoritmo ( versiones completa y semidiscreta ) y su derivación se pueden encontrar en el artículo original (Kurganov y Tadmor, 2000), junto con una serie de ejemplos unidimensionales y bidimensionales. También se puede encontrar información adicional en el artículo anterior relacionado de Nessyahu y Tadmor (1990).
Nota: Kurganov y Tadmor presentaron originalmente este esquema como un esquema de segundo orden basado en la extrapolación lineal . Un artículo posterior (Kurganov y Levy, 2000) demuestra que también puede formar la base de un esquema de tercer orden. En las secciones de reconstrucción parabólica y ecuación de Euler que aparecen a continuación, se muestran un ejemplo advectivo unidimensional y un ejemplo de ecuación de Euler de su esquema, utilizando una reconstrucción parabólica (tercer orden) .
Reconstrucción parabólica por partes
Es posible extender la idea de la extrapolación lineal a una reconstrucción de orden superior, y en el diagrama opuesto se muestra un ejemplo. Sin embargo, para este caso, los estados izquierdo y derecho se estiman mediante la interpolación de una ecuación diferencial de segundo orden sesgada en contra del viento. Esto da como resultado un esquema de reconstrucción parabólica que tiene una precisión de tercer orden en el espacio.
Seguimos el enfoque de Kermani (Kermani, et al., 2003) y presentamos un esquema de tercer orden con sesgo en contra del viento, donde los símbolos y nuevamente representan funciones dependientes del esquema (de las variables de borde de celda reconstruidas limitadas). Pero para este caso se basan en estados reconstruidos parabólicamente, es decir ,
y
Donde = 1/3 y,
y la función limitadora , es la misma que la anterior.
La reconstrucción parabólica es fácil de implementar y se puede utilizar con el esquema de Kurganov y Tadmor en lugar de la extrapolación lineal que se muestra arriba. Esto tiene el efecto de elevar la solución espacial del esquema KT a tercer orden. Funciona bien al resolver las ecuaciones de Euler, ver abajo. Este aumento en el orden espacial tiene ciertas ventajas sobre los esquemas de segundo orden para soluciones suaves, sin embargo, para choques es más disipativo - compare el diagrama opuesto con la solución anterior obtenida usando el algoritmo KT con extrapolación lineal y limitador Superbee. Esta simulación se llevó a cabo en una malla de 200 celdas usando el mismo algoritmo KT pero con reconstrucción parabólica. La integración temporal fue por RK-4, y la forma alternativa del limitador de van Albada, , se usó para evitar oscilaciones espurias.
Ejemplo: ecuaciones de Euler en 1D
Para simplificar, consideramos el caso unidimensional sin transferencia de calor y sin fuerza del cuerpo. Por lo tanto, en forma de vector de conservación, las ecuaciones generales de Euler se reducen a
dónde
y donde es un vector de estados y es un vector de flujos.
Las ecuaciones anteriores representan la conservación de la masa , el momento y la energía . Por lo tanto, hay tres ecuaciones y cuatro incógnitas: (densidad), (velocidad del fluido), (presión) y (energía total). La energía total está dada por:
donde representa la energía interna específica.
Para cerrar el sistema se requiere una ecuación de estado . Una que se adapta a nuestro propósito es
donde es igual a la relación de calores específicos del fluido.
Ahora podemos proceder, como se muestra arriba en el ejemplo 1D simple, obteniendo los estados extrapolados izquierdo y derecho para cada variable de estado. Por lo tanto, para la densidad obtenemos
dónde
De manera similar, para el momento y la energía total , la velocidad se calcula a partir del momento y la presión se calcula a partir de la ecuación de estado.
Una vez obtenidos los estados extrapolados limitados, procedemos a construir los flujos de borde utilizando estos valores. Con los flujos de borde conocidos, ahora podemos construir el esquema semidiscreto, es decir ,
La solución ahora puede proceder mediante integración utilizando técnicas numéricas estándar.
Lo anterior ilustra la idea básica del esquema MUSCL. Sin embargo, para una solución práctica de las ecuaciones de Euler, también se debe elegir un esquema adecuado (como el esquema KT anterior) para definir la función .
El diagrama opuesto muestra una solución de segundo orden al problema del tubo de choque de GA Sod (Sod, 1978) utilizando el esquema central de Kurganov y Tadmor (KT) de alta resolución anterior con extrapolación lineal y limitador de Ospre. Esto demuestra claramente la eficacia del enfoque MUSCL para resolver las ecuaciones de Euler. La simulación se llevó a cabo en una malla de 200 celdas utilizando el código Matlab (Wesseling, 2001), adaptado para utilizar el algoritmo KT y el limitador de Ospre . La integración temporal se realizó mediante un integrador SHK de cuarto orden (rendimiento equivalente a RK-4). Se utilizaron las siguientes condiciones iniciales (unidades del SI ):
presión izquierda = 100000 [Pa];
presión derecha= 10000 [Pa];
densidad izquierda = 1,0 [kg/m3];
densidad derecha = 0,125 [kg/m3];
longitud = 20 [m];
velocidad izquierda = 0 [m/s];
velocidad derecha = 0 [m/s];
duración = 0,01 [s];
lambda = 0,001069 (Δt/Δx).
El diagrama opuesto muestra una solución de tercer orden al problema del tubo de choque de GA Sod (Sod, 1978) utilizando el esquema central de Kurganov y Tadmor (KT) de alta resolución anterior pero con reconstrucción parabólica y limitador de van Albada. Esto ilustra nuevamente la efectividad del enfoque MUSCL para resolver las ecuaciones de Euler. La simulación se llevó a cabo en una malla de 200 celdas utilizando código Matlab (Wesseling, 2001), adaptado para utilizar el algoritmo KT con extrapolación parabólica y limitador de van Albada . La forma alternativa del limitador de van Albada, , se utilizó para evitar oscilaciones espurias. La integración temporal se realizó mediante un integrador SHK de cuarto orden. Se utilizaron las mismas condiciones iniciales.
Se han desarrollado otros esquemas de alta resolución que resuelven las ecuaciones de Euler con buena precisión. Algunos ejemplos de estos esquemas son:
Puede encontrar más información sobre estos y otros métodos en las referencias que aparecen a continuación. Puede encontrar una implementación de código abierto del esquema central de Kurganov y Tadmor en los enlaces externos que aparecen a continuación.
Kermani, MJ, Gerber, AG y Stockie, JM (2003), Predicción de humedad basada en termodinámica utilizando el esquema de Roe, 4.ª Conferencia de la Sociedad Aeroespacial Iraní , Universidad Tecnológica Amir Kabir, Teherán, Irán, 27 al 29 de enero. [1]
Kurganov, Alexander y Eitan Tadmor (2000), Nuevos esquemas centrales de alta resolución para leyes de conservación no lineales y ecuaciones de convección-difusión, J. Comput. Phys. , 160 , 241–282. [2]
Kurganov, Alexander y Doron Levy (2000), Un esquema central semidiscreto de tercer orden para leyes de conservación y ecuaciones de convección-difusión, SIAM J. Sci. Comput. , 22 , 1461–1488. [3]
Lax, PD (1954). Soluciones débiles de ecuaciones hiperbólicas no lineales y su cálculo numérico, Comm. Pure Appl. Math. , VII, págs. 159-193.
Leveque, RJ (2002). Métodos de volumen finito para problemas hiperbólicos , Cambridge University Press.
van Leer, B. (1979), Hacia el esquema diferencial conservativo último, V. Una secuela de segundo orden del método de Godunov, J. Com. Phys. ., 32 , 101–136.
Nessyahu, H. y E. Tadmor (1990), Diferenciación central no oscilatoria para leyes de conservación hiperbólica, J. Comput. Phys. , 87 , 408–463. [4].
Rusanov, VV (1961). Cálculo de la intersección de ondas de choque no estacionarias con obstáculos, J. Comput. Math. Phys. URSS , 1 , págs. 267-279.
Sod, GA (1978), Un estudio numérico de un choque cilíndrico convergente. J. Fluid Mechanics , 83 , 785–794.
Toro, EF (1999), Solucionadores de Riemann y métodos numéricos para dinámica de fluidos , Springer-Verlag.
Wesseling, Pieter (2001), Principios de dinámica de fluidos computacional , Springer-Verlag.
Lectura adicional
Hirsch, C. (1990), Cálculo numérico de flujos internos y externos , vol 2, Wiley.
Laney, Culbert B. (1998), Dinámica computacional de gases , Cambridge University Press.
Tannehill, John C., et al. (1997), Mecánica de fluidos computacional y transferencia de calor , 2.a edición, Taylor y Francis.
Enlaces externos
GEES – Código fuente abierto para resolver las ecuaciones de Euler utilizando el esquema central de Kurganov y Tadmor, escrito en Fortran (autor: Arno Mayrhofer)