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

A mi casi-código SPH le he añadido la libreria Voro++ de Chris Rycroft y ésta va calculando las correspondientes teselaciones de Voronoi tridimensionales correspondientes a las partículas siguiendo movimientos pendulares (sobre planos z=cte) .

Pulsar sobre la imagen para empezar la animación.animate_1

Anuncios

Definimos a los kernels como funciones del tipo:

W_{ab}=W(\boldsymbol{r}_a - \boldsymbol{r}_b,h),

donde a es la partícula en la que está centrada la función y b es una partícula dentro del soporte compacto de la función kernel, controlado éste último por h, la smoothing length (longitud de suavizado).

En este post básicamente pretendo aclarar lo que significa \nabla_a W_{ab} cuando, por ejemplo, tenemos definido el kernel como:

W(q) = \alpha_D \exp (-q^2) con 0 \leq q \leq 2.

Para empezar, \alpha_D es una constante de dimensionalidad, por lo que la fórmula está escrita de manera compacta y sirve para cualquier dimensión. Además, tenemos que q = \frac{r}{h}, siendo r la distancia ente las partículas, por lo que:

r = |\boldsymbol{r}_a - \boldsymbol{r}_b| =^{(3D)} \sqrt{(x_a-x_b)^2 + (y_a-y_b)^2 + (z_a-z_b)^2}.

Si fijamos la posición de la partícula a, la función que nos da la distancia de esta a cualquier punto dentro del soporte compacto es:

r_a (\boldsymbol{r}) = |\boldsymbol{r}_a-\boldsymbol{r}| =^{(3D)} \sqrt{(x_a-x)^2 + (y_a-y)^2 + (z_a-z)^2},

siendo q_a lo mismo añadiendo el factor h.

Por lo tanto, en este caso tenemos, en tres dimensiones y donde b es una partícula en una posición arbitraria (x,y,z):

\nabla_a W_{ab}(q) =^{(3D)} (\partial_x W_{ab}(q_a), \partial_y W_{ab}(q_a), \partial_z W_{ab}(q_a)) =

= \alpha_D \exp(-q^2) (-2q) (\partial_x (q_a), \partial_y (q_a), \partial_z (q_a)) = \alpha_D \exp(-q^2) (-2q) \nabla_a q_a

donde:

\nabla_a q_a = \frac{-1}{h r_a} (x_a-x,y_a-y,z_a-z).

De la misma manera, si tenemos:

W(q) = \alpha_D \begin{cases} 1-\frac{3}{2} q^2 + \frac{3}{4} q^3, 0 \leq q < 1\\ \frac{1}{4} (2-q)^3, 1 \leq q < 2 \\ 0, q \geq 2 \end{cases}

entonces:

\nabla_a W_{ab}(q) = \alpha_D \begin{cases} (-3q + \frac{9}{4}q^2) \nabla_a q_a, 0 \leq q < 1 \\ -\frac{3}{4} (2-q)^2 \nabla_a q_a, 1 \leq q < 2 \\ 0, q \geq 2 \end{cases}

Así pues:

\nabla W(q) = \frac{d}{dq} W(q) \nabla q.

 

En la discretización que hicimos teníamos dos sistemas acoplados, uno para las X^i y otro para las \beta^i. Procedemos ahora a desacoplarlos.

Para empezar, tomamos la divergencia (plana) del sistema:

\Delta X^i = 8 \pi f^{ij} S^*_j - \frac{1}{3}\mathcal{D}^i \mathcal{D}_j X^j

y, teniendo en cuenta que \mathcal{D} conmuta con \Delta (métrica plana), tenemos:

\Delta (\mathcal{D}_i X^i) = 8 \pi \mathcal{D}^j S^*_j - \frac{1}{3} \Delta (\mathcal{D}_j X^j),

por lo que:

\Delta (\mathcal{D}_i X^i) = \frac{3}{4} 8 \pi \mathcal{D}^j S^*_j.

De esta manera, si definimos \Theta_X := \mathcal{D}_i X^i, nos queda:

\Delta \Theta_X = \frac{3}{4} 8 \pi \mathcal{D}^j S^*_j = 6 \pi (\partial_x S^*_x + \partial_y S^*_y +\partial_z S^*_z ),

que discretizado queda:

\frac{(\Theta_X)_{i-1,j,k}-2(\Theta_X)_{i,j,k}+(\Theta_X)_{i+1,j,k}}{h_x^2} +

\frac{(\Theta_X)_{i,j-1,k}-2(\Theta_X)_{i,j,k}+(\Theta_X)_{i,j+1,k}}{h_y^2} +

\frac{(\Theta_X)_{i,j,k-1}-2(\Theta_X)_{i,j,k}+(\Theta_X)_{i,j,k+1}}{h_z^2} =

= 6 \pi (\partial_x S^*_x + \partial_y S^*_y +\partial_z S^*_z )_{i,j,k} ,

donde inicialmente:

(S^*_a)_{i,j,k} = (\psi^6)_{i,j,k}\rho_{i,j,k}h_{i,j,k}w^2_{i,j,k}(v_a)_{i,j,k},

(\partial_x S^*_x + \partial_y S^*_y +\partial_z S^*_z )_{i,j,k} =

\frac{(S^*_x)_{i+1,j,k}-(S^*_x)_{i-1,j,k}}{2h_x} + \frac{(S^*_x)_{i,j+1,k}-(S^*_x)_{i,j-1,k}}{2h_y} + \frac{(S^*_x)_{i,j,k+1}-(S^*_x)_{i,j,k-1}}{2h_z}

y que es lineal.

El primer sistema acoplado de ecuaciones quedaría ahora:

\partial_{xx} X^x + \partial_{yy} X^x + \partial_{zz} X^x = 8 \pi S^*_x - \frac{1}{3} \partial_x \Theta_X \approx

\approx \frac{X^x_{i-1,j,k}-2X^x_{i,j,k}+X^x_{i+1,j,k}}{h_x^2} + \frac{X^x_{i,j-1,k}-2X^x_{i,j,k}+X^x_{i,j+1,k}}{h_y^2} + \frac{X^x_{i,j,k-1}-2X^x_{i,j,k}+X^x_{i,j,k+1}}{h_z^2} =

= 8 \pi (S^*_x)_{i,j,k} - \frac{1}{3} (\partial_x \Theta_X)_{i,j,k},

¡que vuelve a ser lineal!

Continuamos con:

\partial_{xx} X^y + \partial_{yy} X^y + \partial_{zz} X^y = 8 \pi S^*_y - \frac{1}{3} \partial_y \Theta_X \approx

\approx \frac{X^y_{i-1,j,k}-2X^y_{i,j,k}+X^y_{i+1,j,k}}{h_x^2} + \frac{X^y_{i,j-1,k}-2X^y_{i,j,k}+X^y_{i,j+1,k}}{h_y^2} + \frac{X^y_{i,j,k-1}-2X^y_{i,j,k}+X^y_{i,j,k+1}}{h_z^2} =

= 8 \pi (S^*_y)_{i,j,k} - \frac{1}{3} (\partial_y \Theta_X)_{i,j,k}

y, finalmente:

\partial_{xx} X^z + \partial_{yy} X^z + \partial_{zz} X^z = 8 \pi S^*_z - \frac{1}{3} \partial_z \Theta_X \approx

\approx \frac{X^z_{i-1,j,k}-2X^z_{i,j,k}+X^z_{i+1,j,k}}{h_x^2} + \frac{X^z_{i,j-1,k}-2X^z_{i,j,k}+X^z_{i,j+1,k}}{h_y^2} + \frac{X^z_{i,j,k-1}-2X^z_{i,j,k}+X^z_{i,j,k+1}}{h_z^2} =

= 8 \pi (S^*_z)_{i,j,k} - \frac{1}{3} (\partial_z \Theta_X)_{i,j,k},

donde calculamos al principio:

(\partial_x \Theta_X)_{i,j,k} = \frac{(\Theta_X)_{i+1,j,k}-(\Theta_X)_{i-1,j,k}}{2h_x}

(\partial_y \Theta_X)_{i,j,k} = \frac{(\Theta_X)_{i,j+1,k}-(\Theta_X)_{i,j-1,k}}{2h_y}

(\partial_z \Theta_X)_{i,j,k} = \frac{(\Theta_X)_{i,j,k+1} - (\Theta_X)_{i,j,k-1}}{2h_z}

A continuación, discretizamos las siguientes ecuaciones:

\hat{A}^{xx} = 2 \partial_x X^x - \frac{2}{3} (\partial_x X^x + \partial_y X^y + \partial_z X^z) \approx

\approx \frac{2}{3}\frac{X^x_{i+1,j,k}-X^x_{i-1,j,k}}{h_x} -\frac{1}{3} \frac{X^y_{i,j+1,k}-X^y_{i,j-1,k}}{h_y} - \frac{1}{3} \frac{X^z_{i,j,k+1}-X^z_{i,j,k-1}}{h_z}) = \hat{A}^{xx}_{i,j,k},

\hat{A}^{xy} = \hat{A}^{yx}= \partial_x X^y + \partial_y X^x \approx

\approx \frac{X^y_{i+1,j,k}-X^y_{i-1,j,k}}{2h_x} + \frac{X^x_{i,j+1,k}-X^x_{i,j-1,k}}{2h_y} = \hat{A}^{xy}_{i,j,k} = \hat{A}^{yx}_{i,j,k},

\hat{A}^{xz} = \hat{A}^{zx} = \partial_x X^z + \partial_z X^x \approx

\approx \frac{X^z_{i+1,j,k}-X^z_{i-1,j,k}}{2h_x} + \frac{X^x_{i,j,k+1}-X^x_{i,j,k-1}}{2h_z} = \hat{A}^{xz}_{i,j,k} = \hat{A}^{zx}_{i,j,k},

\hat{A}^{yy} = 2 \partial_y X^y - \frac{2}{3} (\partial_x X^x + \partial_y X^y + \partial_z X^z) \approx

\approx -\frac{1}{3}\frac{X^x_{i+1,j,k}-X^x_{i-1,j,k}}{h_x} +\frac{2}{3} \frac{X^y_{i,j+1,k}-X^y_{i,j-1,k}}{h_y} - \frac{1}{3} \frac{X^z_{i,j,k+1}-X^z_{i,j,k-1}}{h_z}) = \hat{A}^{yy}_{i,j,k},

\hat{A}^{yz} = \hat{A}^{zy} = \partial_y X^z + \partial_z X^y \approx

\approx \frac{X^z_{i,j+1,k}-X^z_{i,j-1,k}}{2h_y} + \frac{X^y_{i,j,k+1}-X^y_{i,j,k-1}}{2h_z} = \hat{A}^{yz}_{i,j,k} = \hat{A}^{zy}_{i,j,k},

\hat{A}^{zz} = 2 \partial_z X^z - \frac{2}{3} (\partial_x X^x + \partial_y X^y + \partial_z X^z) \approx

\approx -\frac{1}{3}\frac{X^x_{i+1,j,k}-X^x_{i-1,j,k}}{h_x} -\frac{1}{3} \frac{X^y_{i,j+1,k}-X^y_{i,j-1,k}}{h_y} + \frac{2}{3} \frac{X^z_{i,j,k+1}-X^z_{i,j,k-1}}{h_z}) = \hat{A}^{zz}_{i,j,k}.

Por tanto, la siguiente ecuación:

\Delta \psi = -2 \pi \psi^{-1} E^* - \psi^{-7} \frac{(\hat{A}^{xx})^2+(\hat{A}^{yy})^2+(\hat{A}^{zz})^2+2(\hat{A}^{xy})^2+2(\hat{A}^{xz})^2+2(\hat{A}^{yz})^2}{8}

queda:

\approx \frac{\psi_{i-1,j,k}-2\psi_{i,j,k}+\psi_{i+1,j,k}}{h_x^2} + \frac{\psi_{i,j-1,k}-2\psi_{i,j,k}+\psi_{i,j+1,k}}{h_y^2} + \frac{\psi_{i,j,k-1}-2\psi_{i,j,k}+\psi_{i,j,k+1}}{h_z^2} =

=-2 \pi \psi^{-1}_{i,j,k} E^*_{i,j,k} -

- \frac{\psi^{-7}_{i,j,k}}{8} ( (\hat{A}^{xx}_{i,j,k})^2+(\hat{A}^{yy}_{i,j,k})^2+(\hat{A}^{zz}_{i,j,k})^2+2(\hat{A}^{xy}_{i,j,k})^2+2(\hat{A}^{xz}_{i,j,k})^2+2(\hat{A}^{yz}_{i,j,k})^2 ) ,

con:

\partial_{\psi_{i,j,k}} F(\psi_{i,j,k}) = -2 ( \frac{1}{h_x^2} + \frac{1}{h_y^2} + \frac{1}{h_z^2} ) -2 \pi \psi_{i,j,k}^{-2} E^*_{i,j,k} -

- \frac{7}{8} \psi^{-8}_{i,j,k} ( (\hat{A}^{xx}_{i,j,k})^2+(\hat{A}^{yy}_{i,j,k})^2+(\hat{A}^{zz}_{i,j,k})^2+2(\hat{A}^{xy}_{i,j,k})^2+2(\hat{A}^{xz}_{i,j,k})^2+2(\hat{A}^{yz}_{i,j,k})^2 ),

donde:

E^*_{i,j,k} = \psi^{6}_{i,j,k} (D_{i,j,k}+\tau_{i,j,k})

y la ecuación:

\Delta (\alpha\psi) = (\alpha \psi) (2 \pi \psi^{-2} (E^*+2S^*) +

+ \frac{7}{8} \psi^{-8} ((\hat{A}^{xx})^2+(\hat{A}^{yy})^2+(\hat{A}^{zz})^2+2(\hat{A}^{xy})^2+2(\hat{A}^{xz})^2+2(\hat{A}^{yz})^2) )

como:

\approx \frac{(\alpha\psi)_{i-1,j,k} - 2(\alpha\psi)_{i,j,k}+(\alpha\psi)_{i+1,j,k}}{h_x^2} +

+ \frac{(\alpha\psi)_{i,j-1,k}-2(\alpha\psi)_{i,j,k}+(\alpha\psi)_{i,j+1,k}}{h_y^2} +

+ \frac{(\alpha\psi)_{i,j,k-1}-2(\alpha\psi)_{i,j,k}+(\alpha\psi)_{i,j,k+1}}{h_z^2} =

= (\alpha \psi)_{i,j,k} (2 \pi \psi^{-2}_{i,j,k} (E^*_{i,j,k}+2S^*_{i,j,k}) +

+ \frac{7}{8} \psi^{-8}_{i,j,k} ((\hat{A}^{xx}_{i,j,k})^2+(\hat{A}^{yy}_{i,j,k})^2+(\hat{A}^{zz}_{i,j,k})^2+2(\hat{A}^{xy}_{i,j,k})^2+2(\hat{A}^{xz}_{i,j,k})^2+2(\hat{A}^{yz}_{i,j,k})^2) ),

donde:

\partial_{(\alpha \psi)_{i,j,k}} F((\alpha \psi)_{i,j,k}) = -2 ( \frac{1}{h_x^2} + \frac{1}{h_y^2} + \frac{1}{h_z^2} ) - 2 \pi \psi^{-2}_{i,j,k} (E^*_{i,j,k}+2S^*_{i,j,k}) +

- \frac{7}{8} \psi^{-8}_{i,j,k} ((\hat{A}^{xx}_{i,j,k})^2+(\hat{A}^{yy}_{i,j,k})^2+(\hat{A}^{zz}_{i,j,k})^2+2(\hat{A}^{xy}_{i,j,k})^2+2(\hat{A}^{xz}_{i,j,k})^2+2(\hat{A}^{yz}_{i,j,k})^2) )

con:

S^*_{i,j,k} = \psi^6_{i,j,k}(\rho_{i,j,k}h_{i,j,k}(w^2_{i,j,k}-1) + 3 p_{i,j,k}).

Finalmente, tenemos el otro sistema acoplado:

\Delta \beta^i = \mathcal{D}_j(2 \alpha \psi^{-6} \hat{A}^{ij}) - \frac{1}{3} \mathcal{D}^i(\mathcal{D}_j \beta^j),

con el que procedemos de igual manera que con las X^i:

\Delta(\mathcal{D}_i \beta^i) = \mathcal{D}_i (\mathcal{D}_j (2 \alpha \psi^{-6} \hat{A}^{ij})) - \frac{1}{3} \Delta (\mathcal{D}_i \beta^i),

de manera que:

\Delta \Theta_\beta = \frac{3}{4} \mathcal{D}^i (\mathcal{D}_j (2 \alpha \psi^{-6} \hat{A}^{ij})) =

\frac{3}{2}(\partial_{xx}(\alpha \psi^{-6} \hat{A}^{xx}) + \partial_{yy}(\alpha \psi^{-6} \hat{A}^{yy}) + \partial_{zz}(\alpha \psi^{-6} \hat{A}^{zz}),

con:

\Theta_\beta := \mathcal{D}_i \beta^i,

que discretizada queda:

\frac{(\Theta_\beta)_{i-1,j,k}-2(\Theta_\beta)_{i,j,k}+(\Theta_\beta)_{i+1,j,k}}{h_x^2} +

\frac{(\Theta_\beta)_{i,j-1,k}-2(\Theta_\beta)_{i,j,k}+(\Theta_\beta)_{i,j+1,k}}{h_y^2} +

\frac{(\Theta_\beta)_{i,j,k-1}-2(\Theta_\beta)_{i,j,k}+(\Theta_\beta)_{i,j,k+1}}{h_z^2} =

\frac{3}{2}((\partial_{xx}(\alpha \psi^{-6} \hat{A}^{xx}))_{i,j,k} + (\partial_{yy}(\alpha \psi^{-6} \hat{A}^{yy}))_{i,j,k} + (\partial_{zz}(\alpha \psi^{-6} \hat{A}^{zz})_{i,j,k}),

De esta manera, tenemos:

\Delta \beta^x = \partial_x (2 \alpha \psi^{-6} \hat{A}^{xx}) + \partial_y (2 \alpha \psi^{-6} \hat{A}^{xy}) + \partial_z (2 \alpha \psi^{-6} \hat{A}^{xz}) - \frac{1}{3} \partial_x \Theta_\beta \approx

\approx \frac{\beta^x_{i-1,j,k}-2\beta^x_{i,j,k}+\beta^x_{i+1,j,k}}{h_x^2} + \frac{\beta^x_{i,j-1,k}-2\beta^x_{i,j,k}+\beta^x_{i,j+1,k}}{h_y^2} + \frac{\beta^x_{i,j,k-1}-2\beta^x_{i,j,k}+\beta^x_{i,j,k+1}}{h_z^2} =

= (\partial_x (2 \alpha \psi^{-6} \hat{A}^{xx}))_{i,j,k} + (\partial_y (2 \alpha \psi^{-6} \hat{A}^{xy}))_{i,j,k} + (\partial_z (2 \alpha \psi^{-6} \hat{A}^{xz}) )_{i,j,k} -

- \frac{1}{3} (\partial_x \Theta_\beta)_{i,j,k}.

De la misma manera:

\Delta \beta^y = \partial_x (2 \alpha \psi^{-6} \hat{A}^{yx}) + \partial_y (2 \alpha \psi^{-6} \hat{A}^{yy}) + \partial_z (2 \alpha \psi^{-6} \hat{A}^{yz}) - \frac{1}{3} \partial_y \Theta_\beta \approx

\approx \frac{\beta^y_{i-1,j,k}-2\beta^y_{i,j,k}+\beta^y_{i+1,j,k}}{h_x^2} + \frac{\beta^y_{i,j-1,k}-2\beta^y_{i,j,k}+\beta^y_{i,j+1,k}}{h_y^2} + \frac{\beta^y_{i,j,k-1}-2\beta^y_{i,j,k}+\beta^y_{i,j,k+1}}{h_z^2} =

= (\partial_x (2 \alpha \psi^{-6} \hat{A}^{yx}))_{i,j,k} + (\partial_y (2 \alpha \psi^{-6} \hat{A}^{yy}))_{i,j,k} + (\partial_z (2 \alpha \psi^{-6} \hat{A}^{yz}) )_{i,j,k} -

- \frac{1}{3} (\partial_y \Theta_\beta)_{i,j,k}.

Y, por último:

\Delta \beta^z = \partial_x (2 \alpha \psi^{-6} \hat{A}^{zx}) + \partial_y (2 \alpha \psi^{-6} \hat{A}^{zy}) + \partial_z (2 \alpha \psi^{-6} \hat{A}^{zz}) - \frac{1}{3} \partial_z \Theta_\beta \approx

\approx \frac{\beta^z_{i-1,j,k}-2\beta^z_{i,j,k}+\beta^z_{i+1,j,k}}{h_x^2} + \frac{\beta^z_{i,j-1,k}-2\beta^z_{i,j,k}+\beta^z_{i,j+1,k}}{h_y^2} + \frac{\beta^z_{i,j,k-1}-2\beta^z_{i,j,k}+\beta^z_{i,j,k+1}}{h_z^2} =

= (\partial_x (2 \alpha \psi^{-6} \hat{A}^{zx}))_{i,j,k} + (\partial_y (2 \alpha \psi^{-6} \hat{A}^{zy}))_{i,j,k} + (\partial_z (2 \alpha \psi^{-6} \hat{A}^{zz}) )_{i,j,k} -

- \frac{1}{3} (\partial_z \Theta_\beta)_{i,j,k}.

Parece que, del sistema no lineal acoplado inicial, hemos llegado a un sistema de diez ecuaciones desacopladas donde ocho de ellas son lineales y solo dos son no linales. No pinta mal. Ya escribiremos próximamente sobre las condiciones de contorno…

Vamos a discretizar las ecuaciones que comentamos en este post. Para ello, discretizaremos las derivadas de la siguiente manera:

\partial_x u = \frac{u_{i+1,j,k}-u_{i-1,j,k}}{2h_x},

\partial_y u = \frac{u_{i,j+1,k}-u_{i,j-1,k}}{2h_y},

\partial_z u = \frac{u_{i,j,k+1}-u_{i,j,k-1}}{2h_z},

\partial_{xx} u = \frac{u_{i-1,j,k}-2u_{i,j,k}+u_{i+1,j,k}}{h_x^2},

\partial_{yy} u = \frac{u_{i,j-1,k}-2u_{i,j,k}+u_{i,j+1,k}}{h_y^2},

\partial_{zz} u = \frac{u_{i,j,k-1}-2u_{i,j,k}+u_{i,j,k+1}}{h_z^2},

\partial_{xy} u = \frac{u_{i-1,j-1,k}-u_{i+1,j-1,k}-u_{i-1,j+1,k}+u_{i+1,j+1,k}}{4h_xh_y},

\partial_{xz} u = \frac{u_{i-1,j,k-1}-u_{i+1,j,k-1}-u_{i-1,j,k+1}+u_{i+1,j,k+1}}{4h_xh_z},

\partial_{yz} u = \frac{u_{i,j-1,k-1}-u_{i,j+1,k-1}-u_{i,j-1,k+1}+u_{i,j+1,k+1}}{4h_yh_z}.

El primer grupo de ecuaciones quedaría:

\partial_{xx} X^x + \partial_{yy} X^x + \partial_{zz} X^x = 8 \pi \psi^6 \rho h w^2 v_x - \frac{1}{3} \partial_x (\partial_x X^x + \partial_y X^y + \partial_z X^z) \approx

\approx \frac{X^x_{i-1,j,k}-2X^x_{i,j,k}+X^x_{i+1,j,k}}{h_x^2} + \frac{X^x_{i,j-1,k}-2X^x_{i,j,k}+X^x_{i,j+1,k}}{h_y^2} + \frac{X^x_{i,j,k-1}-2X^x_{i,j,k}+X^x_{i,j,k+1}}{h_z^2} =

= 8 \pi \psi^6_{i,j,k} \rho_{i,j,k} h_{i,j,k} w^2_{i,j,k} v_{x_{i,j,k}} - \frac{1}{3} ( \frac{X^x_{i-1,j,k}-2X^x_{i,j,k}+X^x_{i+1,j,k}}{h_x^2} +

+ \frac{X^y_{i-1,j-1,k}-X^y_{i+1,j-1,k}-X^y_{i-1,j+1,k}+X^y_{i+1,j+1,k}}{4h_xh_y} +

+ \frac{X^z_{i-1,j,k-1}-X^z_{i+1,j,k-1}-X^z_{i-1,j,k+1}+X^z_{i+1,j,k+1}}{4h_xh_z} ),

y además, para los esquemas de relajación no lineales, reescribimos la igualdad anterior como F(X^x_{i,j,k})=0 y entonces tenemos:

\partial_{X^x_{i,j,k}} F(X^x_{i,j,k}) = -2 ( \frac{4}{3}\frac{1}{h_x^2} + \frac{1}{h_y^2} + \frac{1}{h_z^2}).

\partial_{xx} X^y + \partial_{yy} X^y + \partial_{zz} X^y = 8 \pi \psi^6 \rho h w^2 v_y - \frac{1}{3} \partial_y (\partial_x X^x + \partial_y X^y + \partial_z X^z) \approx

\approx \frac{X^y_{i-1,j,k}-2X^y_{i,j,k}+X^y_{i+1,j,k}}{h_x^2} + \frac{X^y_{i,j-1,k}-2X^y_{i,j,k}+X^y_{i,j+1,k}}{h_y^2} + \frac{X^y_{i,j,k-1}-2X^y_{i,j,k}+X^y_{i,j,k+1}}{h_z^2} =

= 8 \pi \psi^6_{i,j,k} \rho_{i,j,k} h_{i,j,k} w^2_{i,j,k} v_{y_{i,j,k}} - \frac{1}{3} ( \frac{X^x_{i-1,j-1,k}-X^x_{i+1,j-1,k}-X^x_{i-1,j+1,k}+X^x_{i+1,j+1,k}}{4h_xh_y} +

+ \frac{X^y_{i,j-1,k}-2X^y_{i,j,k}+X^y_{i,j+1,k}}{h_y^2} +

+ \frac{X^z_{i-1,j,k-1}-X^z_{i+1,j,k-1}-X^z_{i-1,j,k+1}+X^z_{i+1,j,k+1}}{4h_yh_z} ),

con:

\partial_{X^y_{i,j,k}} F(X^y_{i,j,k}) = -2 ( \frac{1}{h_x^2} +\frac{4}{3} \frac{1}{h_y^2} + \frac{1}{h_z^2}).

\partial_{xx} X^z + \partial_{yy} X^z + \partial_{zz} X^z = 8 \pi \psi^6 \rho h w^2 v_z - \frac{1}{3} \partial_z (\partial_x X^x + \partial_y X^y + \partial_z X^z) \approx

\approx \frac{X^z_{i-1,j,k}-2X^z_{i,j,k}+X^z_{i+1,j,k}}{h_x^2} + \frac{X^z_{i,j-1,k}-2X^z_{i,j,k}+X^z_{i,j+1,k}}{h_y^2} + \frac{X^z_{i,j,k-1}-2X^z_{i,j,k}+X^z_{i,j,k+1}}{h_z^2} =

= 8 \pi \psi^6_{i,j,k} \rho_{i,j,k} h_{i,j,k} w^2_{i,j,k} v_{z_{i,j,k}} - \frac{1}{3} ( \frac{X^x_{i-1,j,k-1}-X^x_{i+1,j,k-1}-X^x_{i-1,j,k+1}+X^x_{i+1,j,k+1}}{4h_xh_z} +

+ \frac{X^y_{i,j-1,k-1}-X^y_{i,j+1,k-1}-X^y_{i,j-1,k+1}+X^y_{i,j+1,k+1}}{4h_yh_z} )

+ \frac{X^z_{i,j,k-1}-2X^z_{i,j,k}+X^z_{i,j,k+1}}{h_z^2}

con:

\partial_{X^z_{i,j,k}} = F(X^z_{i,j,k}) = -2 ( \frac{1}{h_x^2} + \frac{1}{h_y^2} + \frac{4}{3} \frac{1}{h_z^2}).

A continuación, discretizamos las siguientes ecuaciones:

\hat{A}^{xx} = 2 \partial_x X^x - \frac{2}{3} (\partial_x X^x + \partial_y X^y + \partial_z X^z) \approx

\approx \frac{2}{3}\frac{X^x_{i+1,j,k}-X^x_{i-1,j,k}}{h_x} -\frac{1}{3} \frac{X^y_{i,j+1,k}-X^y_{i,j-1,k}}{h_y} - \frac{1}{3} \frac{X^z_{i,j,k+1}-X^z_{i,j,k-1}}{h_z}) = \hat{A}^{xx}_{i,j,k},

\hat{A}^{xy} = \hat{A}^{yx}= \partial_x X^y + \partial_y X^x \approx

\approx \frac{X^y_{i+1,j,k}-X^y_{i-1,j,k}}{2h_x} + \frac{X^x_{i,j+1,k}-X^x_{i,j-1,k}}{2h_y} = \hat{A}^{xy}_{i,j,k} = \hat{A}^{yx}_{i,j,k},

\hat{A}^{xz} = \hat{A}^{zx} = \partial_x X^z + \partial_z X^x \approx

\approx \frac{X^z_{i+1,j,k}-X^z_{i-1,j,k}}{2h_x} + \frac{X^x_{i,j,k+1}-X^x_{i,j,k-1}}{2h_z} = \hat{A}^{xz}_{i,j,k} = \hat{A}^{zx}_{i,j,k},

\hat{A}^{yy} = 2 \partial_y X^y - \frac{2}{3} (\partial_x X^x + \partial_y X^y + \partial_z X^z) \approx

\approx -\frac{1}{3}\frac{X^x_{i+1,j,k}-X^x_{i-1,j,k}}{h_x} +\frac{2}{3} \frac{X^y_{i,j+1,k}-X^y_{i,j-1,k}}{h_y} - \frac{1}{3} \frac{X^z_{i,j,k+1}-X^z_{i,j,k-1}}{h_z}) = \hat{A}^{yy}_{i,j,k},

\hat{A}^{yz} = \hat{A}^{zy} = \partial_y X^z + \partial_z X^y \approx

\approx \frac{X^z_{i,j+1,k}-X^z_{i,j-1,k}}{2h_y} + \frac{X^y_{i,j,k+1}-X^y_{i,j,k-1}}{2h_z} = \hat{A}^{yz}_{i,j,k} = \hat{A}^{zy}_{i,j,k},

\hat{A}^{zz} = 2 \partial_z X^z - \frac{2}{3} (\partial_x X^x + \partial_y X^y + \partial_z X^z) \approx

\approx -\frac{1}{3}\frac{X^x_{i+1,j,k}-X^x_{i-1,j,k}}{h_x} -\frac{1}{3} \frac{X^y_{i,j+1,k}-X^y_{i,j-1,k}}{h_y} + \frac{2}{3} \frac{X^z_{i,j,k+1}-X^z_{i,j,k-1}}{h_z}) = \hat{A}^{zz}_{i,j,k}.

Por tanto, la siguiente ecuación:

\Delta \psi = -2 \pi \psi^{-1} (D + \tau) - \psi^{-7} \frac{(\hat{A}^{xx})^2+(\hat{A}^{yy})^2+(\hat{A}^{zz})^2+2(\hat{A}^{xy})^2+2(\hat{A}^{xz})^2+2(\hat{A}^{yz})^2}{8}

queda:

\approx \frac{\psi_{i-1,j,k}-2\psi_{i,j,k}+\psi_{i+1,j,k}}{h_x^2} + \frac{\psi_{i,j-1,k}-2\psi_{i,j,k}+\psi_{i,j+1,k}}{h_y^2} + \frac{\psi_{i,j,k-1}-2\psi_{i,j,k}+\psi_{i,j,k+1}}{h_z^2} =

=-2 \pi \psi^{-1}_{i,j,k} (D_{i,j,k}+\tau_{i,j,k}) -

- \frac{\psi^{-7}_{i,j,k}}{8} ( (\hat{A}^{xx}_{i,j,k})^2+(\hat{A}^{yy}_{i,j,k})^2+(\hat{A}^{zz}_{i,j,k})^2+2(\hat{A}^{xy}_{i,j,k})^2+2(\hat{A}^{xz}_{i,j,k})^2+2(\hat{A}^{yz}_{i,j,k})^2 ) ,

con:

\partial_{\psi_{i,j,k}} F(\psi_{i,j,k}) = -2 ( \frac{1}{h_x^2} + \frac{1}{h_y^2} + \frac{1}{h_z^2} ) -2 \pi \psi_{i,j,k}^{-2} (D_{i,j,k}+\tau_{i,j,k}) -

- \frac{7}{8} \psi^{-8}_{i,j,k} ( (\hat{A}^{xx}_{i,j,k})^2+(\hat{A}^{yy}_{i,j,k})^2+(\hat{A}^{zz}_{i,j,k})^2+2(\hat{A}^{xy}_{i,j,k})^2+2(\hat{A}^{xz}_{i,j,k})^2+2(\hat{A}^{yz}_{i,j,k})^2 ).

y la ecuación:

\Delta (\alpha\psi) = 2 \pi (\alpha\psi)^{-1} ( D + \tau + 2 \rho h (w^2-1) + 6 p) +

+ \frac{7}{8} (\alpha \psi)^{-7} ((\hat{A}^{xx})^2+(\hat{A}^{yy})^2+(\hat{A}^{zz})^2+2(\hat{A}^{xy})^2+2(\hat{A}^{xz})^2+2(\hat{A}^{yz})^2)

como:

\approx \frac{(\alpha\psi)_{i-1,j,k} - 2(\alpha\psi)_{i,j,k}+(\alpha\psi)_{i+1,j,k}}{h_x^2} + \frac{(\alpha\psi)_{i,j-1,k}-2(\alpha\psi)_{i,j,k}+(\alpha\psi)_{i,j+1,k}}{h_y^2} + \frac{(\alpha\psi)_{i,j,k-1}-2(\alpha\psi)_{i,j,k}+(\alpha\psi)_{i,j,k+1}}{h_z^2} =

=2 \pi (\alpha\psi)_{i,j,k}^{-1} (D_{i,j,k}+\tau_{i,j,k} + 2 \rho_{i,j,k} h_{i,j,k} (w^2_{i,j,k}-1)+6p_{i,j,k}) +

+ \frac{7}{8}(\alpha\psi)_{i,j,k}^{-7} ( (\hat{A}^{xx}_{i,j,k})^2+(\hat{A}^{yy}_{i,j,k})^2+(\hat{A}^{zz}_{i,j,k})^2+2(\hat{A}^{xy}_{i,j,k})^2+2(\hat{A}^{xz}_{i,j,k})^2+2(\hat{A}^{yz}_{i,j,k})^2 ) ,

donde:

\partial_{\psi\alpha_{i,j,k}} F(\psi\alpha_{i,j,k}) = -2 ( \frac{1}{h_x^2} + \frac{1}{h_y^2} + \frac{1}{h_z^2} ) +

+ 2 \pi (\psi\alpha)_{i,j,k}^{-2} (D_{i,j,k}+\tau_{i,j,k} + 2 \rho_{i,j,k} h_{i,j,k} (w^2_{i,j,k}-1)+6p_{i,j,k}) -

+ \frac{49}{8} (\psi\alpha)_{i,j,k}^{-8} ( (\hat{A}^{xx}_{i,j,k})^2+(\hat{A}^{yy}_{i,j,k})^2+(\hat{A}^{zz}_{i,j,k})^2+2(\hat{A}^{xy}_{i,j,k})^2+2(\hat{A}^{xz}_{i,j,k})^2+2(\hat{A}^{yz}_{i,j,k})^2 ).

Finalmente, tenemos:

\Delta \beta^x = \partial_x (2 \alpha \psi^{-6} \hat{A}^{xx}) + \partial_y (2 \alpha \psi^{-6} \hat{A}^{xy}) + \partial_z (2 \alpha \psi^{-6} \hat{A}^{xz}) -

- \frac{1}{3} \partial_x (\partial_x \beta^x + \partial_y \beta^y + \partial_z \beta^z) \approx

\approx \frac{\beta^x_{i-1,j,k}-2\beta^x_{i,j,k}+\beta^x_{i+1,j,k}}{h_x^2} + \frac{\beta^x_{i,j-1,k}-2\beta^x_{i,j,k}+\beta^x_{i,j+1,k}}{h_y^2} + \frac{\beta^x_{i,j,k-1}-2\beta^x_{i,j,k}+\beta^x_{i,j,k+1}}{h_z^2} =

= \frac{(\alpha \psi)_{i+1,j,k}^{-6} \hat{A}_{i+1,j,k}^{xx} - (\alpha \psi)_{i-1,j,k}^{-6} \hat{A}_{i-1,j,k}^{xx}}{h_x} +

+ \frac{(\alpha \psi)_{i,j+1,k}^{-6} \hat{A}_{i,j+1,k}^{xy} - (\alpha \psi)_{i,j-1,k}^{-6} \hat{A}_{i,j-1,k}^{xy}}{h_y} +

+ \frac{(\alpha \psi)_{i,j,k+1}^{-6} \hat{A}_{i,j,k+1}^{xz} - (\alpha \psi)_{i,j,k-1}^{-6} \hat{A}_{i,j,k-1}^{xz}}{h_z} -

- \frac{1}{3} ( \frac{\beta^x_{i-1,j,k}-2\beta^x_{i,j,k}+\beta^x_{i+1,j,k}}{h_x^2} +

+ \frac{\beta^y_{i-1,j-1,k}-\beta^y_{i+1,j-1,k}-\beta^y_{i-1,j+1,k}+\beta^y_{i+1,j+1,k}}{4 h_x h_y} +

+ \frac{\beta^z_{i-1,j,k-1}-\beta^z_{i+1,j,k-1}-\beta^z_{i-1,j,k+1}+\beta^z_{i+1,j,k+1}}{4 h_x h_z} ,

con:

\partial_{\beta^x_{i,j,k}} F(\beta^x_{i,j,k}) = -2 ( \frac{4}{3}\frac{1}{h_x^2} + \frac{1}{h_y^2} + \frac{1}{h_z^2}),

\Delta \beta^y = \partial_x (2 \alpha \psi^{-6} \hat{A}^{yx}) + \partial_y (2 \alpha \psi^{-6} \hat{A}^{yy}) + \partial_z (2 \alpha \psi^{-6} \hat{A}^{yz}) -

- \frac{1}{3} \partial_y (\partial_x \beta^x + \partial_y \beta^y + \partial_z \beta^z) \approx

\approx \frac{\beta^y_{i-1,j,k}-2\beta^y_{i,j,k}+\beta^y_{i+1,j,k}}{h_x^2} + \frac{\beta^y_{i,j-1,k}-2\beta^y_{i,j,k}+\beta^y_{i,j+1,k}}{h_y^2} + \frac{\beta^y_{i,j,k-1}-2\beta^y_{i,j,k}+\beta^y_{i,j,k+1}}{h_z^2} =

= \frac{(\alpha \psi)_{i+1,j,k}^{-6} \hat{A}_{i+1,j,k}^{yx} - (\alpha \psi)_{i-1,j,k}^{-6} \hat{A}_{i-1,j,k}^{yx}}{h_x} +

+ \frac{(\alpha \psi)_{i,j+1,k}^{-6} \hat{A}_{i,j+1,k}^{yy} - (\alpha \psi)_{i,j-1,k}^{-6} \hat{A}_{i,j-1,k}^{yy}}{h_y} +

+ \frac{(\alpha \psi)_{i,j,k+1}^{-6} \hat{A}_{i,j,k+1}^{yz} - (\alpha \psi)_{i,j,k-1}^{-6} \hat{A}_{i,j,k-1}^{yz}}{h_z} -

- \frac{1}{3} ( \frac{\beta^x_{i-1,j-1,k}-\beta^x_{i+1,j-1,k}-\beta^x_{i-1,j+1,k}+\beta^x_{i+1,j+1,k}}{4h_xh_y} +

+ \frac{\beta^y_{i,j-1,k}-2\beta^y_{i,j,k}+\beta^y_{i,j+1,k}}{h_y^2} +

+ \frac{\beta^z_{i-1,j,k-1}-\beta^z_{i+1,j,k-1}-\beta^z_{i-1,j,k+1}+\beta^z_{i+1,j,k+1}}{4h_yh_z} ),

con:

\partial_{\beta^y_{i,j,k}} F(\beta^y_{i,j,k}) = -2 ( \frac{1}{h_x^2} + \frac{4}{3} \frac{1}{h_y^2} + \frac{1}{h_z^2}),

\Delta \beta^z = \partial_x (2 \alpha \psi^{-6} \hat{A}^{zx}) + \partial_y (2 \alpha \psi^{-6} \hat{A}^{zy}) + \partial_z (2 \alpha \psi^{-6} \hat{A}^{zz}) -

- \frac{1}{3} \partial_z (\partial_x \beta^x + \partial_y \beta^y + \partial_z \beta^z) \approx

\approx \frac{\beta^z_{i-1,j,k}-2\beta^z_{i,j,k}+\beta^z_{i+1,j,k}}{h_x^2} + \frac{\beta^z_{i,j-1,k}-2\beta^z_{i,j,k}+\beta^z_{i,j+1,k}}{h_y^2} + \frac{\beta^z_{i,j,k-1}-2\beta^z_{i,j,k}+\beta^z_{i,j,k+1}}{h_z^2} =

= \frac{(\alpha \psi)_{i+1,j,k}^{-6} \hat{A}_{i+1,j,k}^{zx} - (\alpha \psi)_{i-1,j,k}^{-6} \hat{A}_{i-1,j,k}^{zx}}{h_x} +

+ \frac{(\alpha \psi)_{i,j+1,k}^{-6} \hat{A}_{i,j+1,k}^{zy} - (\alpha \psi)_{i,j-1,k}^{-6} \hat{A}_{i,j-1,k}^{zy}}{h_y} +

+ \frac{(\alpha \psi)_{i,j,k+1}^{-6} \hat{A}_{i,j,k+1}^{zz} - (\alpha \psi)_{i,j,k-1}^{-6} \hat{A}_{i,j,k-1}^{zz}}{h_z} -

- \frac{1}{3} ( \frac{\beta^x_{i-1,j,k-1}-\beta^x_{i+1,j,k-1}-\beta^x_{i-1,j,k+1}+\beta^x_{i+1,j,k+1}}{4h_xh_z} +

+ \frac{\beta^y_{i,j-1,k-1}-\beta^y_{i,j+1,k-1}-\beta^y_{i,j-1,k+1}+\beta^y_{i,j+1,k+1}}{4h_yh_z} )

+ \frac{\beta^z_{i,j,k-1}-2\beta^z_{i,j,k}+\beta^z_{i,j,k+1}}{h_z^2},

con:

\partial_{\beta^z_{i,j,k}} F(\beta^z_{i,j,k}) = -2 ( \frac{1}{h_x^2} + \frac{1}{h_y^2} + \frac{4}{3} \frac{1}{h_z^2} ).

Ya escribimos al respecto en este post. Aquí lo que haremos es reescribir las expresiones allí introducidas

En primer lugar, teniamos:

 \Delta X^i = 8 \pi f^{ij}S_j^* - \frac{1}{3}\mathcal{D}^i \mathcal{D}_j X^j

donde:

S_j^* := \sqrt{ \frac{\gamma}{f} } S = \psi^6 S_j,

S_j := \rho h w^2 v_j.

En el caso de estar trabajando en cartesianas y teniendo en cuenta todo el trabajo realizado en el artículo, nos queda:

\partial_{xx} X^x + \partial_{yy} X^x + \partial_{zz} X^x = 8 \pi \psi^6 \rho h w^2 v_x - \frac{1}{3} \partial_x (\partial_x X^x + \partial_y X^y + \partial_z X^z),

\partial_{xx} X^y + \partial_{yy} X^y + \partial_{zz} X^y = 8 \pi \psi^6 \rho h w^2 v_y - \frac{1}{3} \partial_y (\partial_x X^x + \partial_y X^y + \partial_z X^z),

\partial_{xx} X^z + \partial_{yy} X^z + \partial_{zz} X^z = 8 \pi \psi^6 \rho h w^2 v_z - \frac{1}{3} \partial_z (\partial_x X^x + \partial_y X^y + \partial_z X^z).

A continuación, y para la siguiente ecuación, necesitamos:

\hat{A}^{ij} = \mathcal{D}^i X^j + \mathcal{D}^j X^i - \frac{2}{3} \mathcal{D}_k X^k f^{ij}

que queda como:

\hat{A}^{xx} = 2 \partial_x X^x - \frac{2}{3} (\partial_x X^x + \partial_y X^y + \partial_z X^z),

\hat{A}^{xy} = \hat{A}^{yx}= \partial_x X^y + \partial_y X^x,

\hat{A}^{xz} = \hat{A}^{zx} = \partial_x X^z + \partial_z X^x,

\hat{A}^{yy} = 2 \partial_y X^y - \frac{2}{3} (\partial_x X^x + \partial_y X^y + \partial_z X^z),

\hat{A}^{yz} = \hat{A}^{zy} = \partial_y X^z + \partial_z X^y,

\hat{A}^{zz} = 2 \partial_z X^z - \frac{2}{3} (\partial_x X^x + \partial_y X^y + \partial_z X^z),

por lo que:

\Delta \psi = -2 \pi \psi^{-1} E^* - \psi^{-7} \frac{f_{il}f_{jm}\hat{A}^{lm}\hat{A}^{ij}}{8}

donde:

E^*:= \sqrt{ \frac{\gamma}{f} } E = \psi^6 E,

E:= D + \tau

es:

\Delta \psi = -2 \pi \psi^{-1} (D + \tau) - \psi^{-7} \frac{(\hat{A}^{xx})^2+(\hat{A}^{yy})^2+(\hat{A}^{zz})^2+2(\hat{A}^{xy})^2+2(\hat{A}^{xz})^2+2(\hat{A}^{yz})^2}{8}.

La siguiente:

\Delta (\alpha\psi) = 2 \pi (\alpha\psi)^{-1} (E^* + 2S^*) + \frac{7}{8} (\alpha\psi)^{-7} (f_{il} f{jm} \hat{A}^{lm} \hat{A}^{ij})

con:

S^*:= \sqrt{ \frac{\gamma}{f} } S = \psi^6 S,

S:= \rho h (w^2-1) + 3 p

queda:

\Delta (\alpha\psi) = 2 \pi (\alpha\psi)^{-1} ( D + \tau + 2 \rho h (w^2-1) + 6 p) +

+ \frac{7}{8}(\alpha\psi)^{-7} ((\hat{A}^{xx})^2+(\hat{A}^{yy})^2+(\hat{A}^{zz})^2+2(\hat{A}^{xy})^2+2(\hat{A}^{xz})^2+2(\hat{A}^{yz})^2)

Y la última:

\Delta \beta^i = \mathcal{D}_j (2 (\alpha\psi)^{-6} \hat{A}^{ij}) - \frac{1}{3} \mathcal{D}^i (\mathcal{D}_j \beta^j),

que escribimos como:

\Delta \beta^x = \partial_x (2 (\alpha \psi)^{-6} \hat{A}^{xx}) + \partial_y (2 (\alpha \psi)^{-6} \hat{A}^{xy}) + \partial_z (2 (\alpha \psi)^{-6} \hat{A}^{xz}) -

- \frac{1}{3} \partial_x (\partial_x \beta^x + \partial_y \beta^y + \partial_z \beta^z)

\Delta \beta^y = \partial_x (2 (\alpha \psi)^{-6} \hat{A}^{yx}) + \partial_y (2 (\alpha \psi)^{-6} \hat{A}^{yy}) + \partial_z (2 (\alpha \psi)^{-6} \hat{A}^{yz}) -

- \frac{1}{3} \partial_y (\partial_x \beta^x + \partial_y \beta^y + \partial_z \beta^z)

\Delta \beta^z = \partial_x (2 (\alpha \psi)^{-6} \hat{A}^{zx}) + \partial_y (2 (\alpha \psi)^{-6} \hat{A}^{zy}) + \partial_z (2 (\alpha \psi)^{-6} \hat{A}^{zz}) -

- \frac{1}{3} \partial_z (\partial_x \beta^x + \partial_y \beta^y + \partial_z \beta^z)

CoCoNuT es un código que permite realizar simulaciones de colapso estelar. Reescribimos las ecuaciones CFC, que son un caso particular de la aproximación FCF haciendo que las h^{ij} sean cero, en terminos de las variables que éste utiliza. Empezamos con una auxilar:

 \Delta X^i = 8 \pi f^{ij}S_j^* - \frac{1}{3}\mathcal{D}^i \mathcal{D}_j X^j

donde:

S_j^* := \sqrt{ \frac{\gamma}{f} } S = \psi^6 S_j,

S_j := \rho h w^2 v_j.

La primera es:

\Delta \psi = -2 \pi \psi^{-1} E^* - \psi^{-7} \frac{f_{il}f_{jm}\hat{A}^{lm}\hat{A}^{ij}}{8}

donde:

E^*:= \sqrt{ \frac{\gamma}{f} } E = \psi^6 E,

E:= D + \tau

La siguiente:

\Delta (\psi \alpha) = 2 \pi \alpha (E^* + 2S^*) + \alpha \psi^{-7} \frac{7 f_{il} f{jm} \hat{A}^{lm} \hat{A}^{ij}}{8}

con:

S^*:= \sqrt{ \frac{\gamma}{f} } S = \psi^6 S,

S:= \rho h (w^2-1) + 3 p

Y la última:

\Delta \beta^i = \mathcal{D}_j (2 \alpha \psi^{-6} \hat{A}^{ij}) - \frac{1}{3} \mathcal{D}^i (\mathcal{D}_j \beta^j).

Además, en CFC, tenemos:

\hat{A}^{ij} = (LX)^{ij} + \hat{A}^{ij}_{TT} \approx (LX)^{ij} = \mathcal{D}^i X^j + \mathcal{D}^j X^i - \frac{2}{3} \mathcal{D}_k X^k f^{ij}

donde L es el operador de Killing conforme actuando sobre la parte longitudinal X^i sin traza y A^{ij}_{TT} es la parte transversal sin traza de la curvatura extrínseca , y de FCF tenemos:

  • la métrica inducida en cada hipersuperficie \gamma_{\mu \nu} := g_{\mu \nu} + n_{\mu} n_{\nu} (o \boldsymbol{\gamma} := \boldsymbol{g} + \boldsymbol{n} \otimes \boldsymbol{n} ) con \boldsymbol{n} = \frac{dt}{|dt|}.
  • la curvatura extrínseca \boldsymbol{K:=-\frac{1}{2}\mathcal{L}_{\boldsymbol{n}} \boldsymbol{\gamma}} (o, con índices, K_{\mu \nu} = -\frac{1}{2} \mathcal{L}_{\boldsymbol{n}} \gamma_{\mu \nu}).

We try to write the formulas in its most general form and fix

First of all, in the simplest form of the HD equations we have

\frac{d\boldsymbol{v}_a}{dt} = -\sum_b m_b (\frac{P_a}{\rho_a^2} + \frac{P_b}{\rho_b^2}) \nabla_a W_{ab} (momentum)

where:

\rho_a = \sum_b m_b W_{ab} (density),

\frac{de_a}{dt} = \frac{1}{2} \sum_b m_b (\frac{P_a}{\rho_a^2} + \frac{P_b}{\rho_b^2}) \boldsymbol{v}_{ab} \cdot \nabla_a W_{ab} (density energy),

P_a and P_b

are obtained from the EoS and \frac{d\boldsymbol{r}_a}{dt} = \boldsymbol{v}_a.

In SHD we have:

d

Finally, the GHD equations are:

d

En el artículo A modified SPH approach for fluids with large density differences de F. Ott y E. Schnetter se explica una nueva modificación del método SPH que permite la interacción entre fluidos con densidades muy diferentes.

En la aproximación del SPH estandar tenemos:

\rho_i = \sum_j m_j W_{ij}

con W_{ij}:=W(\boldsymbol{x}_i - \boldsymbol{x}_j), o también:

\frac{d}{dt}\rho_i = \sum_j m_j \boldsymbol{v}_{ij} \cdot \nabla W_{ij}

donde \boldsymbol{v}_{ij}:=(\boldsymbol{v}_i - \boldsymbol{v}_j) y \nabla W_{ij}:=(\nabla W)(\boldsymbol{x}_i - \boldsymbol{x}_j), que tiene la ventaja de que los valores iniciales para \rho_i pueden escogerse libremente permitiendo la modelización de discontinuidades.

Ellos proponen utilizar,

\frac{d}{dt}\rho_i = m_i \sum_j \boldsymbol{v}_{ij} \cdot \nabla W_{ij}

deducida suavizando la densidad de partículas n_i = \frac{1}{V_i} en lugar de la densidad de masa \rho_i, donde V_i = \frac{m_i}{\rho_i} y \frac{d}{dt}m_i = 0.

justificandolo con tests positivos en tres problemas: ecuación de advección, onda sonora y tubo de choque.

En el review que Rosswog hace sobre el método SPH, en especial en sus aplicaciones a la astrofísica, hay un apartado dedicado al tratamiento de los shocks. El tratamiento de los strong shocks es uno de los puntos débiles del método SPH (penetración de partículas parcialmente resuelto mediante XSPH). En este apartado comenta que, básicamente, exiten dos maneras de tratarlos:

  1. Utilizar soluciones (analíticas) de un problema de Riemann entre dos partículas adyacentes.
  2. Ampliar la discontinuidad hasta una longitud resoluble numéricamente, para lo que se añaden terminos disipativos artificiales extra al flujo.

Comenta que la mayoria de las implementaciones de SPH utilizan esta segunda aproximación, consiguiendo esta amplicación o mediante el método de discretización utilizado, utilizando métodos que modifican las ecuaciones originales de una manera que nos puede venir bien, tipo el esquema Lax; o mediante la adición de terminos ad-hoc, algo bastante habitual, de tipo presión.

Pero también existen métodos que implementan la primera, y cita dos artículos:

  1. El artículo Reformulation of Smoothed Particle Hydrodynamics with Riemann Solver de Inutsuka y Shu-Ichiro.
  2. El artículo Implementations and tests of Godunov-type particle hydrodynamics de Cha y Whitworth.

En el primero comenta como consiguen, tras algunos intentos fallidos,  introducir Riemann Solver en las ecuaciones de manera similar a como se consigue en los métodos basados en malla, manteniendo la formulación conservativa de las mismas. Para empezar, comentan que mientras el error en la mayoria de los métodos basados en malla proviene de la truncación de la serie de de Taylor, en SPH, si consideramos:

\langle f \rangle(\boldsymbol{x}) := \int f(\boldsymbol{x'})W(\boldsymbol{x'}-\boldsymbol{x},h)d\boldsymbol{x'}

como la convolución de la función f(x) con el kernel. Mediante Taylor, de f alrededor de x=x_i, podemos ver que:

\langle f\rangle(\boldsymbol{x}_i) - f(\boldsymbol{x}_i) = \frac{h^2_{ef}}{4} \nabla^2 f(\boldsymbol{x}) + O(h^4_{ef}).

donde

h^2_{ef}:=2 \int x^2 W(\boldsymbol{x},h)d\boldsymbol{x},

por lo que el error proviene de la convolución con la función simétrica W.

También se tiene:

\langle \frac{\partial f}{\partial x} \rangle (\boldsymbol{x}) = \int \frac{\partial f(\boldsymbol{x'})}{\partial x'} W(\boldsymbol{x'}-\boldsymbol{x},h)d\boldsymbol{x'}

= - \int f(\boldsymbol{x'}) \frac{\partial}{\partial x'} W(\boldsymbol{x'}-\boldsymbol{x},h)d\boldsymbol{x'}

= \frac{\partial}{\partial x} \int f(x')W(x'-x,h)dx'

De manera que \nabla f = \nabla \langle f \rangle. Si definimos:

\rho(\boldsymbol{x}):=\sum_j m_j W(\boldsymbol{x}-\boldsymbol{x}_j,h),

entonces:

\sum_j \frac{m_j}{\rho(\boldsymbol{x})}W(\boldsymbol{x}-\boldsymbol{x}_j,h) = 1,

\sum_j m_j \nabla \big [ \frac{W(\boldsymbol{x}-\boldsymbol{x}_j,h)}{\rho(\boldsymbol{x})} \big ]= 0

que son la clave del método, pues ahora se puede hacer:

f_i:= \langle f \rangle (\boldsymbol{x}_i) = \int f(\boldsymbol{x'})W(\boldsymbol{x'}-\boldsymbol{x}_i,h)d\boldsymbol{x'}

= \int \sum_j m_j \frac{f(\boldsymbol{x}')}{\rho(\boldsymbol{x}')}W(\boldsymbol{x}'-\boldsymbol{x}_i,h)W(\boldsymbol{x}'-\boldsymbol{x}_j,h)d\boldsymbol{x}'

\sum_j \Delta f_{i,j}

donde:

\Delta f_{i,j}:=\int m_j \frac{f(\boldsymbol{x}')}{\rho(\boldsymbol{x}')}W(\boldsymbol{x}'-\boldsymbol{x}_i,h)W(\boldsymbol{x}'-\boldsymbol{x}_j,h)d\boldsymbol{x}'.

En el método SPH estandar nos encontramos, por una parte, con:

\sum_j \frac{m_j}{\rho(\boldsymbol{x}_j)}W(\boldsymbol{x}-\boldsymbol{x}_j,h) \approx 1

y por otra, con:

\Delta f_{i,j} \approx m_j \frac{f(\boldsymbol{x}_j)}{\rho(\boldsymbol{x}_j)}W(\boldsymbol{x}_i-\boldsymbol{x}_j,h)

y que, evidentemente, son aproximaciones mas pobres a las anteriormente obtenidas.

A continuación derivan las ecuaciones para la posición:

\ddot{\boldsymbol{x}}:=\int \frac{d\boldsymbol{v}(\boldsymbol{x})}{dt}W(\boldsymbol{x}-\boldsymbol{x}_i,h)d\boldsymbol{x}

= - \int \frac{1}{\rho(\boldsymbol{x})} \nabla P(x) W(\boldsymbol{x}-\boldsymbol{x}_i,h)d\boldsymbol{x}

y para la energía

\dot{u}_i := \int \frac{du(\boldsymbol{x})}{dt} W(\boldsymbol{x}-\boldsymbol{x}_i,h)d\boldsymbol{x}

= \sum_j m_j \int \frac{P}{\rho^2}[\boldsymbol{v}-\dot{\boldsymbol{x}}_i] \cdot [\nabla W(\boldsymbol{x}-\boldsymbol{x}_i,h)W(\boldsymbol{x}-\boldsymbol{x}_j,h) - W(\boldsymbol{x}-\boldsymbol{x}_i,h) \nabla W(\boldsymbol{x}-\boldsymbol{x}_j,h) ] d\boldsymbol{x},

que se modifica a

$latex$

para tener en cuenta la disipación

En el segundo artículo define e implementa el método Godunov-type Particle Hydrodynamics (GPH). Básicamente GPH = SPH + Godunov upwind schema.

En el artículo “MAGMA: a three-dimensional, Lagrangian magnetohydrodynamic code for merger applications” de S. Rosswog y D. Price comentan como introducir campos magnéticos en SPH. Las ecuaciones ya discretizadas quedan:

Ecuación de densidad:

\rho = \sum_b m_b W(r-r_b,h)

Ecuación del momento (conh: “grad-h” term, mag: magnetic force term, g: self-gravity and gravitational softening term)

\frac{d}{dt}v_{a,MHD} = \frac{d}{dt} (v_{a,h}+v_{a,h,dis}+v_{a,g}+v_{a,mag}+v_{a,mag,dis})

donde

\frac{d}{dt}v_{a,h} = -\sum_b m_b \big ( \frac{P_a}{\Omega_a \rho_a^2} \nabla_a W_{ab}(h_a)+ \frac{P_b}{\Omega_b \rho_b^2} \nabla_a W_{ab}(h_b) \big )

\frac{d}{dt}v_{a,h,dis} =

\frac{d}{dt}v_{a,g} = -G \sum_b m_b \big [ \frac{\phi'_{ab}(h_a) + \frac{\phi'_{ab}(h_b)}{}}{2} \big ]\hat{e}_{ab}

-\frac{G}{2} \sum_b m_b \big [ \frac{\zeta_a}{\Omega_a} \nabla_a W_{ab}(h_a) + \frac{\zeta_b}{\Omega_b} \nabla_a W_{ab}(h_b) \big ]

con

\phi'_{ab} = \frac{\partial \phi}{\partial|r_a-r_b|}, \zeta_k := \frac{\partial h_k}{\partial \rho_k} \sum_b m_b \frac{\partial \phi_{kb}(h_k)}{\partial h_k}

y

\Omega_a := \big ( 1- \frac{\partial h_a}{\partial \rho_a} \cdot \sum_b m_b \frac{\partial}{\partial h_a} W_{ab}(h_a) \big )

\frac{d}{dt}v_{a,mag} = - \sum_b \frac{m_b}{\mu_0} \big \{ \frac{B_a^2 / 2}{\Omega_a \rho_a^2} \nabla_a W_{ab}(h_a) + \frac{B_b^2 / 2}{\Omega_b \rho_b^2} \nabla_a W_{ab}(h_b) \}

+ \sum_b \frac{m_b}{\mu_0} \big \{ \frac{B_a(B_a \cdot \overline{\nabla_a W_{ab}}) - B_b(B_b \cdot \overline{\nabla_a W_{ab}})}{\rho_a \rho_b} \big \}

con

\overline{\nabla_a W_{ab}} = \frac{1}{2} \big [ \frac{1}{\Omega_a} \nabla_a W_{ab}(h_a) + \frac{1}{\Omega_b} \nabla_a W_{ab}(hb) \big ]

\frac{d}{dt}v_{a,mag,dis} =

Ecuación de la energía (con h: “grad-h” term, AV: Artificial Viscosity term y C: Condutivity):

\frac{d}{dt}u_{a,MDH} = \frac{d}{dt} (du_{a,h} + du_{a,AV} + du_{a,C})

donde

\frac{d}{dt}u_{a,h} = \frac{1}{\Omega_a} \frac{P_a}{\rho_a^2} \sum_b m_b v_{ab} \cdot \nabla_a W_{ab}(h_a)

\frac{d}{dt}u_{a,AV} =

\frac{d}{dt}u_{a,C} =

En la tesis doctoral “Relativistic simulations of compact object mergers for nucleonic matter and strange quark matter” de A. Bauswein comenta la implementación numérica del modelo físico para simular colisiones de objetos compactas.

La implementación numérica consta de dos partes: una que resuelve las ecuaciones de Euler tridimensionales relativistas que nos dan la evolución hidrodinámica del sistema y la otra que nos da una solución de las ecuaciones de campo de Einstein que nos proporcionan el campo gravitatorio donde se mueve el fluido.

Como ya comentamos en este post, utilizando SPH las ecuaciones de la hidrodinámica se tranforman en un conjunto de ODEs que podemos resolver mediante el método Runge-Kutta de cuarto orden. Las derivadas de los potenciales de la métrica se evaluan sobre una malla superpuesta y se mapean en las partículas. Se utiliza un paso de tiempo adaptativo para satisfacer la condición de Courant-Friedrichs-Levi. Todas las cantidades se representan por sus valores en las partículas, que corresponde con una visión de la hidrodinámica Lagrangiana, y la evolución temporal de las cantidades hidrodinámicas se calcula en la partículas.

Se añade un esquema de viscosidad artificial dependiente del tiempo para manejar choques hidrodinámicos. Este modelo resuelve un problema de Riemann local dado por los dos estados de partículas vecinas y añadiendo la correspondiente contribución a las ecuaciones de conservación del momento y de la energía.

A diferencia del SPH Newtoniano, la gravedad no se puede implementar dentro del esquema de partículas. Debemos resolver las ecuaciones de campo de Einstein, para lo que se utiliza el formalismo 3+1 o ADM que consiste en foliar el espacio-tiempo en hipersuperfícies espaciales con coordenada de tiempo constante. En esta aproximación, las ecuaciones de campo de Einstein se dividen en un conjunto de ecuaciones de ligadura y un conjunto de ecuaciones de evolución formulando de esta manera un problema de Cauchy.

El problema de Riemann es un caso especial de problema de valor inicial (IVP) en el que la PDE es:

u_t + a u_x = 0

y la condicion inicial (IC):

u(x,0) = u_0(x) = \begin{cases} u_L\mbox{ si } x < 0 \\ u_R \mbox{ si } x > 0 \end{cases}

donde u_L y u_R son dos valores constantes, de manera que tenemos una discontinuidad en x=0. En este caso, la solución es sencilla:

u(x,t) = u_0(x-at) = \begin{cases} u_L \mbox{ si } x-at < 0\\ u_R \mbox{ si } x-at > 0 \end{cases}.

Podemos extender el problema a un conjunto de m PDEs hiperbólicas:

U_t + AU_x = 0 con -\infty < x < \infty, t>0

donde la matriz A es de coeficientes constantes. Al asumir la hiperbolicidad de A, tenemos m valores propios reales \lambda_i y m vectores propios independientes K^{(i)}. En este caso, la IC se escribe como:

U(x,0) = U^{(0)}(x) = \begin{cases} U_L, x<0\\ U_R, x>0 \end{cases}

Los Riemann solvers son métodos numéricos que permiten resolver el problema de Riemann.

En [Monaghan 2005] se revisan la teoría y las aplicaciones del SPH desde su aparición hasta la actualidad.

El artículo empieza comentando cinco de los puntos fuertes del método SPH:

  1. La advección se trata de manera exacta: si a una partícula se le asigna un color, y se le especifica una velocidad, el transporte del color por la partícula del sistema es exacto.
  2. Cuando trabajamos con mas de un material, cada uno se describe mediante su propio conjunto de partículas de manera que los problemas de interfaz resultan triviales.
  3. Cierra la brecha entre el contínuo y la fragmentación de una forma natural.
  4. La resolución puede depender tanto de la posición como del tiempo.
  5. SPH tiene la ventaja computacional de que la computación ocurre solo donde está realmente la materia, con la consecuente reducción de calculo y almacenamiento.

En [Monaghan 1992] comenta algunos detalles de implementación del método SPH.

La inicialización de un cálculo SPH comienza especificando la masa, posición, velocidad y energía térmica de cada partícula. Si calculamos la densidad mediante la forma alternativa:

\frac{d}{dt}\rho_a = \sum_b m_b v_{ab} \nabla_a W_{ab}

donde v_{ab} = v_a - v_b y que tiene algunas ventajas (evita la oscilación de los bordes del fluido y es computacionalmente mas eficiente el calculo de la tasa de cambio de todas las variables físicas), necesitamos asignar la densidad inicial de cada partícula. También la mezcla de elementos de cada partícula puede ser necesaria.

Para inicializar las partículas, es conveniente situarlas en una malla regular, que puede ser una malla Cartesiana regular con igual espaciado o una malla cúbica centrada en el cuerpo en tres dimensiones (mejor representación de las integrales por sumas, mejor detección de vecinos y mejor distribución de partículas en caso de compresión planar). Además, si el tamaño de celda asociado con una partícula a es \Delta V_a entonces m_a se toma como \rho \Delta V_a.

Si utilizamos la misma h para cada partícula, entonces basta trabajar con listas enlazadas. Si el kernel esta basado en splines, las celdas enlazadas por la lista deben ser de 2h de ancho de manera que solo las celdas vecinas pueden contribuir en las partículas en una celda dada.

Si cada partícula tiene su propia h entonces el calculo de las sumas del SPH pueden formar parte de un calculo en arbol que es el método natural para calcular las fuerzas autogravitantes de un conjunto de partículas. Este codigo en arbol puede ser vectorizado.

Trabajando en GRMHD nos encontramos el siguiente problema de Cauchy:

  1. Ecuaciones de Einstein: (\gamma,K)_{t=0} \longrightarrow (\gamma,K)_t
  2. Hidrodinámica: (\rho,q,\tau)_{t=0} \longrightarrow (\rho,q,\tau)_t
  3. Electromagetismo: (E,B)_{t=0} \longrightarrow (E,B)_t

donde \gamma es la métrica y K la curvatura extrinseca; \rho es la densidad de masa, q la densidad de momento, \tau la densidad de energia; E es el campo de fuerza eléctrica y B la inducción magnética.

En [Rosswog 2009] tenemos las ecuaciones de la hidrodinámica en forma Lagrangiana discretizadas y en su forma mas básica. Las partículas avanzaran en el tiempo siguiendo las siguientes ecuaciones:

Para empezar, no hay necesidad de resolver la ecuación de continuidad ya que la masa de las partículas permanece fija. Podemos obtener las densidades mediante:

\rho_a = \sum_b m_b W_{ab}

La ecuación del momento queda:

\frac{d}{dt}\vec{v}_a = - \sum_b m_b (\frac{P_a}{\rho_a^2} + \frac{P_b}{\rho_b^2} + \Pi_{ab} ) \nabla_a W_{ab}

La ecuación de evolución para la energía interna específica puede escribirse como:

\frac{d}{dt} u_a = \sum_b m_b (\frac{P_a}{\rho_a^2} + \frac{1}{2}\Pi_{ab}) \vec{v}_{ab} \nabla_b W_{ab}

Rosswog las llama “vanilla ice” SPH.

En las Jornadas sobre los problemas del milenio celebradas en Barcelona del 1 al 3 de junio de 2011, Diego Córdoba dió una charla sobre el problema Clay de las ecuaciones de Navier-Stokes para la que escribió estas notas.

Un fluido es ideal si es incompresible, homogéneo  y perfecto.

La idea del problema consiste en determinar si:

…un fluido incompresible con energía finita puede desarrollar singularidades en tiempo finitos.

Mas formalmente, si consideramos un fluido viscoso (\nu > 0), homogéneo (\rho = 1) e incompresible (\nabla \cdot u = 0), tenemos:

u_t + u \cdot \nabla u = - \nabla p + \nu \Delta u + f con \nu >0, x \in \mathbb{R}^3, t \geq 0

\nabla \cdot u = 0

u(x,0) = u_0

donde cada partícula en el tiempo t está en la posición x = (x_1,x_2,x_3) del dominio que ocupa el fluido \Omega \subset \mathbb{R}^3, u(x,t) = (u_1(x,t), u_2(x,t), u_3(x,t)) es el campo de velocidades, p = p(x,t) son las presiones en el seno del fluido y \rho = \rho(x,t) es la densidad. Además, \nu = cte \geq 0 es la viscosidad y f la fuerza externa.

La fuerza exterior  debe verificar:

|\partial_x^\alpha\partial_t^m f| \leq C_{\alpha,k,m}(1+|x|+t)^{-k} para todo \alpha, k, m > 0

y el dato inicial las siguientes condiciones de regularidad:

|\partial_x^\alpha u_{0i}| \leq C_{\alpha,k}(1+|x|)^{-k} para todo \alpha, k > 0

Las soluciones (u,p) admisibles para x \in \mathbb{R}^3 estan o en C^{\infty}(\mathbb{R}^3 \times [0,\infty)) con decaimiento en el infinito de la presión y de energía finita, es decir, \int_{\mathbb{R}^3} |u|^2 dx < \infty para todo t, o estan en C^{\infty}(\mathbb{T}^3 \times [0,\infty)) periódicas y la presión de media cero.

Sea u_0 satisfaciendo las condiciones de regularidad. El problema del Instituto Clay consiste, entonces, en determinar si siempre existen soluciones admisibles para u_0 o si existe algún caso en que no.

En [Monaghan 1992] comenta el caso del método SPH en relatividad especial.

Para empezar asumimos que el fluido está constituido por bariones, por lo que el tensor de energia-momento es:

T^{\mu \nu} = (n m_0 c^2 + n \tau + P) U^\mu U^\nu + P g^{\mu \nu}

donde los indices griegos van de 0 a 3 y los coeficientes de la métrica se definen

g_{00} = -1 y g_{ij} = 1

En estas ecuaciones, n representa la densidad de bariones, P es la presión, \tau la energía térmica, c la velocidad del sonido, U^\nu la 4-velocidad con U_\nu U^\nu = -1 y m_0 la masa.

Las ecuaciones del momento se siguen de:

\frac{\partial}{\partial x^\nu} T^{i \nu} = 0

que es:

\frac{d}{dt} q = - \frac{1}{N} \nabla P

y en forma SPH queda:

\frac{d}{dt}q_a = -\sum_b \nu_b (\frac{P_a}{n_a^2} + \frac{P_b}{N_b^2}) \nabla_a W_{ab}

donde \nu_b es el número de bariones asociados a la partícula b.

Y la de la energía se sigue de:

\frac{\partial}{\partial x^j} T^{0j} = 0

que es:

\frac{d}{dt} \epsilon = - \frac{1}{N} \nabla \cdot (Pv)

y en foma SPH queda:

\frac{d}{dt} \epsilon_a = -\sum_b m_b (\frac{P_a v_a}{N_a^2} + \frac{P_b v_b}{N_b^2}) \nabla W_{ab}

Según Monaghan en [Monaghan, 1992], uno de los padres del método, buscaban un método con el que fuera fácil trabajar y diera buenos resultados y encontraron eso y mucho mas.

Además de no  necesitar malla para calcular las derivadas, puesto que se calculan de manera analítica a partir de una fórmula de interpolación, las ecuaciones de conservación del momento y la energía pasan a ser un conjunto de ecuaciones diferenciales ordinarias fáciles de entender.

En esencia, el método SPH es un método de interpolación que permite expresar una función a partir de los valores en un conjunto desordenado de puntos a los que llamamos particulas.

La idea es que aproximamos cualquier función A(r) de la siguiente manera:

A_I(r) = \int A(r') W(r-r',h) dr'

donde r es el vector posición, W es una función de pesos o kernel (con las propiedades ya comentadas en otros posts) y h es la longitud de suavizado. Este tipo de aproximaciones reciben el nombre de interpolación integral.

Para trabajar de manera numérica, la interpolador integral se aproxima por:

A_S(r) = \sum_b m_b \frac{A_b}{\rho_b} W(r-r_b,h)

donde el índice del sumatorio b identifica a cada partícula y el sumatorio es sobre todas las partículas. La partícula b tiene masa m_b, posición r_b, densidad \rho_b y velocidad v_b. El valor de la función A en r_b es A_b.

En el artículo [Rosswog 2009], Stephan Rosswog hace un repaso detallado del método SPH centrandose especialmente en sus aplicaciones en astrofísica. Repasamos el apartado que hace referencia a las ecuaciones de la hidrodinámica en forma Lagrangiana.

A diferencia de los metodos basados en malla, que son Eulerianos, es decir, métodos donde  describimos el fluido desde un punto fijo del espacio, el Smoothed Particle Hydrodynamics es totalmente Lagrangiano, por lo que describimos el fluido desde un sistema de coordenadas fijado en una particula del fluido en movimiento.

La derivada Lagrangiana o sustancial respecto del tiempo, \frac{d}{dt}, se relaciona con la derivada Euleriana respecto al tiempo, \frac{\partial}{\partial t} de la siguiente manera:

\frac{d}{dt} = \frac{dx^i}{dt} \frac{\partial}{\partial x^i} + \frac{\partial}{\partial t} = \vec{v} \cdot \nabla + \frac{\partial}{\partial t}

Aplicando esta relación a las ecuaciones en forma Euleriana, las ecuaciones de la hidrodinámica en forma Lagrangiana quedan:

  1. Ecuación de continuidad: \frac{d}{dt} \rho = - \rho \nabla \cdot \vec{v}
  2. Ecuacion del momento: \frac{d}{dt} \vec{v} = -\frac{\nabla P}{\rho} + \vec{f}
  3. Ecuación de la energía: \frac{d}{dt}u = \frac{P}{\rho^2} \frac{d}{dt} \rho = - \frac{P}{\rho} \nabla \cdot \vec{v}
  4. Ecuación de estado, que describe la termodinámica del fluido estelar: P = (\gamma -1) \cdot \rho \cdot \epsilon (ecuación del gas ideal)

La hipótesis de Rieman dice que todos los ceros no triviales de la función zeta de Riemann:

\zeta(z) = \sum_{n=1}^{\infty} \frac{1}{n^z}

tienen parte real \frac{1}{2}.

En el entretenido libro “La música de los números primos” de Marcus de Sautoy, éste comenta la posibilidad de que la teoría de números y la física estén mas estrechamente relacionadas de lo que pensamos:

Los físicos creen que la razón por la que los ceros de Riemann deben situarse todos sobre la recta es que terminarán por ser las frecuencias de un tambor matemático. A un cero que se situara fuera de la recta le correspondería una frecuencia imaginaria prohibida por la teoría.

Y para afianzar esta idea, comenta un problema clásico de hidrodinámica resuelto por Bernhard Riemann argumentando de la misma manera:

El problema se refiere a una esfera de fluido en rotación que se mantiene unida gracias a interacciones gravitacionales recíprocas entre las partículas que la componen. Una estrella, por ejemplo, es una enorme bola de gas giratorio que se mantiene unido por su propia gravedad. La cuestión es: ¿qué sucederá con la bola si se le da una patada?¿Se limitará a temblar ligeramente o se desitegrará? Para responder a estas preguntas es necesario determinar si ciertos números imaginarios determinados están o no alineados. Si lo están, la esfera de fluido en rotación quedará intacta.

Los Nachlass son un conjunto de notas inéditas de Riemann que están en la biblioteca de Gotinga. Cuenta el libro que cuando el físico Jon Keating pidió dos partes de los Nachlass, una correspondiente a sus intentos de demostración de la hipotesis de Riemann y otra a sus estudios en hidrodinámica, el bibliotecario le entrego un único grupo de documentos. De ahí que Marcus escriba:

Una vez más, los Nachlass revelaban hasta qué punto Riemann se adelantó a su tiempo. Es imposible que no fuera consciente del significado que implicaba su solución al problema de dinámica de fluidos. Su método había demostrado por qué cietos números imaginarios que aparecían en su análisis de la esfera de fluido se colocaban en línea recta; y al mismo tiempo -y en los mismos folios- estaba intentando demostrar por qué los ceros de su paisaje zeta se situaban todos sobre la misma línea.

 

Existen diferentes posibilidades a la hora de definir una función kernel:

  1. Gaussiana [Gingold & Monaghan, 1977]:W(r,h) = \alpha_D \cdot e^{-q^2} con 0 \leq q \leq 2 donde q=\frac{r}{h}, r es la distancia entre dos partículas determinadas y \alpha_D, el factor dimensional, que es \frac{1}{\pi h^2} en dos dimensiones y \frac{1}{\pi^{\frac{3}{2}} h^3} en tres.
  2. Cuadrática [Gingold & Monaghan, 1977]: \alpha_D (\frac{3}{16} q^2 - \frac{3}{4} q + \frac{3}{4}) con 0 \leq q \leq 2 y donde \alpha_D es \frac{2}{\pi h^2} en 2D y \frac{5}{4 \pi h^3} en 3D.
  3. B-spline cúbico [Monaghan & Lattancio, 1985]W(r,h) = \alpha_D \begin{cases} 1 - \frac{3}{2} q^2 + \frac{3}{4} q^3 ,& \mbox{if } 0 \leq q \leq 1 \\ \frac{1}{4}(2-q)^3 ,& \mbox{if } 1 < q \leq 2 \\ 0 ,& \mbox{if } q > 2\end{cases} con \alpha_D = \frac{10}{7 \pi h^2} en dos dimensiones y \frac{1}{\pi h^3} en 3 dimensiones.
  4. Quíntica: W(r,h) = \alpha_D (1-\frac{q}{2}^4)(2q + 1) con 0 \leq q \leq 2 y \alpha_D = \frac{7}{4 \pi h^2} en 2D y \frac{7}{8 \pi h^3} en 3D.

Pensando en posibles implementaciones, parece lógico utilizar una classe abstracta Kernel de la que heredaran las anteriores implementando sus métodos, de manera que, las posibles clases que llamen a la misma puedan hacerlo siempre de la misma manera independientement de cual estemos utilizando.

En el artíclo [Gingold & Monaghan, 1977] se introducen

La función kernel, W(\vec{r}-\vec{r'},h), es una función que permite interpolar los valores de cualquier propiedad del fluido en función del valor de las partículas del entorno. Su papel es similar al de los diferentes esquemas de diferencias en el ámbito de las Diferencias Finitas o las funciones de forma en los Elementos Finitos.

Existen diferentes funciones kernel: Gaussiana, cuadrática, spline cúbico, quíntica, etc.

La función kernel debe cumplir:

  1. Positiva: W(r-r',h) \geq 0 dentro del dominio.
  2. Soporte compacto: W(r-r',h) = 0 fuera del dominio.
  3. Normalizada: \int W(r-r',h) dr' = 1.
  4. Comportamiento de función delta: \lim_{h \rightarrow 0} W(r-r',h) dr' = \delta(r-r').
  5. Monotona decreciente.

Para derivar, tomamos la derivada analítica de la suma aproximada. De esta manera, como la derivada de la función kernel es conocida, no necesitamos diferencias finitas y el conjunto de ecuaciones PDE pasa a ser ODE.

\nabla f(\vec{r}) = \sum_b \frac{m_b}{\rho_b} f_b W(\vec{r}-\vec{r'},h)

En el artículo  [Gingold & Monaghan, 1977] se presenta por primera vez el método Smoothed Particle Hydrodynamics. Los autores, originalmente, buscaban un método que permitiera tratar problemas en astrofísica asimétricos (sin simetría esférica, sin simetría axial, etc.) . En estos casos, los métodos de diferencias finitas no se adaptan bien, pues requieren elevar el número de puntos en la malla para seguir con la precisión deseada su evolución, lo cual complica enormemente la evolución de las integrales múltiples.

Lo que pensaron es utilizar la descripción Lagrangiana del flujo del fluido que centra su atención en los elementos del fluido. En la discretización, estos elementos se mueven siguiendo las leyes de Newton con fuerzas debidas a los gradientes de presión y a otras fuerzas de cuerpo como la gravedad, rotación o magnéticas. ¿Qué método utilizar para determinar las fuerzas que actuan en un momento determinado sobre un elemento de fluido?

Para empezar, para elementos de fluido de igual masa, el número de elementos por unidad de volumen debe ser proporcional a la densidad. Además, sin ningún tipo especial de simetría, la posicion de los elementos será aleatoria de acuerdo con la densidad. Para recuperar la densidad de la distribución conocida de los elementos es equivalente a recuperar una distribución de probabilidad a partir de una muestra. Existen dos métodos para conseguir esto que funcionan bien con problemas de fluidos: el smoothing kernel method y  la técnica del spline delta. Ambos métodos se pueden pensar como la aproximación de una integral por el procedimiento de Monte Carlo.

En relatividad especial, o relatividad restringida, podemos generalizar las ecuaciones de Euler clásicas que determinan la evolución de los fluidos.

Para ello, definimos el tensor de energia-impulso de la siguiente manera:

T^{\alpha \beta} = (\rho + \frac{p}{c^2}) u^\alpha u^\beta + p \eta^{\alpha \beta}

De esta manera:

\frac{\partial}{\partial x^\beta} T^{\alpha \beta} = 0

generaliza tanto la ecuación de conservación de la masa como las ecuaciones de Euler.

octubre 2017
L M X J V S D
« Ago    
 1
2345678
9101112131415
16171819202122
23242526272829
3031