Método de Godunov

La versión actual de la página aún no ha sido revisada por colaboradores experimentados y puede diferir significativamente de la versión revisada el 9 de abril de 2016; las comprobaciones requieren 4 ediciones .

El método de Godunov es una implementación de esquemas de conteo que se pueden usar para calcular flujos dinámicos de gas con discontinuidades en los parámetros dentro del dominio computacional. Este esquema fue propuesto por S. K. Godunov en 1959. El método de Godunov es una variante del método del volumen de control . Los flujos por las caras laterales se determinan a partir de la solución del problema del decaimiento de una discontinuidad arbitraria . Vamos a explicar con un ejemplo.

Ejemplo

Considere la construcción del método numérico Godunov de primer orden de precisión usando el ejemplo de resolver el sistema de ecuaciones de dinámica de gas inestable unidimensional, escrito en forma divergente :

{\begin{cases}{\begin{array}{lll}{\frac {\parcial \rho }{\parcial t))+{\frac {\parcial \rho u}{\parcial x))&=&0 \\{\frac {\parcial (\rho u)}{\parcial t))+{\frac {\parcial (p+\rho u^{2})}{\parcial x))&=&0\\{ \frac {\E parcial}{\t parcial}}+{\frac {\u parcial(E+p)}{\x parcial}}&=&0\end{matriz}}\end{casos}}

Aquí:

$\rho$ - densidad
$pags$ - presión
$tu$ - velocidad
$mi$ es la energía total por unidad de volumen

Darse cuenta de:

La primera ecuación expresa la ley de conservación de la masa.
La segunda ecuación expresa la ley de conservación de la cantidad de movimiento
La tercera ecuación expresa la ley de conservación de la energía. $E=e+{\frac {\rho u^{2}}{2}}$ $e=\rho\varepsilon,$ $\varepsilon =\varepsilon (p,\rho ),$ $\varepsilon ={\frac {1}{\gamma -1}}\cdot {\frac {p}{\rho }}$ , donde es el exponente adiabático $\gama$

Forma diferencial

El sistema inicial se puede escribir en una forma más compacta:

{\frac {\q parcial}{\t parcial))+{\frac {\f parcial}{\x parcial))=0

dónde:

$q$ es el vector de variables conservativas $q=\left({\begin{matriz}{c}\rho \\\rho \,u\\E\end{matriz}}\right)$
$F$ es el vector de flujo $f=\left({\begin{matriz}{c}\rho \,u\\p+\rho \,u^{2}\\u(E+p)\end{matriz}}\right)$

Forma integral

En lugar de la forma diferencial de las ecuaciones, derivamos una nueva forma integral de las ecuaciones, más adecuada para representar una solución débil . Aquí, una solución débil es una función generalizada definida por igualdades integrales obtenidas de las correspondientes ecuaciones diferenciales y las condiciones iniciales del problema. Para ello, seleccionamos un cierto volumen de control e integramos el sistema de ecuaciones sobre este volumen. Aplicamos el teorema de Stokes generalizado a la integral de divergencia obtenida (para dos variables independientes, este será el teorema de Green y la fórmula de Ostrogradsky-Gauss en el espacio tridimensional). En este caso, introducimos la dirección de eludir el contorno en sentido antihorario . $\Omega$

Por separado, considerando la ecuación de continuidad , obtenemos:

\iint \limits _{{\Omega }}\left({\frac {\parcial \rho }{\parcial t}}+{\frac {\parcial \rho u}{\parcial x}}\right)d \,xd\,t=\oint \limits _{{\parcial \Omega }}\rho dx-\rho udt

Para todo el sistema de ecuaciones

\iint \limits _{{\Omega }}\left({\frac {\parcial q}{\parcial t}}+{\frac {\parcial f}{\parcial x}}\right)d\,xd \,t=0\Leftrightarrow \oint \limits _{{\parcial \Omega }}(qdx-fdt)=0

Escribiendo el sistema en forma desarrollada:

{\begin{cases}{\begin{array}{lll}\oint \limits _{{\parcial \Omega }}(\rho \;d\,x-\rho u\;d\,t)&= &0\\\oint \limits _{{\parcial \Omega }}(\rho u\;d\,x-(p+\rho u^{2})\;d\,t)&=&0\\\ punto \limits _{{\parcial \Omega }}(E\;d\,xu(p+E)\;d\,t)&=&0\\\end{matriz}}\end{casos}}

Aproximación

Se realiza la transición de la forma diferencial de escribir el sistema de ecuaciones original a la forma integral. La forma integral se escribe como igualdad a cero de las integrales sobre el contorno (el límite del volumen de control seleccionado) de los vectores de flujos y variables conservativas. Representamos la integral de contorno como la suma de las integrales sobre las secciones (intervalos) 1-2 , 2-3 , 3-4 , 4-1 del volumen de control en la figura (que aún no está disponible) y en cada sección aproximamos la integral usando el método de los rectángulos como producto del integrando en el centro del intervalo por la longitud del intervalo de integración:

{\begin{matriz}{c}q_{{12}}\cdot (x_{{2}}-x_{{1}})-f_{{12}}\cdot (t_{{2}}-t_ {{1}})+\\\qquad q_{{23}}\cdot (x_{{3}}-x_{{2}})-f_{{23}}\cdot (t_{{3}} -t_{{2}})+\\\qquad \qquad q_{{34}}\cdot (x_{{4}}-x_{{3}})-f_{{34}}\cdot (t_{ {4}}-t_{{3}})+\\\qquad \qquad \qquad q_{{41}}\cdot (x_{{1}}-x_{{4}})-f_{{41} }\cdot (t_{{1}}-t_{{4}})=0\end{matriz}}\Rightarrow q_{{34}}=q_{{12}}-{\frac {\Delta t} {\Deltax}}(f_{{23}}-f_{{41}})

teniendo en cuenta las igualdades que son válidas para el volumen de control construido sobre la grilla computacional cartesiana:

$x_{3}-x_{2}=0$
$x_{1}-x_{4}=0$
$t_{2}-t_{1}=0$
$t_{4}-t_{3}=0$

Es más:

$x_{2}-x_{1}=\Delta x$
$x_{4}-x_{3}=-\Delta x$
$t_{3}-t_{2}=\Delta t$
$t_{1}-t_{4}=-\Delta t$

encuentre los valores del vector de variables conservativas en el intervalo 3-4 perteneciente a la nueva capa:

q_{{34}}=q_{{12}}-{\frac {\Delta t}{\Delta x}}(f_{{23}}-f_{{41}})\qquad \Rightarrow \qquad q_ {{j}}^{{n+1}}=q_{{j}}^{{n}}-{\frac {\Delta t}{\Delta x}}(f_{{j+{\frac { 1}{2))))-f_{{j-{\frac{1}{2)))))

En este caso, los valores con índices semienteros denotan los flujos de valores almacenados a través de los límites de la celda computacional durante el tiempo o los flujos a través de las caras laterales ( 2-3 y 4-1 ) del control volumen. Si la velocidad del flujo se dirige en la misma dirección que la normal exterior a la cara lateral, entonces el flujo es negativo , es decir, sale del volumen de control y viceversa.

Expandido:

{\begin{casos}{\begin{matriz}{lll}\rho_{{j}}^{{n+1}}&=&\rho_{{j}}^{{n}}-{ \frac {\Delta t}{\Delta x))((\rho u)_{{j+{\frac {1}{2))))-(\rho u)_{{j-{\frac { 1}{2}}}})\\(\rho u)_{{j}}^{{n+1}}&=&(\rho u)_{{j}}^{{n}} -{\frac {\Delta t}{\Delta x))((p+\rho u^{2})_{{j+{\frac {1}{2))))-(p+\rho u^{ 2})_{{j-{\frac {1}{2)))))\\E_{{j}}^{{n+1}}&=&E_{{j}}^{{n} }-{\frac {\Delta t}{\Delta x))((u(p+E))_{{j+{\frac {1}{2))))-(u(p+E)) _{{j-{\frac {1}{2}}}})\end{matriz}}\end{casos}}

Los flujos por las caras laterales, y se determinan a partir de la solución del problema del decaimiento de una discontinuidad arbitraria . $f_{{j+{\frac{1}{2))))$ $f_{{j-{\frac{1}{2))))$

Declaración de condiciones de contorno

Una característica del establecimiento e implementación de condiciones de contorno en los métodos de volumen de control (incluido el método de Godunov) es la necesidad de especificar o calcular flujos a través de la cara del volumen de control que coincide con el límite del dominio computacional. Para la primera y la última celda de la capa computacional, es necesario determinar los flujos de masa, cantidad de movimiento y energía a través de las caras.

A menudo, se introducen celdas computacionales "virtuales" para establecer las condiciones de contorno. Para ello, se introduce una celda adicional más a la izquierda de la primera celda y a la derecha de la última celda, en cada una de las cuales se especifican dichos parámetros de flujo para que se modelen los flujos requeridos en la cara lateral al resolver el Riemann problema

Tipos de procedimientos de límites

Todas las suposiciones se hacen en relación con el borde izquierdo.

Muro rígido fijo

La condición principal es la ausencia de flujo del flujo másico de gas a través de la frontera, lo que corresponde a la condición de velocidad de flujo cero en la cara dada En la celda virtual, se deben configurar los siguientes parámetros de flujo: $tu=0$

{\begin{casos}{\begin{matriz}{lcl}p_{w}&=&p_{1}\\\rho _{w}&=&\rho _{1}\\u_{w}&= &-u_{1}\\\end{matriz}}\end{casos}}

"w" - parámetros en una celda virtual
"1" - parámetros en la primera celda

Los parámetros de flujo obtenidos en el problema de decaimiento de la discontinuidad en la cara lateral dan cuenta de un flujo de masa cero a través de esta cara.

Embalse de capacidad ilimitada

Matemáticamente, este caso corresponde a establecer el valor de presión en la cara . La tasa de flujo de entrada se puede determinar mediante la fórmula ${\ tilde {P}}$

{\tilde {U}}=u_{1}+{\frac {{\tilde {P}}-p_{1}}{c_{1}}}

Donde:

si , entonces ${\tilde {P}}>p_{1}$ $c_{1}={\frac {{\sqrt {\rho _{1}\cdot ((\gamma +1){\tilde {P))+(\gamma +1)p_{1})))} {2}}$
si , entonces ${\tilde {P}}<p_{1}$ $c_{1}={\frac {a_{1}\rho _{1}(\gamma -1)\left(1-{\frac {P}{p_{1))}\right)}{2\ gamma \left(1-\left({\frac {P}{p_{1))}\right)^{({\frac {\gamma -1}{2\gamma ))))\right)))$

Flujo supersónico entrante

Deje que el guión bajo denote los parámetros del flujo supersónico, luego si , entonces ${\bar {U}}>{\bar {c}}={\sqrt {\frac {\gamma {\bar {P}}}{\bar {\rho }}}}$

{\begin{casos}{\begin{array}{lcl}p_{w}&=&{\bar {P}}\\\rho _{w}&=&{\bar {\rho }}\\ u_{w}&=&{\bar {U}}\\\end{matriz}}\end{casos}}

Escapando del flujo supersónico

En este caso, los siguientes parámetros de flujo se establecen en la celda virtual:

{\begin{casos}{\begin{matriz}{lcl}p_{w}&=&p_{1}\\\rho _{w}&=&\rho _{1}\\u_{w}&= &u_{1}\\\end{matriz}}\end{casos}}

Selección de opciones de malla

El paso de la cuadrícula computacional a lo largo de la coordenada de tiempo en el método de Godunov se puede determinar a partir del criterio de estabilidad de Courant-Friedrichs-Levy . Con respecto al régimen considerado, esta condición se formula de la siguiente manera:

Las ondas que surgen en el problema del decaimiento de una discontinuidad arbitraria en el punto no deben llegar a tiempo a las caras laterales y distorsionar la solución autosimilar . $j+{\frac{1}{2}}$ $\Delta t$ $j+{\frac {3}{2}}$ $j-{\frac {1}{2}}$

La aplicación de este principio conduce a las siguientes relaciones:

\Delta t_{j}=r{\frac {\Delta x_{j}}{\max \left(|D_{{j+{\frac {1}{2}}}}^{{L}}|, |D_{{j-{\frac {1}{2}}}}^{{L}}|\right)))

dónde

$D_{{j+{\frac {1}{2}}}}^{{L}}$ es el valor de la velocidad de la onda más a la izquierda en el decaimiento de la discontinuidad;
$D_{{j-{\frac{1}{2}}}}^{{R}}$ es el valor de la velocidad de la onda más a la derecha en el decaimiento de la discontinuidad;

Como resultado, tomamos:

\Delta t=\min_{j}\Delta t_{j}

Literatura

Solución numérica de problemas multidimensionales de dinámica de gases. Álbum / editor Godunov S. K. . — M .: Nauka, 1976. — 400 p. - 6500 copias.
Samarsky A.A., Popov Yu.P. Métodos diferenciales para la resolución de problemas de dinámica de gases. - M. : Nauka, 1992. - 2470 ejemplares.

Enlaces

Métodos para resolver ecuaciones diferenciales

Métodos de cuadrícula

Métodos de elementos finitos	Método de elementos finitos método Galerkin Método de Galerkin discontinuo Método de elementos finitos multiescala método multired
Otros metodos	Método de diferencias finitas Método de volumen finito método de Godunov Método de elementos de contorno

Métodos sin cuadrícula