Finished the viscosity and thermal conductivity of the theory document. The...

Finished the viscosity and thermal conductivity of the theory document. The equation set is now complete.


Former-commit-id: 34b780bc7102b62d441c2eb85c3b84aa77baece6

theory/latex/sph.tex

@@ -29,15 +29,26 @@ Every particle contains the following information:
  Internal energy per unit mass & Primary & $u$ & $[m^2 \cdot s^{-2}]$ \\ 
  Energy derivative & Tertiary & $\frac{du}{dt}$ & $[ m^2 \cdot s^{-3}]$ \\
  Correction term & Secondary & $\Omega$ & $[-]$ \\
- Smoothing length & --- &$h$ & $[m]$ \\
- Smoothing length derivative & --- &$\frac{dh}{dt}$ & $[m\cdot s^{-1}]$ \\
- Time step & --- & $\Delta t$ & $[s]$ \\
+ Smoothing length & Secondary &$h$ & $[m]$ \\
+ Smoothing length derivative & Tertiary &$\frac{dh}{dt}$ & $[m\cdot s^{-1}]$ \\
+ Time step & Secondary & $\Delta t$ & $[s]$ \\
+\hline
+ Artificial viscosity & Primary & $\alpha$ & $[-]$\\
+ Artificial conductivity & Primary & $\alpha_u$ & $[-]$\\
+ Velocity divergence & Secondary & $\nabla\cdot \vec{v}$ & $[s^{-1}]$ \\
+ Internal energy laplacian & Secondary & $\nabla u$ & $[m\cdot s^{-2}]$\\
 \hline
 \end{tabular} 
 \end{table}
 
+Secondary quantities are computed from the primary one in a loop of particle neighbors. Tertiary ones are computed from
+seondary ones in another loop. \\
+
 For optimization purposes, any function of these quantities could be stored. For instance, $1/h$ instead of $h$ or
-$\frac{P}{\rho\Omega}$ instead of $\Omega$ may be options worth exploring.
+$\frac{P}{\rho\Omega}$ instead of $\Omega$ may be options worth exploring. \\
+
+The four quantities in the second half of the table are used in improved state-of-the-art implementations of SPH. In a
+first approximations, they can be neglected.
 
 \section{Kernel function}
 
@@ -131,7 +142,8 @@ physically relevant parameter $\eta$. In 3D this reads
 \end{equation}
 
 We usually use $N_{ngb} = 48$, which corresponds to a sub-optimal value of $\eta=1.127$. The optimal value should be
-$N_{ngb}=57.9$ but this is computationally more expensive and the improvement over $N_{ngb}=48$ is not obvious.\\ 
+$N_{ngb}=57.9$ ($\eta=1.2$) but this is computationally more expensive and the improvement over $N_{ngb}=48$ is not
+obvious.\\ 
 
 To increase the convergence rate, one can use the derivative of the density with respect to the smoothing length in the
 Newton iterations:
@@ -145,7 +157,7 @@ The derivative of the kernel function has to be computed anyway to obtain a valu
 This term is given by
 
 \begin{equation}
-  \Omega_i = 1 + \frac{h_i}{3\rho_i}\sum_b m_b\frac{\partial W(\vec{r}_{ij},h_i)}{\partial h}
+  \Omega_i = 1 + \frac{h_i}{3\rho_i}\sum_j m_j\frac{\partial W(\vec{r}_{ij},h_i)}{\partial h}
 \end{equation}
 
 This concludes the first SPH loop in the standard implementation. More complicated quantities such as
@@ -266,8 +278,6 @@ function. It is just a convenient way to represent entropy.
 
 \section{Improved SPH equations}
 
-\textcolor{red}{WORK IN PROGRESS !!!\\}
-
 The equations \ref{eq:acceleration} and \ref{eq:dudt} correspond to a non-physical system with no viscosity and no
 thermal conduction. The physical model can be improved by adding some terms which vary depending on the authors. We
 follow here, D. Price and W. Dehnen.
@@ -300,12 +310,12 @@ in the equations. A decent place-holder for the second derivative of $W$ is then
 
  Artificial viscosity can be introduced by adding a term to both equation \ref{eq:acceleration} and \ref{eq:dudt}:
 
-\begin{eqnarray*}
+\begin{eqnarray}
   \vec{a_i} &\stackrel{visc}{=}& 2\sum_j m_j \frac{\alpha v_{sig}\left(\vec{v}_i -
-\vec{v}_j\right)\cdot\hat{r}_{ij}}{\left(\rho_i + \rho_j\right)}\hat{r}_{ij}\cdot \bar{F}_{ij} \\
+\vec{v}_j\right)\cdot\hat{r}_{ij}}{\left(\rho_i + \rho_j\right)}\hat{r}_{ij}\cdot \bar{F}_{ij} \label{eq:visc}\\
  \frac{du_i}{dt} &\stackrel{visc}{=}& -\sum_j \frac{m_j}{(\rho_i + \rho_j)} \alpha
 v_{sig}\left[\left(\vec{v}_i-\vec{v}_j\right)\cdot\hat{r}_{ij}\right]^2 \bar{F}_{ij}
-\end{eqnarray*}
+\end{eqnarray}
 
 where $\alpha$ is the dimensionless artificial viscosity and
 
@@ -318,21 +328,53 @@ where $\alpha$ is the dimensionless artificial viscosity and
 \end{equation}
 
 corresponds to the maximal (average) signal speed between pairs of particles.
-GADGET uses $\alpha=1$ (but can be changed) and $\beta=3\alpha$ (fixed). In addition, the Balsara switch is used.
+GADGET uses $\alpha=1$ (but can be changed) and $\beta=3\alpha$ (fixed). In addition, the Balsara switch is used. \\
 
 Modern implementations of SPH use a variable viscosity $\alpha_i$ for each particle. The idea behind this is to switch
 of viscosity in the part of the flows where the fluid is dissipation-less and to switch it on in shocks. This is done by
-using a shock detector and then a slow decay of the viscosity with time. Following Dehen, we use
+using a shock detector and then a slow decay of the viscosity with time. Following Price, we use
 
 \begin{eqnarray}
- \alpha_{loc,i} &=& \alpha_{max} \frac{h_i^2 A_i}{v_{sig,i}^2 + h_i^2 A_i} \\
- \dot\alpha_i   &=& \frac{2lv_{sig,i}\left(\alpha_{loc,i}-\alpha_i\right)}{h_i} \\
-  A_i &=& \xi_i \max\left(-\frac{d}{dt}\nabla\cdot v_i\right) \\
-  v_{sig,i} &=& \max\left(\frac{1}{2}(c_i+c_j) - \min(0, \left(\vec{v}_i-\vec{v}_j\right)\cdot \hat{r}_{ij})\right) \\
+ \frac{d\alpha_i}{dt} &=& \left(\alpha_{max}-\alpha_i\right)\mathcal{S}~+\frac{(\alpha_i-\alpha_{min})c_i\sigma}{h_i}\\
+ \mathcal{S} &=& \max\left(0, -\nabla\cdot \vec{v}_i \right)
 \end{eqnarray}
 
-The constants are usually $\alpha_{max} =2$, $l=0.05$ and $\xi_i$ = 1.
+with (usually) $\alpha_{max} = 2$, $\alpha_{min} = 0.1$ and $\sigma=0.1$. The $\alpha$ term in (\ref{eq:visc}) is then
+replaced by $\bar\alpha = \frac{1}{2}(\alpha_i + \alpha_j)$. The same applies to $\beta$.\\
+The divergence of the velocity field can be computed in the density loop and the exact expression is
+
+\begin{equation}
+ \label{eq:div_v}
+ \nabla\cdot v_i = \frac{1}{\rho_i}\sum_j m_j \left(\vec{v}_j - \vec{v}_i\right)\cdot \nabla_r W(\vec{r}_{ij},h_i) 
+\end{equation}
 
 \subsection{Thermal conductivity}
 
+The thermal conductivity dissipating energy at contact discontinuities can be modelled by adding another term to the
+internal energy equation \ref{eq:dudt}. 
+
+\begin{equation}
+ \frac{du_i}{dt} \stackrel{cond}{=} - \sum_j \alpha_u v_{sig,u}\left(u_i - u_j\right)\bar{F}_{ij}
+\end{equation}
+
+This time, the signal velocity must vanish when we are dealing with a contact discontinuity. A good choice is to used
+
+\begin{equation}
+ v_{sig,u} = \sqrt{\frac{2|P_i-P_j|}{\rho_i+\rho_j}}
+\end{equation}
+
+Once again, the $\alpha$-term is made to decay far from any discontinuity. In this case, the equation reads
+
+\begin{equation}
+ \frac{d\alpha_{u,i}}{dt} = \frac{h_i \nabla^2 u_i}{10} - \frac{\alpha_{u,i}c_i\sigma}{h_i} 
+\end{equation}
+
+where all constants usually take the same value than in the viscosity terms. The laplacian of $u$ can again be computed
+in the density loop and reads
+
+\begin{equation}
+ \nabla^2 u_i = \frac{2}{\rho_i} \sum_j m_j \left(u_j - u_i\right) \frac{F_{ij}}{|\vec{r}_{ij}|}
+\end{equation}
+
+
 \end{document}