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.

Anuncios