The role of time integration in energy conservation in Smoothed Particle Hydrodynamics fluid dynamics simulations

Herein we focus on weakly compressible flows, even though similar analyses can be carried out for incompressible flows, or even for different physical phenomena. Hence, the governing equation for the evolution of density is the conservation equation:

where $\left\langle{\cdot}\right\rangle$ denotes SPH operators, and abusing the notation, any magnitude resulting from the application of the SPH methodology. In the equation above $\rho_{i}$ is the density of an arbitrary $i$-th particle, and $\boldsymbol{u}_{i}$, its velocity.

The evolution of the velocity field is a discrete version of the Navier-Stokes momentum equation,

where $p_{i}$ is the pressure, $\mu$ the viscosity coefficient, $\boldsymbol{g}$ the acceleration due to external forces, and $\gamma$ the Shepard renormalization factor (see, for instance [41]). The extra term with the coefficient $k$, which appears in this discrete version and is absent in the continuum, is explained below.

In WC-SPH, these equations are closed by an equation of state (EOS) relating pressure and density:

where $p_{0}$ is the background pressure, $\rho_{0}$ is the reference density, and $c_{0}$ the speed of sound in the fluid. The latter is customarily set to a value high enough that the fluid behaves almost as if incompressible.

Incidentally, it may be highlighted that every single SPH related operator can be split in 2 terms,

i.e., all SPH operators can be split in the effect of the interactions with other fluid particles and the interactions with the boundary. The latter may adopt a number of forms [40]. For the sake of simplicity, only the first, “bulk”, part of the operators is considered in the main body of this article, and boundary effects are deferred to A.

Regarding the viscous term $\left\langle{\Delta{\boldsymbol{u}}}\right\rangle$, it may take several forms [42, 43], and may even include bulk viscosity [44]. As with boundary effects, we leave this term undefined, for the sake of simplicity and generality.

The momentum equation (2) differs from the traditional SPH formalism (as e.g. in the works of Monaghan [1] and Violeau [45]), by an extra term involving $\left\langle{\nabla{\gamma}}\right\rangle^{\Omega}$, whose form is:

where $m_{j}$ the mass of an arbitrary $j$-th particle, $W_{ij}:=W(|\boldsymbol{r}_{i}-\boldsymbol{r}_{j}|)$ is the value of kernel function for particles $i$ and $j$, and $\nabla W_{ij}:=\nabla W(|\boldsymbol{r}_{i}-\boldsymbol{r}_{j}|)$. This term consistently vanishes in the continuum, but not in the discrete formulation of the equations, except for highly ordered particle distributions, such as those often set up at the beginning of a simulation. This form of momentum equation is introduced by Colagrossi et al. [26], where the consequences of such a term are thoroughly discussed.

Equations (1) and (2) feature discrete divergence and gradient operators. For these, anti-symmetric SPH operators $\left\langle{\nabla\cdot{\boldsymbol{u}}}\right\rangle$ and $\left\langle{\nabla{p}}\right\rangle$ may be considered, of the form:

which are consistent with an $O(h^{2})$ error.

If the value $k$ is set to $2$ in Eq. (2), the resulting equation is effectively equivalent to using the popular symmetric pressure gradient operator [1, 45]

This formulation has been already considered by Colagrossi et al. [26] in the context of setting an optimal initial particle distribution. A symmetric pressure gradient operator grants immediate linear and angular momentum conservation, as well as power conservation [46]. However, we prefer to keep this term and the anti-symmetric gradient operator as two distinct terms, since the interplay between the two will be crucial, as further discussed below. The symmetry properties of these operators also play a main role in multiphase flow simulations and free-surface modelling, as discussed by Colagrossi et al. [47].

Finally, even if a linear EOS is suggested in Eq. (3), other expressions have been traditionally accepted in SPH, provided this more general relationship holds:

Analyses of power balance in SPH have been carried out several times in the past, e.g.: [37, 40, 39]. Every time that the balance is revisited, new terms are included, and a better arrangement of the existing ones is proposed. In this work, power balance is critical, so its development is carried out from the basics in this section. Incidentally, it should be highlighted that in this work “power balance” refers to the instantaneous balance, which only depends on the spatial SPH operators, while “energy balance” involves the time integration scheme.

Thus, temporarily accepting the acceleration term computed in Eq. (2) as the actual rate of variation of the particle velocity, the following kinetic power term can be defined,

which is indeed the discrete version of the continuous kinetic power. This kinetic power can be conveniently split in several terms applying Eqs. (2) and (4),

The first three terms above, Eqs. (12-14), are just discrete versions of well known terms in fluid dynamics. To be precise: $\left\langle{P_{\nabla{p}}}\right\rangle$ corresponds to the work due to pressure, $\left\langle{P_{\mu}}\right\rangle$, to the power due to the viscous forces, and $P_{p}$ to power from conservative external forces.

The extra term (15) is due to the $\left\langle{\nabla{\gamma}}\right\rangle$ operator, and obviously vanishes if $k=0$. Paradoxically, the vanishing of this term is not beneficial for the numerical scheme [46]. This can be easily checked out by rearranging the work due to pressure, $\left\langle{P_{\nabla{p}}}\right\rangle$, as the power due to compressibility, as it is usually done in fluid dynamics. To this end, the following identity, which is often used in order to prove that the anti-symmetric divergence (6) and the symmetric gradient (8) operators are skew-adjoint [48, 37] can be invoked:

See B for a straightforward calculation that shows that any pair of gradient and divergence operators that satisfy this identity will lead to the conservation of mechanical energy.

Eq. (16) can be written in terms of the anti-symmetric gradient (5) as

Using the identity of Eq. (17) in the work due to pressure, Eq. (12), and applying mass conservation, Eq. (1), the work due to pressure can be subsequently rearranged as

where $\left\langle{P_{\mathrm{c}}}\right\rangle$ is the power due to the compressibility.

Finally, all these terms may be introduced in the power balance, Eq. (18), which can then be rewritten

The term at the right hand side is an extra term, whose presence causes a lack of energy conservation (in addition to the term due to viscosity, which has a physical basis). Hence, it is possible to automatically ensure exact power balance if a symmetric pressure gradient operator ($k=2$) is chosen. In other words: the symmetric pressure gradient operator enables the model to non-physically exchange kinetic and compressibility energies, for the sake of total power conservation. Consistently, the same calculation for the conservation of mechanical energy referred to above will arrive at an expression that is not constant in general, unless $k=2$. Such power conservation feature of the symmetric pressure gradient operator was already discussed by Monaghan [46].

A number of works (e.g. [49], [47]) have targeted the SPH operator $\left\langle{\nabla{\gamma}}\right\rangle$. Some authors, [50] [24], proposed to diminish the value of such term by different means. However, its energy balance term, $\left\langle{P_{\nabla{\gamma}}}\right\rangle$, has not been addressed in the past.

The effect of boundaries is avoided in the present work, as this introduces additional terms to the equations, as discussed in A. Also, the SPH model is a simple one: other formulations such as $\delta$-SPH would add additional terms, see C.

In order to analyze a general integration scheme, the actual variation rates are computed a posteriori. These variation rates are defined, within a time step, as

where $\Delta t=t^{n+1}-t^{n}$ is the time step.

According to this, we define

Let us remark that the variation rates defined by (21), (22) and (23) can only be computed a posteriori, that is, after the flow fields at time $n+1$ are obtained by the actual time integration algorithm. Let us also point out that these variation rates serve as an analysis tool to asses the effect of the time integration scheme, but are not involved in the time integration itself.

With this in mind, an “equivalent Euler integration scheme” can be defined a posteriori. This “equivalent scheme” gives, within the n-th time step and using the rates defined by (21), the same result as the considered time integration scheme using the actual variation rates, given by the SPH equations, and applying the corresponding algorithm.

Considering the total kinetic energy of the discrete system, at any time step,

we can write the rate of total kinetic energy variation during consecutive time steps as

Developing this term one gets

Substituting $\boldsymbol{u}_{i}^{n+1}$ by its value, according to the equivalent Euler scheme, gives

Let us now compute the deviation of this power away from that of Eq. (10), caused by the errors in the velocity field due to the time integration scheme:

for any time $t\in(t_{n},t_{n+1})$. Substituting (27) and (10) in (28) we get

On the other hand, the velocities computed with the fully discretized SPH scheme do not agree with the ones resulting from solving the semidiscrete scheme, as the time integrations involve an error.

Furthermore, the SPH approximation of the derivative, for $t\in(t_{n},t_{n+1})$ can be split into two terms:

where $\boldsymbol{\delta}^{n}(t)$ is the vector difference, at time $t\in(t_{n},t_{n+1})$, between the SPH and the Equivalent Euler a posteriori variation rates. According to this, (29) can be rewritten as

The two terms in the sum above represent, respectively:

• 1.

  The error due to approximating $\boldsymbol{u}_{i}(t_{n})$ by $\boldsymbol{u}_{i}^{n}$, due to considering each $\left(\dfrac{\Delta{\boldsymbol{u}}}{\Delta t}\right)_{i}^{{m}^{\ast}}$ as the variation rate of $\boldsymbol{u}$ during the $(t_{m},t_{m+1})$ time step.

• 2.

  The error committed in the variation rate, as expressed by (30).

Let us now return to equation (20). It should be noticed that this energy balance is only valid if all the magnitudes involved are considered at the same time instant.

In order to evaluate the error using equation (29), we need to evaluate $\displaystyle{\left\langle{\frac{\mathrm{d}\boldsymbol{u}}{\mathrm{d}t}}\right\rangle_{i}(t^{\ast})}$ at some $t^{\ast}\in(t_{n},\,t_{n+1})$ This is therefore, computationally, a suboptimal approach. On the other hand, this methodology can be applied to any time integration scheme.

In order to close the balance after applying Eq. (20), similar error terms should be defined for the power due to the compressibility. To that aim, let us examine the nature of equation (27). The variation rate of the kinetic Energy defined by that expression can be seen as an extended kinetic power, built from the discrete velocity field. In effect, that equation represents an exact evaluation of the kinetic power if we consider a linearly extended velocity field. Using the variation rate for the velocity defined by (23), we define:

with $\tau\in[0,\Delta t]$. In effect, equation (27) corresponds to the exact kinetic power of a set of particles, moving in the way described by the field (32), evaluated at $t^{\ast}=t_{n}+\Delta t/2$.

Therefore, a way of building an error term for the power due to compressibility, $P_{\mathrm{d}\rho/\mathrm{d}t}(t)$, is to assume linearly extended fluid fields evaluated at $t^{\ast}=t_{n}+\Delta t/2$. With that in mind, consider the following extended fields according to the variation rates defined by (22) (23):

Therefore, defining

and inserting our definition for $P_{\mathrm{c}}^{n^{\ast}}$,

finally yields

This allows to introduce a final balance:

where $P_{\Delta t}$ is a combined deviation term:

Again, all terms at right hand side of Eq. (38) have to be considered as extra terms that cause a deviation from constant energy. It is therefore patent that $P_{\Delta t}$ enters as an additional deviation term, whose expression depends on the particular time integration scheme.

Regarding the time discretization, several schemes have been considered in this work:

1. 1.

   Explicit Euler scheme: it is the simplest integration scheme, in which values of the derivative at some time step are used to evaluate values at the next step. The scheme is first order and not A-stable. A-stability is a category of ODE solvers for which the solution is always stable, see [56].

2. 2.

   Explicit Heun method: also known as improved [56] or modified [57] Euler’s method, it is an explicit midpoint method. Simple predictor-corrector in which the explicit Euler method is used to obtain a prediction, which is then used in the correction. It is second order and not A-stable.

3. 3.

   Implicit midpoint method: a predictor-corrector scheme in which the mid point is guessed, and then iteratively refined. It is second order and A-stable.

In the particular case of the implicit midpoint method, the time in which the fields are evaluated is $t^{\ast}=t_{n+1/2}$. In order to solve the implicit equation at every time step, a solution strategy must be defined, such as an iteration to find a fixed solution. To this end, initial guesses $\lambda^{0}$ and $\mu^{0}$ for the velocity and density variation rates, respectively, are (we use superscripts as iteration indices):

which start an iterative process. The values of the extended velocity and density fields at each iteration are noted by $\hat{\boldsymbol{u}}^{m}$ and $\hat{\rho}^{m}$, respectively. The SPH approximations to these fields at each iterations are noted $\left\langle{\frac{\mathrm{d}\hat{\boldsymbol{u}}}{\mathrm{d}t}}\right\rangle_{i}$ and $\left\langle{\frac{\mathrm{d}\hat{\rho}}{\mathrm{d}t}}\right\rangle_{i}$. With this in mind, one has

In the expressions above $f^{m}$ is the relaxation factor at iteration $m$, which has been set to

where $M$ is the maximum number of iterations (either $10$ or $30$, in this practical application), and $\sigma_{f}$ a value tuned to get a small initial and final relaxation factor, $f^{0}=f^{M}=10^{-2}$.

This practical application has already been considered in order to assess energy conservation in SPH by Marrone et al. [58].

The initial condition is schematically depicted in Fig. 1.

Two identical water jets of dimensions $2H\times L$ each, with the same constant velocity $U$, but opposite directions, collide along the $y=0$ plane at the initial time instant, $t=0$. Inviscid flow with a initially constant density field, $\rho=\rho_{0}$, is considered for both jets. Also, no background pressure, $p_{0}=0$, or volumetric forces, $\boldsymbol{g}=0$ are considered. Thus, the whole phenomenon is characterized by the Mach number, $\text{Ma}=U/c_{0}$.

When incompressible flow is considered, the impact causes a sudden loss of a fraction of the initial energy dependent on the parameters chosen [59]. Along this line, Marrone et al. [58] demonstrated that applying $\delta$-SPH the simulation converges to the expected energy loss, while however not being an instantaneous process anymore. In contrast, in this paper an energy conserving impact is sought, in such a way that kinetic and compressibility energies are exchanged without any net loss, due to the lack of dissipative terms.

The Taylor-Green Vortex [62] is a popular practical application in which vortices are simulated in a periodic domain. The initial condition is schematically depicted in Fig. 13.

A square domain of size $L\times L$ contains four counter-rotating vortices, described by the initial velocity field

where $U$ is the characteristic velocity. Periodic boundary conditions are considered.

The corresponding pressure field is

with $\rho_{0}$ areplacednd $p_{0}$ the reference density and background pressure. The time evolution is analytically solvable, with the same expressions as above, but with $U_{0}$ replaced by $U_{0}\exp(-2\nu(2\pi/L)^{2}t)$, where $\nu$ is the kinematic viscosity coefficient.

In this particular case a Reynolds number, $\text{Re}:=LU_{0}/\nu=1000$, is chosen. As discussed below, larger Reynolds numbers can safely be considered when the implicit time integrator is applied. However, increasing the Reynolds number for the Euler explicit time scheme would require a extremely large number of particles, as well as an extremely small time step.