How is the divergence and curl approximated in Smoothed Particle Hydrodynamics?

Question

I was reading through Monaghan, SPH, 1992 paper and he states that any quantity $A$ at $r_b$ can be approximated by $A(r) = \sum_b m_b \frac{A_b}{\rho_b}W(r-r_b,h)$ and so the gradient of that quantity can be computed as $\nabla A(r) = \sum_b m_b \frac{A_b}{\rho_b} \nabla W(r-r_b,h)$.

So if we want to approximate the density we have $\rho (r) = \sum_b m_b W(r-r_b,h)$, so far so good.

Then there is the equation to approximate the divergence of the velocity as $\nabla \cdot v = \sum_b m_b v_b \cdot \nabla W(r-r_b,h)$ which I don't understand exactly where it comes from. Is the velocity $v$ playing the role of the approximated quantity, $A_b = v$ ?

Now, then the

second golden rule of SPH is to rewrite formulae with density placed inside operators

which gives $\nabla \cdot v = [\nabla \cdot (\rho v) - v \cdot \nabla \rho] / \rho$. How is this obtained? why is the density appearing in the formula and how is related to the above formula for the divergence of the velocity?

I guess something similar happens when deriving the estimator for the curl of the velocity.

Can someone explain a little more in details how are the formulas derived?

Answer 1

When you approximate a quantity, for example $A(r) = \sum_b m_b \frac{A_b}{\rho_b}W(r-r_b,h)$, then what you are saying is that you are assigning a scalar value $A_b$ to every nodal point $b$. At any point $r$, this formula then simply provides a way to interpolate between the points $r_b$ in the neighborhood. Because the masses and nodal values $A_b$ are just fixed scalars, the gradient of your approximation is simply $\nabla A(r) = \sum_b m_b \frac{A_b}{\rho_b}\nabla W(r-r_b,h)$.

You can do the same construction with vector fields. If you assign a vector $v_b$ to every nodal point $b$, then your approximation becomes $\mathbf v(r) = \sum_b m_b \frac{\mathbf v_b}{\rho_b}W(r-r_b,h)$ and its divergence -- because the $\mathbf v_b$ are just fixed vectors, independent of $r$ -- is $\nabla \cdot \mathbf v(r) = \sum_b m_b \frac{\mathbf v_b}{\rho_b} \cdot \nabla W(r-r_b,h)$. The same can of course be done to the curl: $\nabla \times \mathbf v(r) = \sum_b m_b \nabla \times \left[\frac{\mathbf v_b}{\rho_b} W(r-r_b,h)\right] = \sum_b m_b \left[\frac{\mathbf v_b}{\rho_b} \times \nabla W(r-r_b,h)\right]$.

Answer 2

Your second question is answered by expanding the divergence of $\rho v$ and rearranging terms: $$ \nabla\cdot({\rho v}) = \nabla\rho\cdot v + \rho \nabla \cdot v$$ so $$ \nabla \cdot v = \frac{\nabla\cdot({\rho v}) - \nabla\rho\cdot v}{\rho} $$

Answer 3

Given the fundamental idea of the SPH interpolant with an arbitrary function $A$ (which could be the velocity $\textbf{v}$):

$A(r)=\int_{-\infty}^{\infty} A(r')\delta(r-r')dr'$,

the convolution of an arbitrary $A(r)$ with the Dirac delta $\delta$. Obviously, is not applicable in numerical simulations due to two reasons:

• $\delta$ has infinitesimally small radius, and infinite height
• The integral can be evaluated only analytically.

To solve both issues, first we change $\delta$ to a mollifier function (kernel) to achieve a finite (or sometimes infinitely large) influence radius:

$A(r)\approx\int_{\Omega} A(r')W(r-r',h)dr'$,

where $h$ is analogous to the standard deviation in case of Gaussian function and $\Omega$ is the domain of the kernel. $W(r-r',h)$ is usually a compact bell shaped function of unit volume. Secondly, change the integral to the summation form:

$\langle A(r_i)\rangle = \sum_{j}{A(r_j)\frac{m_j}{\rho_j}W(r_i-r_j,h)}$.

Now $r_i$ is corresponding to the $i^{th}$ discrete point often referred to as particle. This is can be considered as a discrete convolution.

Now apply the continuum convolution to the derivative of A:

$\nabla_{r} A(r)=\int_{\Omega} \nabla_{r'} A(r')W(r-r',h)dr'$,

where $\nabla_{r}$ denotes the derivative at $r$.

Using the rules of differentiation, it also can be written as:

$\nabla_{r} A(r)\approx\int_{\Omega} \nabla_{r'} \big[A(r')W(r-r',h)\big]dr'-\int_{\Omega}A(r')\nabla_{r'} W(r-r',h)dr'$.

Transform the first term in the RHS from volume integral to surface integral using the Gauss-Ostrogradsky theorem:

$\nabla A(r)\approx\int_{\partial\Omega}A(r')W(r-r',h)dS-\int_{\Omega}A(r')\nabla_{r'} W(r-r',h)dr'$.

Since $W(r-r',h)$ is compact and has a zero value at its boundary, the first term on RHS is now zero, so:

$\nabla A(r)\approx-\int_{\Omega}A(r')\nabla_{r'} W(r-r',h)dr'$,

where we can use the identity

$\nabla_{r'}W(r-r',h)=-\nabla_{r}W(r-r',h)$,

to obtain the discrete form

$\langle \nabla A_i\rangle = \sum_{j}{A_j\frac{m_j}{\rho_j}\nabla_i W(r_i-r_j,h)}$.

Unfortunately, this gives a poor estimation on the derivative, because it is extremely sensitive to the particle distribution. To overcome this problem, you can apply the discrete convolution to the formula you've shown:

$\langle\nabla A_i\rangle=\frac{\langle \nabla(\rho_i A_i)\rangle-A_i\langle\nabla \rho_i\rangle}{\rho_i} = \frac{1}{\rho_i}\sum_{j}{(\rho_jA_j)\frac{m_j}{\rho_j}\nabla_i W(r_i-r_j,h)}-\frac{1}{\rho_i}\sum_{j}{(\rho_jA_i)\frac{m_j}{\rho_j}\nabla_i W(r_i-r_j,h)}$,

which is finally:

$\langle\nabla A_i\rangle=\frac{1}{\rho_i}\sum_{j}{(A_j-A_i)m_j \nabla_i W(r_i-r_j,h)}$.

It is important, that this is only one of the several possible forms of differential operators in SPH. They have different consistency and conservation properties forming a trade-off. Which to use? It depends on the equation you want to solve.

Nevertheless, you are free to change the product to use the obtained operator for curl:

$\langle\nabla\times \textbf{v}_i\rangle=\frac{1}{\rho_i}\sum_{j}{(\textbf{v}_j-\textbf{v}_i)\times\nabla_i W(r_i-r_j,h)m_j}$