A Stable Hybrid Potential–SPH Technique to Enforce the Fluid Incompressibility

A Stable Hybrid Potential–SPH Technique to Enforce the Fluid Incompressibility JUAN J. PEREA University of Seville School of Computer Engineering Avda. Reina Mercedes, s/n, 41012 Sevilla SPAIN [email protected] JUAN M. CORDERO University of Seville School of Computer Engineering Avda. Reina Mercedes, s/n, 41012 Sevilla SPAIN [email protected] Abstract: The SPH method has extensively used in fluid flow simulation. Through SPH the fluid is modelled by a particles system whose mutual interaction is weighed by a function, named kernel function, whose limits define the neighbouring of each particle. In spite of the high capabilities of SPH for simulate complex environments, it shows shortcomings specially if the fluid is subject to high changes in the pressure, the velocity and the density as occur in phenomena such as shock–tube, blast–wave or in the boundary and discontinuities where the number of neighbour particles is relatively low. In this case, the pressure gradient is inaccurate. As consequence, the simulation is instable with an erratic behaviour of particles. To avoid this problem, we propose a hybrid technique. This one consists in formulating the pressure gradient from a potential defined on each particles pair. Thus, the pressure gradient is immune to the low number of neighbour particles. Also, our proposal allows enforcing the fluid incompressibility. To show the improvements obtained we will carry out a set of simulations. Key–Words: fluid simulation, SPH, stability, accuracy, potential–based. 1 Introduction The fluid flow simulation plays an essential role in several disciplines both the scientific and technical. From an analytic point of view, it is formulated by Fluid Mechanics that set the partial differential equations describing the fluid flow. However, its complexity only allows obtaining analytic results for a restricted set of cases with high spatial symmetry. In most cases, the numerical methods must be used. These ones are studied by Computational Fluid Dynamic (CFD). Through CFD, the partial differential equations, are transformed into a set of algebraic equations, whose solution allow simulating the fluid flow. One of these methods is Smoothed Particle Hydrodynamics (SPH) which is widely used for its adaptability and versatility [21, 32]. SPH operates on fluid discretization through particles that interact among themselves. Thus, it is based on the Lagrangian formulation. Descriptively, each particle depicts a discrete portion of the continuum fluid and moves with the flow. In this methodology, the dynamic quantities are fit on each particle and its mutual interaction is weighed by a scalar function dependent on distance. This function, known as kernel function, must satisfy a set of analytical properties such as it must be monotonically decreasing, it must be normalized and its domain, named supported domain, should be limited. Each particle has got an associated kernel function, whose supported domain limits the number of neighbour particles. The SPH method has highlight features such as the conservation of mass [13, 19, 21], its adaptability and versatility to simulate complex environments [19] and the simplification of the formulation to solve the fluid dynamic equations [9, 13]. Nevertheless, despite these advantages, the simulation by SPH can show stability problems that affect the realism of simulation [3, 6, 28]. There are several instabilities types: the local mixing instability (LMI) [25], the heating instability [25], but the most remarkable is related to pressure gradient [20, 24, 25]. The inaccuracy of pressure gradient induces a tensile instability [14, 18, 21] or pairing of particles [2, 3, 24]. There is a consensus that relates the instability of pressure gradient and the number of neighbouring particles, specially when this number is relatively low or the particles distribution is not smooth [3, 10, 21, 24, 28]. According the number of neighbouring particles decreases, the inaccuracy of pressure force is increased. Consequently, the particles are scattered and show erratic behaviour. Usually, it has been regarded that increasing the radius of the supported domain the inaccuracy is reduced. Nevertheless, the increase in the supported domain re- WSEAS TRANSACTIONS on FLUID MECHANICS Juan J. Perea, Juan M. Cordero E-ISSN: 2224-347X 50 Volume 13, 2018 duces the resolution and, as consequence, the accuracy [3, 21, 32]. Furthermore, there are regions, such as the boundaries of fluid or its discontinuities, where the number of neighbour particles cannot be increased [28]. Other methodologies have been developed to avoid the problems associated to the low number of neighbouring. These techniques either reformulate the dynamic equations modelled by SPH [21, 25] or positioning control points to reduce the effects of the external loads on the pressure gradient [4, 33]. Nevertheless, these proposals show shortcomings since can violate the conservation of momentum [3, 16], particularly in boundaries and discontinuity regions. Also, they are only effective under specific conditions, if these conditions are not satisfied, an erroneous overdamping is induced [3]. To avoid these shortcomings related to the pressure force and the incompressibility enforcement, we propose a hybrid technique potential–SPH. Succinctly, our proposal is based on a formulation of a conservative potential that allow us to derive the pressure force. The constraints to obtain this force has the aim the improving the stability and accuracy specially in the discontinuities and boundary regions. On the other hand, the formulated pressure force allow us to enforce incompressibility without the need to use a predictor–corrector process what improves the efficiency. To show the improvements that our proposal allows obtaining, we will carry out a set of test usually used to quantify the stability and accuracy. 2 SPH Method The SPH method, initially developed by Gingold and Monaghan [5] and Lucy [11], was formulated from the integral formulation of the delta Dirac function [13]. f(r) = Zf(r0)δ(r−r0)dr0,(1) where fdepicts any scalar function and δis the delta Dirac function. According to the interpolation theory, the function δcan be replaced by a function with spatial extension. Thus, the integral expression (1) can be reformulated by the below equation. f(r) = ZΩ f(r0)W(r−r0, h)dr0+O(h2),(2) where Wis named the kernel function, Ωrefers to the volume of the definition domain of Wand his the radius of the supported domain. Regarding the kernel function, this one must satisfy, at least, the following constraints: 1. It must be normalized and tend to δwhen htends to zero, i.e: lim h→0W(r−r0, h) = δ(r−r0); ZV WdV 0= 1. (3) 2. It must be positive and decrease continuously with (r−r0). 3. It must be symmetric with respect to (r−r0). Quantitatively, the modelling of a continuous medium using a particles system entails transforming the integral, equation (2), by a sum where the mass of each particle depicts a volume element whose mass is ρdV , i.e: hf(r)i=ZΩ f(r0) ρ(r0)W(r−r0, h)ρ(r0)dr0 ≈X j∈N(i) mj fj ρj W(rj−ri, h),(4) where hi refers to the approximated value of function fsince the second order term has not been considered, N(i)depicts the neighbour particles jof each particle iand ρjis the mass density allocated to the particle j. From the symmetry property (item 3), the gradient of any scalar magnitude f(r)can be simplified as: h∇f(r)i=∂ ∂r Zf(r0) ρ(r0)W(r−r0, h)ρ(r0)dr0 ≈X j∈N(i) mj fj ρj∇W(rj−ri, h).(5) Generalizing the equation (5), it is possible to formulate a valid equation for any differential order. This generalized equation is: h∇lf(ri)i=X j∈N(i) mj fj ρj∇lW(rj−ri, h),(6) where lrefers to the differential order equation. It is noted that the equation (6) does not regard the second order term [21]. Also, it violates the conservation of angular momentum [13]. As consequence, although the equation (6) satisfies the SPH premises, the accuracy of the obtained results can be compromised. To prevent this issue, this equation must be modified to be applied to fluid flow simulation, according is described by [13]. Nevertheless, this symmetrical formulation can increase the instability [2, 7, 21]. WSEAS TRANSACTIONS on FLUID MECHANICS Juan J. Perea, Juan M. Cordero E-ISSN: 2224-347X 51 Volume 13, 2018 3 Dynamic Equations The dynamic of an incompressible and inviscid fluid, whose flow evolves adiabatically with high Reynolds number is described by the below PDEs system: ∇·v= 0,(7) Dv Dt =−1 ρ(∇p) + F, (8) where D/Dt =∂/∂t +v·∇refers to the convective derivative, ρis the mass density, vis the velocity, pis the pressure and Fdepicts the body force vector per unit volume. To close this system of equations, it is necessary to regard a relationship between the pressure and the mass density. This relationship is set by the equation of state (EOS). Although there are several formulations of the EOS, in SPH is usually used Tait’s equation expressed as: p= (γ−1)ρ, (9) where γis the adiabatic index. The modelling of the equations (7) and (8) by SPH shows problems of stability and accuracy. To solve these problems, several proposals have been developed [15, 21, 23, 30]. The most outstanding are the research developed by Price [21], appointed as MPM, and the technique described by Inutsuka [7] known as GPSH. Regarding the variational approach, Price [21] models the equation (8) as: dvi dt =−X j∈N(i) mj(χi∇Wij(hi) + χj∇Wij(hj)) + X j∈N(i) αmi ρij vsvij ·ˆ Rij∇˜ Wij +F(10) where the subscript iand jrefers to particle iand j respectively, χ= (Ω−1p/ρ2),Ωis a correction term related to the smoothing length, ˆ Rij is the unit vector defined by the line joining the particles iand j,˜ Wij is the averaged kernel and vs=ci+cj−vij ·ˆ Rij. In Price [21] is shown a detailed explanation of the derivation of this equation. On the other hand, the technique GSPH is seeking to restore consistency for the discrete density estimate by the kernel function convolutions using a Riemann solver. According to the GSPH technique, the equation (8) is formulated as: dvi dt =−X j∈N(i) mjp∗ ij hV2 ij(hi)∇Wij(hi√2)+ V2 ij(hj)∇Wij(hj√2)i ≈ − X j∈N(i) mjp∗ ij 1 ρ2 i∇˜ Wij 1 ρ2 j∇˜ Wij!+F(11) where p∗is the intermediate pressure arising from the solution of a Riemman solver particularized to the particles pair i, j and Vij are the intermediate velocities related to the solution of Riemann solver particularized to two particles [24]. The papers [7, 8] show a detailed explanation of the derivation of the equation (11). 4 Proposed Model In this section, it will be described the features of our proposal to improve the stability and accuracy of SPH simulations. Firstly, the fundamental hypothesis of the model will be set. From this hypothesis will be formulated our proposal to obtain the pressure force. Next, it will be described the analytic features that should be satisfied by the pressure force to guarantee stability. These features are set from the most outstanding studies in stable simulations scope by SPH. From these constraints and the fundamental hypothesis, it will be formulated the pressure force. Its formulation will allow enforcing the incompressibility. Once the incompressibility is enforced, it will be explained the process to calculate the density using both the SPH formulation and the pressure force. This joined formulation allows fulfilling the incompressibility constraint. Analytically, when external forces acting on a fluid, the distances between adjacent regions change. Then, into the fluid appears a force that tends to restore the balance conditions [34]. This internal force is quantified by the pressure gradient. Extrapolating this criterion to the modelled fluid by a particle system, we can formulate our fundamental hypothesis, namely: the pressure force is due to change of the distance between particles. Taking into account this premise, we avoid the problems of stability related to modelling of the pressure gradient by SPH, especially when the number of neighboring particles is low [3, 22]. Furthermore, we can enforce the incompressibility constraint on the pressure force and the density. On the other hand, as the force only depends on distance, we can formulate it from a conservative potential [34]. Nevertheless, WSEAS TRANSACTIONS on FLUID MECHANICS Juan J. Perea, Juan M. Cordero E-ISSN: 2224-347X 52 Volume 13, 2018 any conservative potential is not valid to obtain stable fluid simulations. To set a suitable formulation must be taking into account the constraints that guarantee the stability. These constraints it will be described below. 4.1 Pressure Force Several research have been carried out to set the relationship between stability–accuracy and the pressure gradient. The most highlighting studies assert that the pressure force: 1. It should ensure stability for a wide range of the number of neighbour particles. This feature can be deduced from the studies developed by [3, 12, 21, 28]. These research assert that if the number of particles is too low, the simulations show an erratic scattering of particles, mainly in the boundary of fluid. On the other hand, if the number of particles is high, a clustering of particles can appear. 2. It should be concave in respect of the abscissa axis. This constraint is deduced from the research developed by [1, 3, 24], where is described that a kernel function whose gradient, that appear in pressure force, shows an inflection point and concavity, improve the stability. 3. It should show an asymmetric profile with regard to the minimum value of the force. Thus, two requirements are satisfied: firstly, the repulsive force is increased when the distance is decreased, as it can be deduced from [9] and lastly, the interaction can be weighed by the distance, as is suggested by [3, 21]. A harmonic force can only satisfy the first condition but it would not ensure the last one [21]. 4. It should allow us to enforce the incompressibility [12, 28], namely, its formulation should limit the minimum distance among particles. Thus, the formulated force should show an asymptotic behavior with distance. Its aim would be to avoid that the distance between particles will be incompatible with incompressibility constraint. From these features, we formulate the conservative potential, from which, we will obtain the pressure force. The proposed potential is: Vi(ri) = 21 24πhiX j∈N(i) 2 (Rij −0.9ε)2− 7π2 2ζ(Rij −(0.85ε)) + (3π2σ)!,(12) where Rij =rj−ri,εis named as asymptotic factor, that allows us to control the incompressibility rate, ζ is appointed as depth factor, by this magnitude, we can avoid the particles clustering and σis named as tail coefficient, by means of this coefficient, it is possible to control the potential at the boundary of the supported domain. Empirically, we have evaluated the value range, of these parameters, that offer the best results, these values are: 0.3h≤ε≤h,0.5≤ζ≤2.0 and 0.8≤σ≤2.25. The Figure 1 shows the profile of the formulated potential. Figure 1: Graph of the proposed potential. The parameter values are h= 3,ε= 0.6,ζ= 1.25 and σ= 1.0. Analytically the potential, equation (12), only depends on position. As consequence, the force that is derived can be written as ~ Fpi=−∇V(ri), i.e: ~ Fpi=−5.25 6πhiX j∈N(i) 3.5π2 ζ(Rij −(0.85ε))2− 4 (Rij −0.9ε)3!(h−Rij)~ Rij Rij ,(13) From the Figure 1, we can describe the features of our proposed potential. Firstly, it shows a concavity in respect of the abscissa, as it was established in item 2. Secondly, it shows an asymmetric profile with respect to minimum of force, according to the item 3. Lastly, as a consequence of asymptotic limit, for low values of r, we can control the incompressibility, as it was required in item 4. It is noted, as consequence of the pressure force is derived from this potential, this features will be inherited. On the other hand, analyzing the equation (13) we can say that the force is saturated with the distance [34]. Furthermore, this saturation is weighed by the distance. Thus, the pressure force satisfies the constraint suggested by [3, 13, 17, 32]. Also, WSEAS TRANSACTIONS on FLUID MECHANICS Juan J. Perea, Juan M. Cordero E-ISSN: 2224-347X 53 Volume 13, 2018 the force formulated in equation (13) satisfies the conservation of momentum. Then, once the pressure force has been obtained, we are going to describe our proposal to control the incompressibility, taking advantage of the asymptotic limit shown by the formulated pressure force. Due to the fact that the pressure force is obtained from a potential, our technique will hereinafter be referred to as P–SPH. 4.2 The Enforcement of the Incompressibility Qualitatively, the incompressibility is related to the maximum force that can be borne by the fluid, without changing its local volume [34]. According to this fact, the incompressibility in a fluid modelled by a particle system can be explained as the minimum distance that the particles can approach. In this context, the force calculated by the equation (13), allows us to model the incompressibility by the asymptotic value when R is low. Nevertheless, if the incompressibility coefficient is lower than 2%, we must introduce a control process to reduce the numerical errors accumulated in each time step. Our proposal to control the incompressibility when its coefficient lower than 2% can be divided into two stages. In the first stage, the pressure force is calculated by the equation (13). In the last stage, we check the particles whose pressure force violates the incompressibility constraint. If the pressure force of any particle violates the incompressibility constraint, then the force excess is quantified and both the dynamic and the position of these particles are modified to satisfy the incompressibility constraint. To implement this process, firstly we calculate the maximum compression force, namely, the value from which the incompressibility constraint is violated. To this end, we particularize the equation (13) with a distance, whose value is in line with the rest density fluid. According to this description, the obtained equation is: ~ FC j|i=−5.25 6πh 3.5π2 ζ(ξ−(0.85ε))2−4 (ξ−0.9ε)3!~ R R, (14) where ~ FC j|iis the maximum of compression force that the particle jcan acquire in respect to particle iand ξ is the minimum distance that the particle jcan approximate of particle iwithout violates the incompressibility. The equation that satisfies this criteria is: ξ=3 s3m 4π(1 −η)ρ0 ,(15) where ηis the compressibility coefficient and ρ0is the rest density. If |~ FC j|i|>|~ Fpj|i|, then, we have to control the incompressibility. Thus, we change the force of the particle jand its position, to the maximum tolerated values, i.e.: ~ Fpj|i=~ FC pj|iy ~r0 j=ξ~ R R,(16) where ~ Fpj|iis the pressure force of particle jin respect to particle iand ~r0 jis its new position.This change only is valid if the difference between |FC pj|i|and |FC pj|i|is lower than 5%. Although this limitation seems restrictive, as the pressure force control the incompressibility, it allows us that the difference between two forces is lower than this value. 5 Results In this section, we will carry out a set of tests to show the improvements that offer our proposal. To that end, we will simulate a fluid whose dynamic is described by the equations (7) and (8). To solve this system of equations, it will be regarded the techniques MPM [21], GSPH [7] and P–SPH. The aim is to compare the obtained results from each technique to show the improvement that can be obtained by the proposed technique. It will be implemented three tests commonly used to quantify the stability and accuracy [24, 26, 27], these are: settling of a random particles distribution,Sod’s shock–tube and blast–wave. In each test, we will carry out an analysis of the range of the number of neighbour particles where each technique is stable. We have regarded the most of the values of simulation parameters suggested by [24]. Each test has different values that be specified in each simulation. Nevertheless, in all them, the value of the adiabatic parameter γ, equation (9), will be γ= 1,46. Regarding our proposed pressure force, equation (13), the selected values of εand ζare ε= 0.45 and ζ= 1.0. Likewise, all simulations will be implemented using a Wendland C4as kernel function. We have selected this function due to its analytic features that favour the stability and accuracy, according is highlighted by [3]. 5.1 Settling of a Random Particles Distribution Several research show a consensous which regards the settling of a random particle distribution as a suitable test to quantify the stability in SPH [3, 21, 29]. Descriptively, in this test, a particles system evolve from WSEAS TRANSACTIONS on FLUID MECHANICS Juan J. Perea, Juan M. Cordero E-ISSN: 2224-347X 54 Volume 13, 2018 P-SPH Figure 2: Evolution from a random particles distribution to a settled particles distribution. a random distribution to a state that shows a level of regularity. The achieved regularity allows quantifying the stability degree of any SPH simulation [21, 24]. To implement this test, we have used 980 particles whose mass is m= 0.00213. The simulation limits are [50 ×30] and the time step used ∆t= 1.5 10−4s and h= 0.023. From this values, we have obtained the simulation shown in Figure 2. To analyse the range of the number of neighbour particle that allow obtaining stable and accuracy simulation, we have changed from the lowest number to highest. In this test, the used values are: MPM GSPH P–SPH Ni120–375 95–375 15–375 Ni180 165 25 Table 1: Range of the number of neighbour particles. Nirefers to range that allows stability and Niis the value regarded to obtain the simulations shown in Figure 2. From the Figure 2, we can say that the best result is obtained from P–SPH. This one allows obtaining the results of the highest regularity. To quantify this lattice degree, we have formulated an order parameter ςthat satisfies: ς=Ps k=1(|dl−dk ij|) s,(17) where s= (m−1)(n−1), and mand nare the particles number in each space direction, dlis the lattice length, dk ij is the distance between each pair of neighbour particles. From the formulation of the parameter ς, it can be said that the lower value, the more lattice degree is achieved and, consequently, the more stability is obtained. The obtained values of this parameter are shown in Table 2. MPM GSPH P–SPH ς21.7 16.4 4.7 Table 2: Obtained values of the order parameter from Figure 2. 5.2 Sod’s Tube–Shock Test Sod’s shock–tube test [27] allow quantifying the accuracy and stability of SPH [10, 14, 26]. Descriptively, the shock–tube test is built by a tube where a non– porous membrane is introduced. This membrane defines two sections that are filled with gas. In a section, that is named as high section, the gas is around 10 times denser and more confined than in the other section which is named low section. From these initial constraints, the membrane is instantly taken away, at t= 0. At this moment, due to differences of pressure and density between high–low sections, a shock wave is generated. This one evolves into the low section. At the same time, a rarefaction wave appears into the high section. Between each wave, a contact discontinuity is arisen. Quantitetively, from the obtained accuracy in the simulation of these three regions, it is possible to quantify the suitability of the used technique. To carry out this test we have used a domain [-0.5,0.5]. The low section is limited between [-0.5,0.0)and the high section is limited by WSEAS TRANSACTIONS on FLUID MECHANICS Juan J. Perea, Juan M. Cordero E-ISSN: 2224-347X 55 Volume 13, 2018 (a) (b) (c) (d) (e) (f) Figure 3: Comparative graphs of the exact solution and numerical obtained values from MPM, GSPH and P-SPH for Sod’s shock–tube test. (0.0,0.5]. The remaining used paremeters are the suggested by [24]. From these initial values, we have obtained the results shown in Figure 3 at time t=0.15. It is regarded the density and the velocity to evaluate the obtain accuracy, according to the suggestion given by [24, 26]. The Figure 3 shows a comparison between the exact solution, solid line, and the numerical solution shows as dots. In each graphs, it is possible to identify three region: the contact discontinuity,the rarefaction wave and the shock wave. The contact discontinuity appears around xcd = 1.004, the rarefaction wave appear from x<xcd and the shock wave from x>xcd. Although the position of the contact discontinuity is the same in all graphs, both the spatial width of the contact discontinuity and its waving, shows significant discrepancies. To quantify the obtained accuracy, we have selected the numeric values of each part of the flow, i.e., rarefaction wave, contact discontinuity and shock wave, specifically the width that are appointed as: Wrt the width of the rarefaction tail, Wcd the width of contact discontinuity and Wsw the width of the shock wave. The obtained values are shown in Tables 3 and 4. 5.3 Blast–Wave Test The blast–wave test is a more extreme version of Sod’s shock–tube test and is formulated for a very MPM GSPH P–SPH r20.9651 0.9708 0.9891 Wrt 0.105 0.073 0.0205 Wcd 0.078 0.061 0.026 Wsw 0.128 0.104 0.051 Table 3: The most relevant values associated with the density graphs of Sod’s shock–tube test, Figs. 3a, 3b and 3c. MPM GSPH P–SPH r20.9651 0.9688 0.9891 Wrt 0.126 0.099 0.039 Wcd 0.101 0.075 0.043 Wsw 0.089 0.065 0.037 Table 4: The most relevant values associated with the velocity graphs of Sod’s shock–tube test, Figs. 3d, 3e and 3f. high Mach number, around M= 200. It was developed by Woodward and Colella [31] to highlight the behaviour of the contact discontinuity, mainly for the velocity, at high density and pressure values. In blast–wave test, the flow evolution shows the same features that Sod’s shock–tube test, namely, a contact discontinuity will appear inserted between the shock wave and the rarefaction wave. WSEAS TRANSACTIONS on FLUID MECHANICS Juan J. Perea, Juan M. Cordero E-ISSN: 2224-347X 56 Volume 13, 2018 (a) (b) (c) (d) (e) (f) Figure 4: Comparative graphs of the exact solution and numerical obtained values from MPM, GSPH and P-SPH for blast–wave test. As a result of the similarities with Sod’s shock– tube test, we have only changed the initial values of some simulation parameters. The new values of mass density, pressure and thermal energy are: (ρh, ph, uh) = (1.0,1 103,0.0) and (ρl, pl, ul) = (1.0,0.01,0.0), according to the suggestions given by [26, 31]. The new time step is ∆t= 8 10−5. Using these values, we have obtained the results that are shown in Figure 4 at time t= 7.5 10−4, where the exact solutions and numeric solutions are depicted by solid line and dots, respectively. Similarly to Sod’s shock–tube test, the three regions that arise when the membrane is taken away can be located in the graphs of Figure 4: the contact discontinuity position xcd = 0.1081, the rarefaction wave x<xcd and the shock wave x>xcd. Likewise, in each graph of 4 the width of contact discontinuity changes for technique obtaining the lowest width when our proposal is implemented. Analogously as the previous section 5.2 to shown the obtained accuracy, we have selected the most relevant values associated to the graphs of 4. These values are shown in Tables 5 and 6. 6 Conclusions In this paper, we have reviewed some of the most relevant features that influence on the accuracy and stability of flow simulation by SPH. From these fea- MPM GSPH P–SPH r20.9618 0.9791 0.9903 Wrt 0.094 0.081 0.046 Wcd 0.035 0.018 0.009 Wsw 0.094 0.079 0.038 Table 5: The most relevant values associated with the density graphs of blast–wave shock test, Figs. 4a, 4b and 4c. MPM GSPH P–SPH r20.9621 0.9788 0.9901 Wrt 0.0945 0.0813 0.0585 Wcd 0.0826 0.0582 0.0348 Wsw 0.087 0.026 0.009 Table 6: The most relevant values associated with the velocity graphs of blast–wave test, Figs. 4d, 4e and 4f. tures, we have proposed a hybrid potential–SPH based method. To show the improvements obtained by the proposed technique, we have implemented a set of tests using an incompressible and inviscid fluid. In each implemented test we have compared our proposal with the MPM and GSPH techniques which are formulated to avoid the instabilities related to the pressure gradient. Particularly, analysing the tests’ results, WSEAS TRANSACTIONS on FLUID MECHANICS Juan J. Perea, Juan M. Cordero E-ISSN: 2224-347X 57 Volume 13, 2018 we can conclude: 1. The P-SPH technique allows us to obtain a high stability. This conclusion can be deduced from the obtained results in settling test, Figure 2, where the obtained results by P-SPH show greater lattice degree than the obtained ones by the other two techniques. Quantitatively, this improvement is shown by the results of the parameter ς, Table 2, where the obtained values by P- SPH are the lowest ones. 2. The P-SPH allow us to obtain highly accurate results in simulations where the shock wave and contact discontinuity appear. This affirmation is based on the obtained results in Sod’s shock–tube test. The obtained graphs of the density and velocity, Figure 3, are more accurately simulated by P-SPH than by MPM or GSPH. This accuracy is more apparent in the rarefaction tail, the contact discontinuity and the shock wave. To reinforce this fact, we have selected the values associated to these regions that are shown in Tables 3 and 4. All obtained values evidence that the best accuracy is offered by P-SPH. 3. The improvements obtained in Sod’s shock–tube test is also accomplished in a more extreme version of this test. This conclusion can be obtained from the results shown in the blast–wave test. Qualitatively, the obtained results are in agreement with the exact solutions. Nevertheless, the most accurate results are obtained by the P-SPH, as can be seen in the graphs of Figure 4. Quantitatively, the improvement is shown in the values of Tables 5 and 6. 4. The P-SPH allow us to obtain stable simulations for a wide range of the number of neighbour particles. Also, the stability is achieved with the lowest number, as consequence. The P–SPH not only favour the stability but also the efficiency, from the computational cost point of view. This conclusion can be deduced from Table 1. From these results, we can conclude that the P- SPH improves both the stability and the accuracy of fluid flow simulations by SPH Acknowledgements: This research has been supported by the Pololas project (TIN2016-76953-C3-2- R) of the Spanish Ministry of Economy and Competitiveness. References: [1] R. M. Cabez´ on, D. Garc´ ıa-Senz, and A. Relano. A one–parameter family of interpolating kernels for smoothed particle hydrodynamics studies. Journal of Computational Physics, 227(19):8523–8540, Oct 2000. [2] S. J. Cummins and R. M. An improved sph method for modeling liquid sloshing dynamics. Journal of Computational Physics, 152:584– 607, Jul 1999. [3] W. Dehnen and H. Aly. Improving convergence in smoothed particle hydrodynamics simulations without pairing instability. Monthly Notices of the Royal Astronomical Society, 425(2):1068– 1082, Apr 2012. [4] C. T. Dyka, P. W. Randles, and R. P. Ingel. Stress points for tension instability in sph. International Journal for Numerical Methods in Engineering, 40(13):2325–2341, Jul 1997. [5] R. Gingold and J. Monaghan. Smoothed particle hidrodynamics: Theory and application to non–spherical stars. Royal Astronomical Society, 181:375–389, 1977. [6] M. Ihmsen, J. Orthmann, B. Solenthaler, A. Kolb, and M. Teschner. Sph fluids in computer graphics. In Eurographics 2014 - State of the Art Reports. The Eurographics Association, 2014. [7] S. Inutsuka. Reformulation of smoothed particle hydrodynamics with riemann solve. Journal of Computational Physics, 179(1):238–267, Jun 2012. [8] K. Iwasaki and S. Inutsuka. Smoothed particle magnetohydrodynamics with a riemann solver and the method of characteristics. Monthly Notices of the Royal Astronomical Society, 418(3):1668–1688, Dec 2011. [9] G. Liu and M. Liu. Smoothed Particle Hydrodynamics:A Meshfree Particle Method. World Scientific, 2003. [10] M. Liu, G. Liu, and K. Lam. Constructing smoothing functions in smoothed particle hydrodynamics with applications. Journal of Computational and Applied Mathematics, 155(2):263– 284, Jun 2003. [11] L. B. Lucy. A numerical approach to the testing of the fission hypothesis. Astronomical Journal, 82(12):1013–1024, Dec 1977. WSEAS TRANSACTIONS on FLUID MECHANICS Juan J. Perea, Juan M. Cordero E-ISSN: 2224-347X 58 Volume 13, 2018