You are currently browsing the category archive for the ‘Uncategorized’ category.

No deja de ser curioso que el método de Henyey o la matriz de Henyey estén totalmente vinculados a la resolución de las ecuaciones de estructura y evolución estelar. Se trata de un sistema de ODEs de primer orden con términos no lineales que modelan un problema de frontera. Obviamente, necesitamos de métodos numéricos para abordar su resolución. Sin embargo, no debe ser el único sistema de este tipo del que se necesite conocer su solución (¿o si que lo es?).

Uno de los métodos estándar para este tipo de problemas es el shooting method. Sin embargo, se argumenta que no es el mejor para este problema. Como ya hemos comentado en algún post, el problema de las primera derivadas es que en su discretización solo transportan información en uno de las dos direcciones. Además, resulta que en este problema la información conocida en las fronteras es parcial: conocemos la mitad de las condiciones en un extremo y la otra mitad en el otro (en concreto, se desconocen, por una parte, los valores de P y de T en el centro, y por otra, los valores de r y de l en la superficie). Ésto hace que las soluciones estén parametrizadas por los valores desconocidos. Y por si fuera poco, ademas de este desconocimiento parcial, existen factores con exponentes cuárticos, tanto en T como en r, que hacen que todo sea muy sensible a pequeñas variaciones.

En este caso, aunque dicen que no resulta práctico, se va integrando desde ambos extremos hasta un punto intermedio hasta que el fitting sea suave. Esto se podrá conseguir gracias a la variación de los parámetros libres (esta claro que se desconocen los valores de algunas de las funciones en el centro o en la superficie. Sin embargo, por cuestiones de simetria en el centro o de decaimiento en la superficie, igual si que se pueden imponer condiciones sobre las derivadas. O sea, aunque no conozcamos su valor Dirichlet, igual si que podemos conocer su valor Neumann, ¿no?).

Es aquí cuando entra en juego el método de Henyey. Vamos a tratar de entender el método con un ejemplo práctico, por lo que vamos a tratar de resolver el mismo problema que en este post, del que existe solución analítica, pero suponiendo que nuestro conocimiento en las fronteras es parcial, es decir, solo conoceremos uno de los dos valores en cada uno de los extremos y de manera complementaria. Recordemos que el sistema es

\frac{d}{dm}u = u^2 v, \,\,\,\,\,\, \frac{d}{dm}v = u v^2

en [0,1] donde ahora vamos a suponer que solo conocemos u(0)=1 en un extremo y v(1)=1/e en el otro. Los otros dos valores en los extremos se van a ir ajustando para conseguir el fit central. Los representaremos mediante u(1)=u_1 y v(0)=v_0. En el método del disparo esto habría que hacerlo manualmente, mientras que se supone que en Henyey está automatizado.

Anuncios

Como ya hemos comentado en algún post, MAESTRO es un código de libre distribución para fluidos con número de Mach bajo. Además, durante la instalación de mismo, ya comentamos que se necesita instalar previamente el repositorio AMReX que, entre otras cosas, suministra el mallado para AMR y un resolvedor Multigrid (MG).

La diferencia de resolución entre mallas del MG es un factor 2. Se encuentran el la carpeta F_MG dentro de LinearSolvers. Además, implementa dos tipos diferentes de solvers, uno centrado en nodos y otro centrado en celdas. No solo eso, también soporta mallas tridimensionales con tamaños de celda distintos en cada una de las dimensiones. Y para rematar, stencils en cruz de 5 (7) puntos, de 9 (12) puntos; y densos de 9 (27) puntos en 2D (3D). No esta nada mal.

Entre muchas otras cosas configurables, una de ellas es el tipo de smoother que se utilizará, teniendo disponibles tanto un Jacobi, MG_SMOOTHER_JACOBI, como un Red-Black Gauss-Seidel, MG_SMOOTHER_GS_RB.

Por último, basta recordar que nuestro Chebyshev-Jacobi es un método que acelera a Jacobi, por lo que puede utilizarse también en combinación con un Multigrid en esta parte de suavizado. Basta calcular los esquemas óptimos para las diferentes mallas, ordenarlos para evitar overflows y añadir un peso \omega delante del código correspondiente a Jacobi

smoother

y que se va leyendo de los esquemas óptimos, y que si lo fijamos a 1 recuperamos el esquema clásico que tenemos ahora implementado.

Se podría medir el impacto que tiene el utilizar este tipo de smoother en los métodos de proyección utilizados en problemas astrofísicos reales a números de Mach bajos.

 

Recordemos que hemos trabajado con las ecuaciones hidrodinámicas tanto en 1D, para estructura estelar, como en 3D, con Navier-Stokes.

Buscando las similitudes y diferencias en ambos casos, en esencia, parece que en 1D, en las ecuaciones de estructura estelar, no tenemos ecuaciones sobre las velocidades y todo se reduce a una especie (debido a la no linealidad) de proyección en la que todo depende, en esencia, de las fronteras (que además cumple principio del máximo, por estudios de unicidad que existen sobre las soluciones de las ecuaciones de estructura estelar).

Todo cambia en 3D cuando se toma la aproximación incompresible de las ecuaciones de la hidrodinámica, por cuestiones razonablemente físicas. La introducción de la restricción sobre el campo de velocidades, su divergencia nula, añade algo que con simples proyecciones no podemos arreglar. De hecho, la simple presencia de un campo de velocidades ya hace que, fundamentalmente, el problema parece no ser elíptico ni en el estado estacionario (o por lo menos no se cumple el principio del máximo), puesto que no esta garantizada la solución por los valores en la frontera, sino por condiciones de radiación en la misma de tipo Sommerfeld. Si se opta por evolucionar el mismo, el papel del campo de velocidades inicial juega un papel determinante.

Existen muchas aproximaciones de las ecuaciones de la hidrodinámica para bajas velocidades. Todos estos métodos tienen en común un ecuación de restricción sobre el campo de velocidades.

MAESTRO es un código de libre distribución que permite, mediante métodos pseudoincompresibles y una reformulación de las ecuaciones del momento que conserva la energía, trabajar con este tipo de fluidos.

Realiza dos proyecciones para garantizar la restricción sobre las velocidades: una a mitad de paso de tiempo centrada en interfaces y otra en el nuevo paso de tiempo centrada en celdas.

Para instalar el código hay que bajarse tres repositorios: el propio MAESTRO junto con dos auxiliares llamados AMRex, para el mallado adaptativo, y Microphysics para la microfísica. Una vez creadas las variables de entorno necesarias para que MAESTRO los encuentre, ya podemos compilar todo y empezar a utilizarlo.

El código viene con muchos tests de ejemplo. Simplemente he ejecutado el reacting_buble y automáticamente genera una serie de ficheros para este problema determinado.

Entre todos los ficheros de salida, aparecen los plt***, que en realidad son directorios,  que contienen los datos de diferentes variables en los distintos periodos de tiempo de la simulación. Estos datos se pueden visualizar de diferentes maneras. Una de ellas es VisIt, que es la por la que he optado puesto que ya lo he utilizado en muchas otras ocasiones.

En los siguientes dos gráficos se muestran el número de Mach M y la velocidad v_y en el instante final de la simulación. Necesitamos cargar el fichero headers para que VisIt sepa como tratarlos y los hemos mostrado mediante pseudocolor con dos color tables distintas para cada una de las variables

Acabamos de ver en el post anterior que las ecuaciones de Navier-Stokes para fluidos incompresibles sin fuerzas externas en \mathbb{R}^3 quedan:

\begin{cases} \frac{\partial}{\partial t} v_x + \frac{\partial}{\partial x} p = \nu \Delta v_x - v_x \frac{\partial}{\partial x} v_x - v_y \frac{\partial}{\partial y} v_x - v_z \frac{\partial}{\partial z} v_x \\ \frac{\partial}{\partial t} v_y + \frac{\partial}{\partial y} p = \nu \Delta v_y - v_x \frac{\partial}{\partial x} v_y - v_y \frac{\partial}{\partial y} v_y - v_z \frac{\partial}{\partial z} v_y \\ \frac{\partial}{\partial t} v_z + \frac{\partial}{\partial z} p = \nu \Delta v_z - v_x \frac{\partial}{\partial x} v_z - v_y \frac{\partial}{\partial y} v_z - v_z \frac{\partial}{\partial z} v_z \\ \frac{\partial}{\partial x} v_x + \frac{\partial}{\partial y} v_y + \frac{\partial}{\partial z} v_z = 0 \end{cases}

Sabemos que empezamos con v^0 = (v_x^0, v_y^0, v_z^0) tal que \nabla \cdot v^0 = 0, por lo que ya cumple la última condición.

Para estado estacionario, podríamos encontrar con cualquier integrador un campo escalar p tal que verificara las tres primeras ecuaciones:

\begin{cases} \frac{\partial}{\partial x} p = \nu \Delta v_x^0 - v_x^0 \frac{\partial}{\partial x} v_x^0 - v_y^0 \frac{\partial}{\partial y} v_x^0 - v_z^0 \frac{\partial}{\partial z} v_x^0 \\ \frac{\partial}{\partial y} p = \nu \Delta v_y^0 - v_x^0 \frac{\partial}{\partial x} v_y^0 - v_y^0 \frac{\partial}{\partial y} v_y^0 - v_z^0 \frac{\partial}{\partial z} v_y^0 \\ \frac{\partial}{\partial z} p = \nu \Delta v_z^0 - v_x^0 \frac{\partial}{\partial x} v_z^0 - v_y^0 \frac{\partial}{\partial y} v_z^0 - v_z^0 \frac{\partial}{\partial z} v_z^0 \end{cases}

Llamamos a este campo que acabamos de encontrar p^0. Ahora ya podemos evolucionar un paso de tiempo, también con un integrador del ODEs, de manera que:

\begin{cases} \frac{\partial}{\partial t} v_x = \nu \Delta v_x^0 - v_x^0 \frac{\partial}{\partial x} v_x^0 - v_y^0 \frac{\partial}{\partial y} v_x^0 - v_z^0 \frac{\partial}{\partial z} v_x^0 - \frac{\partial}{\partial x} p^0 \\ \frac{\partial}{\partial t} v_y = \nu \Delta v_y^0 - v_x^0 \frac{\partial}{\partial x} v_y^0 - v_y^0 \frac{\partial}{\partial y} v_y^0 - v_z^0 \frac{\partial}{\partial z} v_y^0 - \frac{\partial}{\partial y} p^0 \\ \frac{\partial}{\partial t} v_z = \nu \Delta v_z^0 - v_x^0 \frac{\partial}{\partial x} v_z^0 - v_y^0 \frac{\partial}{\partial y} v_z^0 - v_z^0 \frac{\partial}{\partial z} v_z^0 - \frac{\partial}{\partial z} p^0 \end{cases}

Y encontramos una nueva v*^1. Pero este campo vectorial de velocidades no es de divergencia nula, o sea, que tenemos:

\nabla \cdot v*^1 \neq 0

Sin embargo, por Helmholtz, podemos proyectar. De esta manera, nos quedaría

v^1 = v*^1 - \nabla \varphi

donde

\Delta \varphi = \nabla \cdot v*^1

con las propiedades de frontera adecuadas. Llegados a este punto, y siguiendo la misma estrategia que antes, podemos encontrar p^1. De esta manera, cualquier pareja (v^n, p^n) cumple las ecuaciones de Navier-Stokes en estado estacionario. Pararemos la evolución cuando converjamos a una solución, es decir, cuando la distancia entre dos campos de velocidades sucesivos sea menor que una tolerancia dada.

Así describe Charles L. Fefferman (Medalla Fields en 1978) de manera oficial uno de los siete Problemas del Milenio propuesto por el Clay Mathematics Institute en el año 2000.

Las ecuaciones de Navier-Stokes describen el movimiento de un fluido en \mathbb{R}^n, \,n \in \{2,3\}.  Las incógnitas son el campo vectorial de velocidades v(t,x) = (v_i(t,x))_{1 \leq i \leq n} \in \mathbb{R}^n y el campo escalar de presiones p(t,x) \in \mathbb{R}, definidos para x \in \mathbb{R}^n y t \geq 0.

El problema se restringe al caso de fluidos incompresibles sobre todo \mathbb{R}^n. Estas ecuaciones, deducidas a partir de elementos de fluidos, quedan como;

\nu \sum_{j=1}^{n} \frac{\partial^2}{\partial x_j^2} v_i - \frac{\partial}{\partial x_i} + f_i = \frac{\partial}{\partial t} v_i + \sum_{j=1}^{n} v_j \frac{\partial}{\partial x_j} v_i

que no es mas que la segunda ley de Newton F = m a donde f(t,x) = (f_i(t,x))_{1 \leq i \leq n} es una fuerza externa que actúa junto con las fuerzas internas de presión y de rozamiento, siendo \nu el coeficiente de viscosidad; y como:

\sum_{i=1}^{n} \frac{\partial}{\partial x_i} v_i = 0

que expresa la incompresibilidad del fluido. Necesitamos además la condición inicial

v(0,x) = v^0(x) \in C^\infty(\mathbb{R}^n)

con x \in \mathbb{R}^n y \sum_{i=1}^{n} \frac{\partial}{\partial x_i} v^0_i = 0, o sea, ya incompresible.

El problema del que se pide demostración vuelve a restringirse una vez mas y solo tiene en cuenta el caso f(t,x)=0. Además, para evitar trabajar con todo \mathbb{R}^n, existe una formulación del mismo en el que solo se buscan soluciones periódicas, es decir, tales que

v^0(x+e_j)=v^0(x) para 1 \leq j \leq n

donde e_j es el j^{\mathrm{th}} vector unitario en \mathbb{R}^n. Para este caso, basta asumir que v^0 es suave i diremos que tenemos una solución si

v, p \in C^{\infty}([0,\infty) \times \mathbb{R}^n)

donde, además, tenemos que

v(t,x) = v(t,x+e_j) \in [0,\infty) \times \mathbb{R}^n

Finalmente, como queremos que v(t,x) no crezca de manera descontrolada cuando |x| \rightarrow \infty, también deseamos que la solución tenga energía acotada:

\int_{\mathbb{R}^n} |u(t,x)|^2 dx < K para todo t \geq 0

Como ya hemos comentado en algún post anterior, recuperamos Euler simplemente haciendo que \nu =0 y el problema también está abierto y es importante.

Para mas información, aquí.

Recordemos que Navier-Stokes no es mas que Euler donde hemos añadido viscosidad \mu, o equivalentemente pero al revés, Euler no es mas que Navier-Stokes donde hemos eliminado la viscosidad (\mu = 0). La viscosidad, de hecho, solo modifica las ecuaciones del momento.

¿Y como quedan estas en el caso de fluidos incompresibles (asumimos, como hasta ahora, que \rho = 1) cuando se encuentran en estado estacionario? Pues son las siguientes (expandidas, que es como las necesitamos al tratarlas numéricamente). Para el momento, nos queda:

\frac{\partial}{\partial x} (v_x^2 + p) + \frac{\partial}{\partial y} v_x v_y + \frac{\partial}{\partial z} v_x v_z = \mu \Delta v_x

\frac{\partial}{\partial x} v_x v_y + \frac{\partial}{\partial y} (v_y^2 + p) + \frac{\partial}{\partial z} v_y v_z = \mu \Delta v_y

\frac{\partial}{\partial x} v_x v_z + \frac{\partial}{\partial y} v_y v_z + \frac{\partial}{\partial z} (v_z^2+p) = \mu \Delta v_z

Y la de continuidad queda:

\frac{\partial}{\partial x} v_x + \frac{\partial}{\partial y} v_y + \frac{\partial}{\partial z} v_z = 0

Una forma de proceder sería la indicada en Solución numérica a las ecuaciones de Euler incompresible n-D (estado estacionario).

Sin embargo, el proceder con Navier-Stokes era para tener operadores elípticos. Podríamos pensar en resolver numéricamente el sistema

\begin{cases} \Delta v_x = \frac{1}{\mu} [ \frac{\partial}{\partial x} (v_x^2 + p) + \frac{\partial}{\partial y} v_x v_y + \frac{\partial}{\partial z} v_x v_z ] \\ \Delta v_y = \frac{1}{\mu} [\frac{\partial}{\partial x} v_x v_y + \frac{\partial}{\partial y} (v_y^2 + p) + \frac{\partial}{\partial z} v_y v_z] \\ \Delta v_z = \frac{1}{\mu} [\frac{\partial}{\partial x} v_x v_z + \frac{\partial}{\partial y} v_y v_z + \frac{\partial}{\partial z} (v_z^2+p) ] \\ 0 = \frac{\partial}{\partial x} v_x + \frac{\partial}{\partial y} v_y + \frac{\partial}{\partial z} v_z \end{cases}

pero existen una serie de problemas.

Para empezar, la última ecuación sigue siendo con derivadas primeras. Además, las tres ecuaciones que aparentemente nos deberían llevar a encontrar v=(v_x, v_y, v_z) presentan, en sus respectivos términos fuente, la también desconocida p, de manera que, tal y como esta escrito, no podemos abordar el problema numéricamente (¿qué valores ponemos en p). Si en la ecuación de continuidad apareciera p, podríamos escribir ésta en función del campo de velocidades y sustituir en las ecuaciones del momento, pero tampoco es el caso, pues en esta última solo se explicita una condición que debe de cumplir v. Por tanto, la p también tiene que aparecer en la parte izquierda de las ecuaciones

\begin{cases} \Delta v_x - \frac{\partial}{\partial x} p  = \frac{1}{\mu} [ \frac{\partial}{\partial x} v_x^2 + \frac{\partial}{\partial y} v_x v_y + \frac{\partial}{\partial z} v_x v_z ] \\ \Delta v_y - \frac{\partial}{\partial y} p= \frac{1}{\mu} [\frac{\partial}{\partial x} v_x v_y + \frac{\partial}{\partial y} v_y^2 + \frac{\partial}{\partial z} v_y v_z] \\ \Delta v_z - \frac{\partial}{\partial z} p = \frac{1}{\mu} [\frac{\partial}{\partial x} v_x v_z + \frac{\partial}{\partial y} v_y v_z + \frac{\partial}{\partial z} v_z^2 ] \\ 0 = \frac{\partial}{\partial x} v_x + \frac{\partial}{\partial y} v_y + \frac{\partial}{\partial z} v_z \end{cases}

La discretización del término de la izquierda de la primera ecuación quedaría:

\frac{v_{x_{i-1,j,k}}-2v_{x_{i,j,k}}+v_{x_{i+1,j,k}}}{h_x^2} + \frac{v_{x_{i,j-1,k}}-2v_{x_{i,j,k}}+v_{x_{i,j+1,k}}}{h_y^2}

+ \frac{v_{x_{i,j,k-1}}-2v_{x_{i,j,k}}+v_{x_{i,j,k+1}}}{h_z^2} - \frac{p_{i+1,j,k}-p_{i-1,j,k}}{2h_x}

Puesto que en el esquema iterativo solo puede quedar una incógnita en la parte izquierda de la iteración ¿Estamos consiguiendo algo? Parece que no… Ya que, como hemos dicho, una de las ecuaciones solo incluye velocidades, tendríamos que reconvertir alguna de las del Laplaciano en una iteración para la p. Sin embargo, junto con el orden de la discretización para incluir el punto actual,  el problema que tenemos es que la presión aparece con primera derivadas…

Como ya comentamos en el post Hidrodinámica topológica y las ecuaciones de Euler incompresible en el caso estacionario, los campos de Beltrami (o force-free fields en otras áreas) mas famosos son los campos ABC.

Recordemos que estos campos son de la forma

v(x,y,z) =

(A \sin z + C \cos y, B \sin x + A \cos z, C \sin y + B \cos x)

con A, B, C \in \mathbb{R} y que son solución de las ecuaciones de Euler tridimensionales para fluidos incompresibles en estado estacionario:

\omega = \nabla \times v =

( \partial_y (C \sin y + B \cos x) - \partial_z (B \sin x + A \cos z),

\partial_z (A \sin z + C \cos y) - \partial_x (C \sin y + B \cos x),

\partial_x (B \sin x + A \cos z) - \partial_y (A \sin z + C \cos y) ) =

(C \cos y - (-A \sin z), A \cos z - (- B \sin x), B \cos x - (-C \sin y)) = v

\nabla \cdot v =

\partial_x (A \sin z + C \cos y) +  \partial_y (B \sin x + A \cos z) + \partial_z (C \sin y + B \cos x) = 0

Para ver mas el aspecto que tiene, podemos pintar el campo particular con A = B = C = 1 y una sección del mismo:

 

 

Dos consideraciones al respecto del campo. Por una parte, la autosemejanza que posee a muchas y diferentes escalas, por poner un ejemplo:

 

ABC1_2d_2

 

La autosemejanza es una propiedad intrínseca de la geometría fractal y ésta, a su vez, esta relacionada con el comportamiento caótico.

Y la otra, que para poder aplicar Helmholtz a un campo necesitamos de un decaimiento mas rápido que 1/r^2 en el infinito.

Ya hemos visto en Hidrodinámica topológica y las ecuaciones de Euler incompresible en el caso estacionario. que las soluciones en el caso tridimensional pueden ser tan complicadas como queramos.

Generalizando el caso unidimensional visto en Euler incompresible en 1D, parece que podríamos llegar a construir soluciones numéricas, también muy complicadas, de la ecuación en cuestión.

Recapitulemos: tenemos el sistema con n+1 PDEs:

(v \cdot \nabla)v + \nabla P = 0, \,\,\,\,\,\, \nabla \cdot v = 0

con las n+1 incognitas v=(v_1, ..., v_n) y P. ¿Cómo podemos construir soluciones? Básicamente aplicaremos dos veces el método de proyección.

Dado un campo de velocidades v^* cualquiera, por el teorema de Helmholtz siempre podemos escribirlo como una cantidad libre de divergencia, a la que vamos a llamar v, ya que será nuestra solución, y el gradiente de un escalar:

v^* = v + \nabla \varphi

Tomando la divergencia de la ecuación, y teniendo en cuenta que por construcción tenemos que \nabla \cdot v = 0, nos queda la ecuación elíptica

\Delta \varphi = \nabla \cdot v^*

que, una vez resuelta, al conocer \varphi nos permite el cálculo

v = v^* - \nabla \varphi

Acabamos de construir un campo de velocidades v que satisface la ecuación de restricción \nabla \cdot v = 0.

Para encontrar el valor de P, procedemos de la misma manera. Utilizamos la ecución (reescrita) que nos queda

\nabla P = - (v \cdot \nabla) u

de la que podemos volver a tomar divergencias, volviendo a encontrar otra ecuación tipo Poisson con fuente conocida

\Delta P = - \nabla \cdot \left [ (v \cdot \nabla) v \right ]

Y ya estaria. Lo que falta aquí es encontrar las condiciones de frontera apropiadas en ambos casos, que para el caso elíptico es fundamental.

Como ya hemos visto en Hidrodinámica incompresible, las ecuaciones de Euler para éste caso quedan:

\frac{\partial}{\partial t}v + (v \cdot \nabla)v + \nabla P = 0, \,\,\,\,\,\, \nabla \cdot v =0

De manera que, para el caso estacionario, que es el que nos interesa, tenemos:

(v \cdot \nabla)v + \nabla P = 0, \,\,\,\,\,\, \nabla \cdot v =0

Tenemos un sistema de n+1 PDEs con n+1 incognitas: el campo de velocidades v=(v_1, ... ,v_n) y la presión interna del fluido P.

Para el caso n=3, existen teoremas en el campo de la hidrodinámica topológica que nos dicen como son las lineas de corriente (trayectorias del campo v) y de vorticidad (trayectorias del campo \omega = \nabla \times v) de sus soluciones. De hecho, en función de si los campos de velocidad y de vorticidad son o no colineales, tendremos dos grandes familias de soluciones en función de la complejidad topología de sus lineas de corriente.

Aunque existen otras hipótesis adicionales, lo fundamental es que para el caso de que v \times w \neq 0, es decir, que los campos de velocidad y vorticidad  no son colineales, el Teorema de Estructura de Arnold nos dice que las trayectorias del fluido son bastante ordenadas y se distribuyen sobre superficies regulares: toros o cilindros.

Para el caso en que u y \omega son colineales, el Teorema de Realización de Enciso y Peralta-Salas nos dice que ahora las lineas de corriente y vorticidad pueden tener topologías arbitrariamente complejas: cualquier enlace (link), una colección de nudos (knots) libre de intersecciones, localmente finito es realizable.

El ejemplo paradigmático de campos con esta última propiedad en el espacio euclídeo tridimensional son  los campos de Beltrami, campos tales que \nabla \times v = \lambda v, siendo los campos ABC, campos con

v(x_1,x_2,x_3) =

(A \sin x_3 + C \cos x_2, B \sin x_1 + A \cos x_3, C \sin x_2 + B \cos x_1)

y A, B, C \in \mathbb{R}, los mas estudiados.

En ambos casos, las demostraciones utilizan la forma equivalente de las ecuaciones de Euler

v \times (\nabla \times v) = \nabla B, \,\,\,\,\,\, \nabla \cdot v =0

siendo B=P + \frac{|v|^2}{2}. Se utiliza el hecho de que

(v \cdot \nabla)v = (\nabla \times v) \times v + \nabla \left ( \frac{|v|^2}{2} \right )

(Para mas información: aquí y aquí)

El sistema de ecuaciones para el caso estacionario que vimos en el post Hidrodinámica incompresible queda, en 1D y coordenadas cartesianas, de la siguiente manera:

\frac{d}{dx}(v^2 + p) = 0,\,\,\,\,\,\,\,\,\frac{d}{dx}v = 0

o, equivalentemente, con la misma estructura que los sistemas vistos en Sistema de ODEs acopladas no lineal:

\begin{cases} \frac{d}{dx}p = -2v \\ \frac{d}{dx}v =0 \end{cases}

En este caso, tenemos solución analítica, pues la segunda ecuación obliga a que v(x) = k_1, mientras que la primera nos dice que p(x) = -v(x)^2 + k_2, que sustituyendo el valor de v, nos queda p(x) = -k_1^2 + k_2 = k_3. De esta forma, tenemos que para densidad constante, la velocidad y la presión también lo son.

Al tratar de aplicar nuestro método numérico tenemos un problema con el hecho de que una de las fuentes ya sea cero, ya que al derivarla numéricamente de nuevo con nuestro método, estamos perdiendo información.

Existen muchas formulaciones de las ecuaciones de la hidrodinámica para bajas velocidades (low Mach number (M) flows, con M = v / c, donde u es la velocidad del flujo y c es la velocidad del sonido en el medio). La aproximación mas sencilla es el límite M \rightarrow 0, que se corresponde con la hidrodinámica incompresible.

En hidrodinámica incompresible, junto con las ecuaciones del movimiento (Navier-Stokes), aparece una ecuación de restricción sobre la velocidad:

\mathrm{div}\,v = 0

de manera que equilibra de manera instantanea los flujos, filtrando las ondas de sonido.

Junto a la restricción en la divergencia del campo de velocidades, tenemos que resolver las ecuaciones de la hidrodinámica incompresible, que son un conjunto de ecuaciones de tipo advección. Por ejemplo, para el caso de densidad constante, tendríamos:

v_t = - (v \cdot \nabla) v - \nabla p

Si nos centramos en el caso estacionario, lo que nos quedan son las ecuaciones:

\begin{cases} (v \cdot \nabla) v + \nabla p = 0 \\ \mathrm{div}\,v = 0 \end{cases}

Resulta que es un caso concreto de aplicabilidad del Teorema de Realización demostrado en su publicación Knots and links in steady solutions of the Euler equation en Annals of Mathematics (por cierto, revista numero uno a nivel mundial para publicaciones en matemáticas, es decir, solo publican los mejores del mundo) por Alberto Enciso y Daniel Peralta-Salas, cuyo abstract dice:

Given any possibly unbounded, locally finite link, we show that there exists a smooth diffeomorphism transforming this link into a set of stream (or vortex) lines of a vector field that solves the steady incompressible Euler equation in \mathbb{R}^3. Furthermore, the diffeomorphism can be chosen arbitrarily close to the identity in any C^r norm.

Una manera de poder imponer la condición r(0)=0 es trabajar con centros de celda en lugar de con centros de malla. Para esto, lo que hacemos es añadir dos celdas fantasma a ambos lados de la malla original de manera que en la interfaz entre ambas el valor corresponda con la condición que queremos imponer.

Si queremos imponer una condición de Dirichlet homogénea basta hacer

u_{-1} = -u_0

mientras que si lo que queremos es que sea Neumann homogénea, basta con imponer

u_{-1} = u_0

La explicación a partir de las discretizaciones. Por ejemplo, si queremos el valor de la derivada en el punto u_{-1/2} (donde estaría nuestro r=0) sea a, la aproximación a segundo orden es

\frac{u_{0}-u_{-1}}{2h} = a

de donde recuperamos la expresión que habíamos puesto simplemente despejando e imponiendo a=0 (el razonamiento es idéntico mutatis mutandis para Dirichlet).

Con esto, la gráfica que obtenemos para r(0)=0 es

stellarA00005

Hemo

En las ecuaciones de evolución estelar (The differential equations of Stellar Evolution) tenemos dos sectores claramente diferenciados: las dos primeras ecuaciones hacen referencia a la parte mecánica y las dos últimas a la parte termoenergética. Y el único vinculo entre ellas se realiza a través de la densidad \rho(P,T).

Así pues, podemos centrarnos solo en la parte mecánica si asumimos densidades que no dependan de T. Un ejemplo un tanto artificial de esto lo conseguimos al asumir \rho = cte. Y éste, a su vez, no sería mas que un caso particular de \rho(m).

Por simplificarlo al máximo, vamos a asumir \rho = \frac{1}{G}, y vamos a cambiar las escalas para que \frac{1}{4 \pi G} sea la unidad. Con esto, las ecuaciones quedan simplemente:

\begin{cases} \frac{d}{dm} r = \frac{1}{r^2} \\ \frac{d}{dm} P = - \frac{m}{r^4} \end{cases}

donde 0 \leq m \leq M. Si imponemos r(0)=0, \frac{d}{dm}r|_{m=M}=0 (?), \frac{d}{dm}P|_{m=0}=0 y P(M)=0 en las fronteras, entonces nuestro código nos da:

star1

Al aparecer r(m) en el denominador de ambos términos fuentes, existen algunos problemas para r=0. En realidad, hemos tenido que imponer r(0)=0.35 en la frontera e inicializar el guess inicial a 1.

En todos los últimos posts hemos visto el buen comportamiento de nuestro método. Sin embargo, recordemos que estamos derivando numéricamente las fuentes de nuestras ecuaciones para poder tener esquemas con derivadas segundas. ¿Qué ocurre si estas fuentes no son derivables (o ni siquiera continuas)?

Supongamos que queremos encontrar la solución para

\frac{d}{dm}u = \mathrm{sgn}(m-1/2)

en el dominio m \in [0,1], cuya solución es u(m) = |m-1/2|.

En este ejemplo, tenemos un término fuente que no es ni continuo. Sin embargo, nuestro método tratará de derivarlo numéricamente (y recordemos que, en términos de distribuciones, si es posible hacerlo siendo \delta_{1/2}(m)).

La primera gráfica del final del post muestra el término fuente, en azul, la solución analítica, en naranja, y la solución que encontramos numéricamente, en verde.

Si utilizamos la función escalón (desplazada), H(m-1/2) como fuente para que la solución sea la función rampa (desplazada), R(m-1/2), lo que obtenemos aparece en la segunda gráfica.

Finalmente, si tomamos como fuente una delta de Dirac en el centro del dominio, \delta_{1/2}(m), el resultado es la última de las gráficas.

Por lo tanto, vemos que el método es lo suficientemente robusto para tratar también con este tipo de términos fuente, aunque las soluciones que vamos a encontrar van a ser versiones suaves de las soluciones exactas.

Dado que la discontinuidad esta en la fuente, y sabemos su localización, si pudiéramos intuir el valor de la solución en ella, podríamos dividir el problema frontera en dos y la solución sería mucho mejor.

 

 

En Sistema de ODEs acopladas no lineal hemos visto como resolver problemas de valor frontera para sistemas de ODEs de primer orden con fuentes suaves utilizando esquemas tipo Jacobi.

¿Podemos acelerar estos métodos tipo Jacobi mediante la utilización de esquemas óptimos de tipo Chebyshev-Jacobi utilizados para ecuaciones elípticas? Calculamos dicho esquema óptimo para el problema resuelto en el post del principio. Solo necesitamos saber que todas las fronteras son de tipo Dirichlet, que la malla tiene 50 puntos y que queremos bajar el residual hasta 0.000001.

Las gráficas muestran la evolución del residual, proporcional al error, en función del número de iteraciones. Se muestra la evolución de éste para las diferentes ejecuciones exteriores del bucle que trata la no linealidad. A la izquierda, el esquema Jacobi clásico. A la derecha, el esquema CJM.

 

De esta manera, confirmamos que la técnica también funciona para este tipo de sistemas. Recordemos que en el ejemplo solo trabajamos con N=50 (consiguiendo esquemas  \approx 20 veces mas rápidos) pero que la mejora de CJM con respecto Jacobi es mayor a medida que crece el valor de N.

 

Siguiendo con las pruebas numéricas, suponemos el siguiente sistema de ODEs:

\begin{cases} \frac{d}{dm}u = u^2 v  \\  \frac{d}{dm}v = u v^2 \end{cases}

que es un sistema sencillo  con fuentes no lineales acopladas entre si del que conocemos su solución analítica:

u(m) = e^m, \,\,\,\, v(m) = e^{-m}

Tratamos de resolverlo tal y como hicimos en el post  Derivando sistemas de ODEs: derivando en ambos lados del sistema:

\begin{cases} \frac{d^2}{dm^2}u = \frac{d}{dm}(u^2 v)  \\  \frac{d^2}{dm^2}v = \frac{d}{dm}(u v^2) \end{cases}

De esta manera, utilizamos una discretización en diferencias finitas de las derivadas segundas de los lados izquierdos de las ecuaciones del sistema; mientras que derivamos numéricamente los términos fuente del lado derecho del sistema. Además, por ser no lineal, tenemos un bucle externo adicional en el que vamos actualizando los términos fuente cada vez que encontramos una solución para las fuentes que teníamos en el paso anterior. Las fronteras son las mismas que teniamos antes de derivar:

u(0)=1,\,\, u(1)=e,\,\, v(0)=1,\,\, v(1)=1/e

La siguiente secuencia de imágenes muestra como nos acercamos correctamente a la solución tras muy pocas iteraciones (tres) de la parte debida a la no linealidad:

Tal y como acabamos de ver en el post anterior, podemos derivar en ambos lados de una ODE para incrementar el orden de las derivadas. Además, hemos realizado la derivación del término fuente de manera numérica.

¿Qué pasa si en lugar de una única ODE tenemos un sistema de ODEs? Pues que todo sigue funcionando como se puede apreciar en la siguiente figura (el sistema de ecuaciones que se está resolviendo es trivialmente inferible):

system

Finalmente, si el sistema no es lineal, añadimos una iteración exterior que, dado que estamos asumiendo suavidad, convergerá en un número finito de iteraciones.

Otra opción que tenemos, si la fuente de la ecuación es lo suficientemente suave, es aplicar lo que comentamos en el post Derivando derivadas.

Para el mismo problema resuelto en Probando esquemas tipo Euler, lo que nos queda es \frac{\partial^2}{\partial m^2}f = \frac{\partial}{\partial m} \cos(m) con las mismas condiciones de frontera f(0)=0 y f(\pi) = 0.

Ahora tenemos un problema donde es aplicable un método iterativo, puesto que tenemos derivadas segundas. Realizamos la derivada del termino fuente de manera numérica. Las soluciones que obtenemos son (para 7 y 32 puntos respectivamente):

Volvemos a recuperar la solución f(m) = \sin (m) y con muchos menos puntos pues tenemos orden cuadrático.

Tal y como comenté hace pocos posts, las discretizaciones de ecuaciones con primeras derivadas llevan a esquemas tipo Euler (que, por tanto, podrían sustituirse por otros de mas alto orden, tipo Runge-Kutta por ejemplo).

Si resolvemos numéricamente el problema \frac{d}{dm}f = \cos(m) en el intervalo [0,\pi] con f(0) = 0 y f(\pi)=0, obtenemos las siguientes soluciones en función de si utilizamos esquemas hacia adelante, hacia atrás o mixto empezando por ambos lados (solución analítica en azul y numérica en amarillo):

El último es el único que asegura los valores conocidos de la solución en la frontera, dejando el error para el centro del dominio. Además, éste se reduce con el número de puntos de la malla utilizada (10, 50 y 1000 puntos):

Así pues, estábamos en lo cierto y, efectivamente, todo funciona como comentamos en el post Esquemas tipo Euler para ODEs en problemas de frontera

Partiendo del mismo problema tratado en otro post

\frac{d}{dm}f = g

con f(0) = a y f(M)=b, suponiendo que 0 \leq m \leq M. Suponiendo que las funciones sean lo suficientemente suaves, entonces podemos derivar a ambos lados de la expresion

\frac{d^2}{dm^2} f = \frac{d}{dm}g

Podriamos derivar analíticamente \frac{d}{dm}g. Sin embargo, si discretizamos, y utilizamos una derivación numérica del lado derecho de la expresion del mismo orden que la discretización, nos queda

\frac{f_{i-1}-2f_i+f_{i+1}}{h^2} = \frac{g_{i+1}-g_{i-1}}{2h}

Los valores de la frontera son condiciones sobre la solución, por lo que podemos utilizar las mismas. Solo hemos cambiado la manera de obtener ésta. ¿Es todo esto correcto?

 

Tratando de discretizar el problema

\frac{d}{dm}f = k_1

sobre una malla de cuatro puntos con f(0) = k1 y f(3) = 3k_1+k_2

nos encontramos con el esquema

f_i = h k_1 + f_{i-1}

si utilizamos la derivada hacia atrás, y con el esquema

f_i = f_{i+1} - h k1

si utilizamos la derivada hacia adelante.

Observamos que la información fluye de una frontera hacia la otra o viceversa. Si empezamos con el valor de la frontera f(0) nos llega f(3)=4k_1 \neq 3k_1+k_2 y si empezamos con la frontera f(3) = 3k_1+k_2 con el otro esquema, nos llega f(0)=k2 \neq k1.

En realidad, estamos obteniendo un esquema de Euler para ODEs. Son problemas de valor inicial. Allí solo sabemos el valor en t_0, en una de las fronteras, y el mayor error se produce en el tiempo final t_f, la otra frontera. Como sabemos el orden del método, sabemos el orden de este error producido. Ahora, como es un problema de frontera, sabemos el valor en ambos lados. Podríamos empezar dos esquemas de este tipo, uno desde cada frontera, y parar a la mitad, donde habría una discontinuidad, pero del orden del error del método de Euler. Sin embargo, ¿como se realizan estos recorridos en mas de una dimensión? (¿O mes mas sencillo propagar todo y realizar la media?)

Supongamos que tenemos la ecuación

\frac{d}{dx} = f(u,x)

y que queremos hacer un estudio de su estabilidad de von Neumann. Para empezar, necesitamos la discretización con factor de relajación y, además, asumimos que el termino fuente no afecta al mismo. Lo que nos queda es

u_{i}^{k+1} = (1-\omega) u_i^k +  \omega u_{i-1}^k

de donde tenemos que

\varepsilon_{i}^{k+1} = \varepsilon_i^k - \omega ( \varepsilon_i^k - \varepsilon_{i-1}^k)

Como lo que nos interesa es el factor de amplificación G = \frac{\varepsilon_i^{k+1}}{\varepsilon_i^k}, lo que nos queda es

G = 1 - \omega (1-e^{-i k h})

Siendo k los números de onda. Así como en las discretizaciones para el Laplaciano derivaba en expresiones donde podíamos combinar las exponenciales para formar funciones trigonométricas, aquí no pasa eso.

Si añadimos discretizaciones de mas puntos, no centradas para obtener el punto actual, seguimos teniendo el mismo problema. Por ejemplo, para el caso de no centrada de orden dos, lo que nos queda es

G = 1- \omega (-2 + 2e^{-i k h} + \frac{1}{2}- \frac{1}{2} e^{-i k 2h})

 

En el post anterior hemos discretizado un sistema de ecuaciones donde únicamente aparecen primeras derivadas.

Para que quede claro lo que queremos decir, vamos a trabajar con la ecuación genérica

\frac{d}{dx}u = f(u,x)

Ya en el punto de la discretización necesitamos comentar dos cosas. Por una parte, dado que queremos escribir el método iterativamente, necesitamos que en la discretización aparezca un término con el punto actual. En caso contrario, no podemos despejar la variable en el punto actual. De esta manera, para la discretización de la primera derivada a orden uno hacia atrás, tenemos:

\frac{u_{i} - u_{i-1}}{h} = f(u,x) \Leftrightarrow u_i = h f(u,x) + u_{i-1}

Y en diferencias hacia adelante, nos encontramos con:

\frac{u_{i+1} - u_i}{h} = f(u,x) \Leftrightarrow u_i = u_{i+1} - h f(u,x)

Sin embargo, si utilizamos la aproximación a orden dos centrada, tenemos

\frac{u_{i+1} - u_{i-1}}{2h} = f(u,x)

nos encontramos con que no podemos despejar u_i para implementar la iteración ya que este termino no aparece en la ecuación.

Pero no solo eso. Además, la discretización que utilicemos debería dar lugar a una matriz con todos los autovalores con el mismo signo, para que una iteración tipo Jacobi pueda aplicarse. La primera discretización y la segunda cumplen esta propiedad (dan lugar a una matriz con 1 en la diagonal y -1 arriba o abajo de la diagonal en función de la discretización escogida). La última, sin embargo, da lugar a una matriz con todos los autovalores complejos (0 en la diagonal, 1 y -1 arriba y abajo de esta).

 

Asumiendo una estrella con simetría esférica y en complete equilibrium (hydrostatic equilibrium + thermally adjusted), las ecuaciones diferenciales básicas que nos darán su estructura son:

\frac{\partial}{\partial m}r = \frac{1}{4 \pi r^2 \rho}

\frac{\partial}{\partial m}P = - \frac{G m}{4 \pi r^4}

\frac{\partial}{\partial m}l = \varepsilon_n + \varepsilon_\nu

\frac{\partial}{\partial m}T = - \frac{G m T}{4 \pi r^4 P} \nabla

donde \nabla = \frac{d \ln T}{d \ln P}. Queremos encontrar funciones solución para la distancia radial r(m), para la presión P(m), para la energía l(m) y para la temperatura T(m) donde 0 \leq m \leq M, siendo M la masa total de la estrella. Necesitamos conocer los valores en la frontera de las funciones incógnita, esto es, r(0), P(0), l(0), T(0),  r(M), P(M), l(M), T(M). A través de la relación de operadores

\frac{\partial}{\partial m} = \frac{1}{4 \pi r^2 \rho} \frac{\partial}{\partial r}

deberíamos de recuperar las ecuaciones de la hidrodinámica en steady-state expresadas en su forma habitual.

Podemos discretizarlas a primer orden:

\frac{r_{i} - r_{i-1}}{\Delta m} = \frac{1}{4 \pi r_{i/2}^2 \rho_i}

\frac{P_{i} - P_{i-1}}{\Delta m} = - \frac{G m_i}{4 \pi r_{i/2}^4}

\frac{l_{i} - l_{i-1}}{\Delta m} = (\varepsilon_n)_i + (\varepsilon_\nu)_i

\frac{T_{i} - T_{i-1}}{\Delta m} = - \frac{G m_{i} T_{i/2}}{4 \pi r_{i/2}^4 P_{i/2}} \nabla_{i/2}

Si despejamos las variables en el punto i de la discretización, obtenemos:

r_{i}^{*} = \Delta m \frac{1}{4 \pi r_{i/2}^2 \rho_i} + r_{i-1}

P_{i}^{*} = - \Delta m \frac{ G m_i}{4 \pi r_{i/2}^4} + P_{i-1}

l_{i}^{*} = \Delta m [ (\varepsilon_n)_i + (\varepsilon_\nu)_i] + l_{i-1}

T_{i}^{*} = - \Delta m \frac{G m_{i} T_{i/2}}{4 \pi r_{i/2}^4 P_{i/2}} \nabla_{i/2} + T_{i-1}

y con esto, podemos obtener esquemas con factor de relajación:

r_{i}^{k} = (1 - \omega) r_{i}^{k-1} + \omega r_{i}^{*}

P_{i}^{k} = (1 - \omega) P_{i}^{k-1} + \omega P_{i}^{*}

l_{i}^{k} = (1 - \omega) l_{i}^{k-1} + \omega l_{i}^{*}

T_{i}^{k} = (1 - \omega) T_{i}^{k-1} + \omega T_{i}^{*}

 

diciembre 2017
L M X J V S D
« Ago    
 123
45678910
11121314151617
18192021222324
25262728293031