Lagrangian methods and Smoothed Particle Hydrodynamics
The classical Lagrangian method is the so-called Smoothed Particle Hydrodynamics method (SPH; [78] ,[79] ), which solves the Lagrangian form of the Euler equations (see Fig. 2) and can achieve good spatial resolution in high-density regions.
The basic idea of SPH is to discretize the fluid by mass elements (i.e. particles), rather than by volume elements as in Eulerian methods (see Fig. 2). It is evident that the mean interparticle distance in collapsed objects is smaller than in underdense regions. The scheme will thus be adaptive in spatial resolution by keeping the mass resolution fixed. For a comprehensive review with rigorous derivations of the SPH formalism, see [47] .
Fig. 5 summarizes the main characteristics of SPH. The general definition of a kernel smoothing method is the first step to build continuous fluid quantities for an arbitrary variable $X$. Here, the kernel depends on the distance modulus only, and is in addition required to be monotonic and differentiable (see upper part of Fig. 5). In the discretized version, we replace the volume element of the integration, $d {\vec x} = d^3 x$, with the ratio of the mass and density $m_j/\rho_j$ of the particles. Although this equation holds for any position ${\vec x}$ in space, we are only interested here in the fluid representation at the original particle positions ${\vec x}_i$, which are the only locations where we will need the fluid representation later on. It is important to note that for kernels with compact support (i.e., $W({\vec x},h)=0$ for $|{\vec x}|>h$), the summation does not have to be done over all the particles but only over the particles within the sphere of radius $h$ (see Fig. 2), namely, the neighbors around the particle $i$ under consideration. Originally, the most frequently used kernel is the $B_2$-Spline, but modern schemes invoke an entire new family of kernels like the HOCT kernels [80] or the so-called Wendland kernels [81] which show better stability and higher accuracy (see upper right panel of Fig. 5 and also recent review by [82] and references therein).
Derivatives can be calculated using the analytically known derivatives of the kernel. Conservation laws are numerically better achieved when pairwise symmetric formulations are used: this can be obtained by exploiting some identities when using derivatives of $\rho X$ or $X/\rho$ (see left middle part of Fig. 5). To profit from the adaptive nature of the SPH method, the smoothing length $h_i$ is typically allowed to vary for each individual particle $i$ and is determined by encompassing either a fixed number of neighbours or a fixed mass within the kernel. However, two complications often stem from the abovementioned individual smoothing lengths: first, a pairwise symmetric formulation needs to define an averaged kernel for each particle pair, and historically various ways to build averages have been discussed. The other is that the derivatives of the kernel -- being the smoothing length spatial dependent -- come with correction terms ($\partial h_i/\partial \rho_i$), which originally have been always ignored because they can't be computed directly. Starting from an entropy formulation, [83] derived for the first time an SPH formulation including the proper correction terms for the varying smoothing length from a Lagrangian formalism, which are listed in the middle right part of Fig. 5. Note that this formalism then also specifies the way in which the averaged kernel for a particle pair has to be constructed.
This can be further generalized, as shown in [84] . There, such a generalized formulation was obtained starting from an x-weighted volume average $$ \bar{y} = y_i = \sum_j x_j W_{ij}(h_i), $$ the (x-weighted) volume element $\Delta\nu_i\equiv x_i/y_i$ and the generalized relation $$ P_i = (\gamma-1)u_i\frac{m_i}{\Delta\nu_i} = A_i\left(\frac{m_i}{\Delta\nu_i}\right)^\gamma $$ between the pressure $P_i$, the internal energy per unit mass $u_i$ and the entropic function $A_i$. This leads to the set of generalized SPH equations $$ m_i \frac{d {\vec v}_i}{d t} = - \sum_j x_i x_j \left(f_{ij}\frac{P_j}{y_j^2}{\vec \nabla}_i W_{ij}(h_i) + f_{ji}\frac{P_i}{y_i^2}{\vec \nabla}_i W_{ij}(h_j) \right)\!, $$ and $$ f_{ij} \equiv 1 - \frac{\tilde{x_i}}{x_j} \left(\frac{h_i}{3\tilde{y_i}} \frac{\partial y_i}{\partial h_i}\right) \left[1 + \frac{h_i}{3\tilde{y_i}}\frac{\partial \tilde{y_i}}{\partial h_i}\right]^{-1}. $$ For the choice of $x_i=\tilde{x_i}=m_i$, which implies $y_i=\tilde{y_i}=\bar{\rho_i}$ and $\Delta\nu_i=m_i/\rho_i$, and following the entropy $A_i$ (e.g., $P_i=A_i\bar{\rho_i}^\gamma$), this set of equations will result in the entropy-conserving formulation of SPH as presented in [83] . For the choice of $x_i=\tilde{x_i}=(\gamma-1)m_i u_i$, implying $y_i=\bar{P_i}$ and $\Delta\nu_i=(\gamma-1)m_i u_i/P_i$, this results in a pressure--energy formulation of SPH, as presented in [85] . As now pressure is a kernel weighted quantity, contact discontinuities are properly treated. A third possibility is to choose $x_i=m_iA_i^{1/\gamma}$ which leads to a pressure--entropy formulation. For more details, see [84] .
The SPH method, however, is not performing as well in low-density regions as in the higher density ones. It also suffers from degraded resolution in shocked regions due to the introduction of a sizeable artificial viscosity. In its classical implementation, discretization errors introduce spurious pressure forces on particles in regions with steep density gradients in particular near contact discontinuities. This results in a boundary gap of the size of an SPH smoothing kernel radius over which interactions are severely damped. The standard implementation typically does not involve an explicit mixing term, which can compensate this effect. Therefore, the classical implementation does not resolve and treat dynamical instabilities in the interaction of multi-phase fluids, such as Kelvin--Helmholtz or Rayleigh--Taylor instabilities. Both these shortcomings can be overcome by a modern implementation as described at the end of this section. In addition, in the cosmological context, the adaptive nature of the SPH method, its simple way to couple to gravity and the possibility to have individual time steps often compensate for such shortcomings, thus making SPH still one of the most commonly used methods in numerical hydrodynamical cosmology.
$$ \frac{d {\vec v}_i}{d t} = - \sum_j m_j \left(f_j\frac{P_j}{\rho_j^2}{\vec \nabla}_i W_{ij}(h_j)+f_i\frac{P_i}{\rho_i^2}{\vec \nabla}_i W_{ij}(h_i) + \Pi_{ij} {\vec \nabla}_i {\bar W}_{ij} \right), $$ 以及 $$ \frac{d A_i}{d t} = \frac{1}{2} \frac{\gamma-1}{\rho_i^{\gamma-1}} \sum_j m_j \Pi_{ij} \left({\vec v}_j - {\vec v}_i\right){\vec \nabla}_i {\bar W}_{ij}. $$
This term describing an artificial viscosity is usually needed to capture shocks and its construction is similar to other hydrodynamical schemes. Popular formulations are those proposed by Monaghan and Gingold [86] and Balsara [87] , which includes a bulk viscosity and a von Neumann--Richtmeyer viscosity term, supplemented by a limiter reducing angular momentum transport in the presence of shear flows at low particle numbers [88] . Modern schemes implement a form of the artificial viscosity as proposed by [89] , based on an analogy with Riemann solutions of compressible gas dynamics. To reduce this artificial viscosity, at least in those parts of the flows where there are no shocks, a posibility is to follow the idea proposed by Morris and Monaghan [90] : every particle carries its own artificial viscosity, which eventually decays outside the regions which undergo shocks. A detailed study of the implications on the ICM of such an implementation can be found in [91] . There are various further improvements on the implementation a higher order artificial dissipation term [92] ,[93] . Modern SPH formulations also make use of higher order calculation schemes for velocity gradients, see lower part of Fig. 5 as well as related discussion in [94] ,[95] . Even better suppression of the artificial viscosity can be reached by following [92] in combination with such higher order calculation schemes for velocity gradients, as shown in [96] .
Considerable effort has been made in the last decade to combine the advantages of Langrangian and Eulerian methods. On one hand, significant progress has been achieved in involving Godunov methods into the SPH formalism, see [97] ,[98] ,[99] ,[100] , which recently led to the so called meshless methods. On the other hand, large effort has been made to extend Eulerian methods to moving mesh methods. Both these new methods are briefely described in the next two sub-sections.
评论