Validation of robust SPH schemes for fully compressible multiphase flows

The present study examines the ability of the SPH number-density scheme to treat multiphase problems of the fully compressible regime. The number-density scheme is extended to the fully compressible regime, using the standard variational SPH framework and incorporate artificial diffusion coming from a generic formula. Aiming at robust schemes, we adopt the differential form of mass conservation. The performance of this scheme is studied with the help of two benchmark tests. It is shown that the standard variational framework of SPH may treat multiphase processes in the fully compressible regime, without reverting to non-standard formulations. The SPH solutions are compared to solutions coming from the Arbitrary Lagrangian Eulerian method and are validated against exact solutions.


INTRODUCTION
Multiphase problems in SPH have received substantial attention, due to SPH's straightforward way of introducing more than one fluids in the computational domain. Especially for the weakly compressible regime, multiphase algorithms have been extensively studied [1,2,3,4,5,6] and remedies have been pinpointed and often fixed. Validated results have been reported involving density ratios of up to 1,000 and sound-speed ratios of almost 10, by simply using standard SPH algorithms [6].
In all these schemes, it is common practice to use particles of different masses for each phase, such that particle mass ratios correspond to the initial density ratios. However, for the simulation of processes in the fully compressible regime, particles of equal masses are advised, so that the computed density depends on the local number density [7,8].
In order to assign equal masses to all particles, the ratio of the initial volumes per phase should be the reciprocal ratio of the phases' densities. Thus, the resolution in the lowestdensity region regulates the discretization length and hence the number of particles in the system. Therefore, simulations involving high density ratios-like gas-liquid shock tubesare either computationally implausible or are bound to be coarse.
Aiming at furnishing robust schemes, in the present work we explore the possibility of using particles of unequal masses within the standard SPH framework of [5,8] for the fully compressible regime. Additionally, we choose to examine the differential form of mass conservation of the number-density scheme, as a computationally efficient alternative to the integral form. In the integral form of mass conservation, the coupling of density with the smoothing length requires the iterative solution of a non-linear algebraic equation at every timestep [8]. This procedure does not necessarily converge when high density ratios are treated.
Finally, we examine the developed scheme employing tests which refer to shock propagation. The first test is the classical one-dimensional shock-tube, while the second test refers to a liquid-gas shock tube, with pressures described by the stiffened-gas equation of state. Results are shown also for a two-dimensional setup, in what can be referred to as a shockchamber test. Initially, air in a square geometry at high pressure is encapsulated within air of lower pressure. The SPH solutions are validated against the corresponding analytical solutions for the one-dimensional tests and the solutions obtained with an Arbitrary Lagrangian-Eulerian (ALE) strategy which utilizes fluid mixing theory, as the reference solution.

EQUATIONS
Casting the number density estimate into the SPH variational framework of [5,8] for the fully compressible regime: (1) we obtain its evolution equation: (2) and the momentum equation: Finally, the evolution of each particle's internal energy is provided by Gibbs fundamental thermodynamic relation: (4) with pressures P i = P(ρ i , e i ) given by an equation of state.

ARTIFICIAL DISSIPATION
Artificial viscosity compensates for errors caused by the subtle assumption of a differentiable Lagrangian function [9,8]. When using the differential form of mass conservation, artificial dissipation should also be used in order to capture shocks in density and smear out instabilities in the distrbution of pressure across contact discontinuities [10]. A generic form of artificial diffusion for each conserved variable appears in [11,5,9]. The following diffusive mass flux term is added in the equation of mass conservation: while the standard artificial viscosity is added to the momentum equation: Both are switched on only for approaching particles (v ij ⋅ x ij < 0). Signal velocities are and the optimal parameters are chosen as α ρ = 0.5, β ρ = 1.0, α v = 1.0 and β v = 2.0. Note that we eschew the mass-flux term in the calculation of the evolution of thermal energy Eqn (4). By that alone, the necessity for an artificial term in the evolution of energy is relaxed, at least for the Air-Air shock-tube. We observed that including the dissipative mass-flux term in the evolution of energy created a lag in the propagation of the shock. Finally, we also examine the artificial conductivity term [8]: which in the present study is added only to approachng particles, with the v sig, e from the formula above for parameters α e and β e .

TIME INTEGRATION AND INITIAL PARTICLE CONFIGURATION
Integration in time is achieved with a leap-frog scheme, using a sufficiently small constant time-step Δt in one-dimensional computations and a varying time-step in two-dimensional computations. The properly scaled Gaussian [8] is used as the kernel function.

ARBITRARY LAGRANGIAN EULERIAN METHOD
A brief description of the ALE (Arbitrary-Lagrangian-Eulerian) formulation used in this paper is presented, additional details can be provided in [12]. To solve multiphase flow problems, an ALE multi-material formulation can be used where both species can be mixed in the same element. Such an element is referred to as a mixed element. In our application, an element may contain two different materials; fluids of high and low pressure, as shown in Figure 1.
A mixture theory must be used to partition the material inside each element and compute the volume weighted stress from the constitutive model of each material as described by [13]. In the ALE description, an arbitrary referential coordinate is introduced in addition to the Lagrangian and Eulerian coordinates. The material derivative with respect to the reference coordinate can be described in Eqn (8). Thus configuration time derivative leads to the ALE equations: where X i is the Lagrangian coordinate, x i the Eulerian coordinate, w i is the relative velocity. Let us denote the material velocity by v and the mesh velocity by u, thus the governing equations for the ALE formulation are given by the following conservation equations: • Mass Equation: (9) • Momemtum Equation: (10) • Energy Equation: (11)

TESTS AND DISCUSSION
Two shock-tube tests in the one-dimensional domain x ∈ [0, 1] allow for an initial evaluation of the developed schemes performance in describing the propagation of shocks. The first test is the classical shock-tube benchmark, which has been addressed since the early development of SPH [14] and its solution within the standard variational SPH framework is extensively examined in [9]. The second test involves a multiphase medium, with discontinuous initial density distribution and fluid parameters as well. In both problems the stiffened-gas equation of state delivers the pressure: using appropriate ratio of heat coefficients γ and reference pressure P ref for each fluid. They are both Riemann problems and admit exact solutions, which are found via a procedure described in [15]. In addition to the analytical solution, a reference solution is obtained with Validation of robust SPH schemes for fully compressible multiphase flows the ALE multi-material formulation which is described above. In order to achieve equal "volumes" for the SPH system, the unit domain is divided in 400 intervals of equal length Δx o and one particle is placed in the middle of each interval. The mass of each particle is assigned according to the initial density distribution, as: m i = ρ o, i Δx o . In Figure 2 the solution at t = 0.2 is presented for the Air-Air shock tube. A constant timestep of Δt = 10 −4 is used. The analytical solution (black line) shows a constant pressure through the contact discontinuity. The suggested artificial dissipation terms manage to supress the singularities (inset plot for pressure in Figure 2) and create a continuous hump instead. Note that we use both artificial mass-flux and artificial conductivity with α e = 0.5 and β e = 1.0.
The second test is a Liquid-Gas shock-tube and involves an initial density ratio of 1/20 and a pressure ratio of 1/100. Additionally, note that the fluid parameters of the equation of state are discontinuous. In Figure 3, we see the results at t = 0.0023, obtained with a timestep of Δt = 10 −7 . The results in the plots are scaled with respect to the highest value of Int. Jnl. of Multiphysics Volume 9 · Number 3 · 2015 229  In the inset plot, the behavior of the SPH solution without the use of artificial massflux is shown each magnitude. Similarly to the Air-Air shock-tube, the solution for pressure suffers from the appearance of a hump. This effect is magnified by a factor of two in the inset plots of pressure. It is evident that the SPH solution manages to capture the left-running rarefaction and the shock which propagates through the low-density fluid. Nevertheless, there seems to be a tendency of the SPH solution to run ahead of the analytical solution. Apart from the one-dimensional experiments, the ability of the developed SPH scheme to describe patterns emerging from shock propagation is examined in a domain of two spatial dimensions, as well. To this end, both previous tests are performed in a twodimensional setup, in what can be referred to as a shock-chamber test. Initially, fluid of high density and high pressure is at rest in the square S in (x, y) = {|x| < 0.5 ∩ |y| < 0.5} and it is encapsulated within a fluid of lower density and pressure. The whole problem domain is the square S large (х, у) = {|x| < 3 ∩ |y| < 3}, where particles outside the square S(x, y) = {|x| < 2.5 ∩ |y| < 2.5} are boundary particles and at every timestep their velocity is kept fixed at zero. The use of the square geometry instead of a smooth circular one, makes the test more demanding.
In all figures, the SPH solution is plotted at particle positions with blue dots, the red dashed line refers to the ALE solution and the black solid line is the exact solution.
In order to realize the initial particle configuration, a number of particles per unit length  is chosen. This number defines a Cartesian grid for the square problem domain. Particles are placed in the middle of the grid's cells. The three spatial resolutions, which are examined, correspond to  ∈ {25, 50, 100} particles per unit length. The abrupt expansion of the highdensity fluid generates a shock which propagates through the low-density fluid and a rarefaction wave which moves towards the center of the chamber. After the expansion wave 230 Validation of robust SPH schemes for fully compressible multiphase flows The inset plot is a zoom by a factor of two of the solution around the contact discontinuity collapses to the centre of the chamber, it reverses its direction and starts moving towards the boundaries of the chamber. Thus, it provides a distinct wave pattern, which we investigate at t = 0.5. Lacking an analytical solution, we consider the solution obtained with the ALE multimaterial methodology as the reference solution. In order to focus our investigation on the aritificial mass-flux term, we neglect the term of artificial conductivity and we adopt the optimal parameters α ρ = 0.5 and β ρ = 1.0, similarly to the one-dimensional case.
In Figure 4, the upper triad of plots present the SPH solution for density, pressure and thermal energy respectively, with the finest resolution of 100 particles per unit length. The lower triad of plots show the ALE solution for the corresponding magnitudes. By and large, the SPH solution captures the wave pattern described by the ALE soultion. Nevetheless, the SPH solution tends to overestimate the zones of high pressure downstream the contact discontinuity.
For closer examination of the previous issue, we plot the problem's solution along the positive x-semiaxis of the domain, in Figures 5-7. The magnitudes of denisty, pressure and thermal energy appear from left to right respectively. The first triad of plots ( Figure 5), shows that the SPH solution does not diverge as the number of particles per unit length  = {25, 50, 100} increases. Besides this observation, we note that although the SPH solution follows the behavior of the ALE solution, there is an overestimation of all magnitudes to the left of the contact discontinuity (approximately at x = 1.0). In our tests, we found that if we include the mass-flux term in the evolution of internal energy Eqn (4) and for α ρ = 1.0, we are able to capture the correct magnitude at the expense of mispredicting the location of the shock.
The effect of the artificial mass-flux term is studied in Figure 5. The reason for this investigation is that the addtion of the artificial mass-flux term to the evolution of density, influences the calculated speed of sound c s = (∂P/∂p) s , via the computed thermal energy. A large value for α ρ might lead to misprediction of the sound speed. We examine the chosen value α ρ = 0.5 against the extreme values α ρ = {0, 1.0}. The results suggest that the predicted position of the shock is independent of the value of α ρ . Additionally, one may discern the catastrophic kinks in the plots of density and pressure, which appear if no dissipation for desnity is used.
Finally, the last triad of plots (Figure 7) shows the effect of adding a small amount of artificial conductivy (α e = 0.1 and β e = 0.2), along with the artificial mass-flux term for the base-case α ρ = 0.5 and β ρ = 1.0. Noteably, the extra artificial conductivity treats the small kink appearing in the plot of internal pressure, at the contact discontinuity.

CONCLUSION
The present study investigates the use of the number-density scheme extended to the fully compressible regime, via the variational framework of [5,8]. In order to endow robustness to the scheme, we adopt the differential form of mass conservation. Due to this choice, an artificial mass-flux term is necessary to counteract the spurious oscillations. The optimal values for the parameters of artificial dissipation terms are benchmarked with the help of the classical shock-tube test. Moreover, the developed shceme seems able to resolve wave structures which arise in a multiphase Liquid-Gas one-dimensional test. Nevertheless, there 232 Validation of robust SPH schemes for fully compressible multiphase flows Figure 6: The effect of α p , for the Air-Air shock chamber at t = 0.5 Figure 7: The coupled effect of artificial mass-flux and artificial conductivity, for the Air-Air shock chamber at t = 0.5 appears a tendency of the SPH solution to mispredict the analytically-derived location of the shock. Regarding the two-dimensional Air-Air shock-chamber problem, is reliable at capturing the corresponding wave patterns. Apart from an overestimation of the magnitudes downstream the cotact discontinuity, the overall behavior of the solution is good, compared to the reference solution, which is provided by the ALE methodology. Our results show that this behavior seems to be independent of the artificial mass-flux term, which manages to supress the emerging instabilities at the contact discontinuity. The next step in this research is to examine the corresponding two-dimensional Liquid-Gas test and evaluate the applicabilty of the developed scheme in more complicated problems of shock-propagation.