A physically consistent particle method for incompressible fluid flow calculation

In general, mechanical energy monotonically decreases in a physically consistent system, constructed with conservative force and dissipative force. This feature is important in designing a particle method, which is a discrete system approximating continuum fluid with particles. When the discretized system can be fit into a framework of analytical mechanics, it will be a physically consistent system which prevents instability like particle scattering along with unphysical mechanical energy increase. This is the case also in incompressible particle methods. However, most incompressible particle methods do not satisfy the physical consistency, and they need empirical relaxations to suppress the system instability due to the unphysical energy behavior. In this study, a new incompressible particle method with the physical consistency, moving particle full-implicit (MPFI) method, is developed, where the discretized interaction forces are related to an analytical mechanical framework for the systems with dissipation. Moreover, a new pressure evaluation technique based on the virial theorem is proposed for the system. Using the MPFI method, static pressure, droplet extension, standing wave and dam break calculations were conducted. The capability to predict pressure and motion of incompressible free surface flow was presented, and energy dissipation property depending on the particle size and time step width was studied through the calculations.


Introduction
Particle methods are widely used to calculate the complex motion of free surface flows in various engineering fields.
Smoothed particle hydrodynamics (SPH) for weakly compressible free surface flow was proposed by Monaghan [1] as an extension from astrophysics, while moving particle semi-implicit (MPS) was developed by Koshizuka and Oka [2,3] to calculate strictly incompressible free surface flows in the nuclear engineering field.
In designing a numerical methodology for physical simulation, it is important to take fundamental physics into consideration. In general, continuum mechanics, e.g., fluid dynamics, satisfies the fundamental laws of physics such as the second law of thermodynamics, which claims a monotonic decrease in mechanical energy. However, it is not always satisfied in a discrete system approximating the continuum equations. When a discrete particle system does not satisfy the second law of thermodynamics, the mechanical energy may increase and cause instability like particle scattering. Therefore, it is important to satisfy the fundamental laws of physics in formulating interaction forces in particle methods.
The physical consistency is taken care of using analytical mechanical frameworks in various scales of calculation. For example, molecular dynamics [4] and astrodynamics [5], where the energy dissipation is negligible, are constructed following the classical analytical mechanics, and the mechanical energy is conserved in their systems. Besides, dissipative particle dynamics (DPD) [6,7], which is for the mesoscale simulation where the thermal fluctuation is taken into consideration, is formulated based on the general equation for the nonequilibrium reversible-irreversible coupling (GENERIC) framework [8,9], and the system satisfies the first and the second laws of thermodynamics. As an extension of DPD to the larger-scale problems, Espanol and Revenga [10] proposed a thermodynamically consistent version of SPH, smoothed dissipative particle dynamics (SDPD), and they also applied the GENERIC framework [8]. These systems in various scales stand on the frameworks of analytical mechanics and satisfy the fundamental laws of physics.
The fundamental laws of physics are to be satisfied also in a discrete particle system for the continuum calculation. In the continuum scale, thermal fluctuation is usually negligible, and only the conservative force and dissipative force are to be considered. In fact, the governing equation for fluid dynamics, Navier-Stokes equation, usually consists of conservative pressure term and dissipative viscosity term. For the system having dissipative force as well as a conservative force, a simple analytical mechanical framework, which is an extension of Lagrangian mechanics [11], is available.
However, it is not easy to apply the theoretical frameworks for the physical consistency together with incompressibility, which is an effective approximation in many industrial situations where the Mach number is small. In fact, the DPD [6,7] and SDPD [10], which are the physically consistent mesoscale particle methods, do not satisfy the incompressible condition, even though the Mach number tends to be small in the smaller scale. The SPH method [1] is also not strictly incompressible, and it is usually expressed as weakly compressible. Conversely, most of the particle methods for strictly incompressible flows do not have the physical consistency, and currently, there are only a limited number of methods having such features. They are incompressible SPH (ISPH) proposed by Ellero et al. [12] and Hamiltonian MPS (HMPS) proposed by Suzuki et al. [13], respectively. They constructed incompressible particle methods within the framework of the energy-conserving system with geometric constraints. They used symplectic time integrator [14], SHAKE and RATTLE, respectively. However, their algorithm needed to solve nonlinear equations in every time step, which is not favorable for numerical efficiency.
Because of this difficulty in building up a particle method with both the strict incompressibility and the physical consistency, most of the incompressible particle methods scarify the physical consistency. The MPS method [2], which is the first strictly incompressible particle method, chose to solve a discretized version of pressure Poisson equation (PPE) in the semi-implicit algorithm which was similar to the one adopted in the finite difference method [15]. Subsequently, several researchers proposed the various PPE-based incompressible particle methods [16][17][18][19]. However, the PPE-based methods often suffered specific instability due to the ignorance of the physical consistency, and a lot of stabilization techniques are proposed such as the ones with respect to PPE formulation [20][21][22][23][24], particle regularization [25][26][27] and surface detection [28][29][30][31]. However, these PPE-based incompressible particle methods [2,[16][17][18][19] might suffer instability like particle scattering in some condition even with such stabilization techniques [20][21][22][23][24][25][26][27][28][29][30][31] because it is difficult to fit their formulations into the framework of analytical mechanics [11], and the systems may allow unphysical increase in mechanical energy.
In this study, a new incompressible particle method with the physical consistency, moving particle full-implicit (MPFI) method, is developed. The interaction forces can be fit into the extended Lagrangian mechanics [11], which is an analytical mechanical framework for the system with dissipation. Therefore, the discretized system in this method is physically consistent although the mechanical energy conservation as in the previous studies [12,13] is not satisfied. Instead of directly solving the constant density constraint, a control equation which constrains velocity divergence was adopted to keep the density constant. The control equation could be solved in a linear matrix system, where velocity and pressure are calculated implicitly at the same time. In addition, a new pressure evaluation methodology for the physically consistent system is proposed based on the virial theorem [32]. Using the MPFI method, static pressure, droplet extension, standing wave and dam break calculations are conducted. In the static pressure calculation, the results are compared with those of the MPS method [2] and the SPH method [1]. Moreover, the energy dissipation property of the MPFI method depending on the particle size and time step width is studied. In the droplet extension and standing wave calculation, the fluid motion is compared with the theoretical solutions. Finally, the fluid motion and pressure in the dam break calculation are compared with the experiment [33].

Governing equation
The Navier-Stokes (NS) equation with a bulk viscosity term is described as and where ρ, u, Ψ, g, λ and κ are the density, velocity, pressure, gravity, bulk viscosity and bulk modulus, respectively. Here, the shear viscosity term is not included in the NS equation because only the calculation cases without shear viscosity are going to be shown in this study. When an infinite value is assigned to the coefficients, λ and κ, Eq. (2) will be changed to

Discretization
In this section, the partial differential operators in the governing equations [Eqs. (1) and (2)] are firstly discretized. In particle methods, the partial differential operators are replaced with particle interaction models expressed by weighted averages using effective radius and weight function. Since this manner is common in the previous particle methods, the SPH method [1] and the MPS method [2], a similar interaction model is adopted in this study.
In the SPH method [1], a bell-shaped kernel function is used to define the continuum field, and the partial differential operators, such as the gradient (∇) and divergence (∇•), are derived from the field. As a result, the partial differential operators are replaced by the interaction models including the differential of the kernel function. The formulation with the kernel differential can be used to make a physically consistent system. However, the interaction model in the SPH method has a disadvantage on the instability with respect to particle clustering. Since it yields only small repulsive force when particles are so close to each other, the particles easily agglomerate. This is why the SPH method had to adopt some smoothing treatment such as the artificial viscosity [1].
On the other hand, in the MPS method [2], the particle interaction models for partial differential operators are directly formulated with a weight function. The weight function in the MPS method [2] gives a larger value when two particles get closer. This helps to avoid the clustering. However, the interaction model of the MPS method does not fit to the framework of analytical mechanics [11]. It is because the interaction model in the MPS method does not include the differential of the weight function, which will emerge when potential energy is differentiated following the analytical mechanical framework [11].
Therefore, in this study, a similar but new particle interaction model is adopted, taking the good points of the models in the SPH method [1] and the MPS method [2]. Here, the following normalized weight function w ij is used: where r e is the effective radius and d ij is the distance between particles i and j. The normalization parameter N 0 is calculated at the particle i around which the neighbor particles are well arranged. Using the weight function, the particle interaction models for the gradient and divergence operators are formulated as where ϕ, A and r ij are an arbitrary scalar, an arbitrary vector and the relative position between particles i and j, respectively. In Eq. (9), S ∇ is a normalization parameter and w ij ′ is the differential of the weight function with respect to particle distance d ij . These formulations include the differential of the weight function as in the SPH method [1]. This choice has an advantage in constructing a particle method within the analytical mechanical framework [11]. However, the shape of the weight function is more like that adopted in the MPS method [2]. The gradient of the weight function will be larger, when the particle distance d ij is smaller. This feature plays an important role in avoiding particle clustering because larger repulsive force emerges between closer particle pairs. Since the formulations in Eq. (9) are the weighted averages very similar to the ones used in the SPH and MPS (8) methods [1,2], it is expected that it can also evaluate the partial differentials when the normalization parameter S ∇ is calibrated properly. In fact, the parameter S ∇ can be calibrated using the virial theorem [32], which will be shown later [Eq. (32) in Sect. 2.6]. In addition to the particle interaction models for discretizing the partial differential operators, a parameter called particle number density is used to evaluate the fluid density. It is a summation of the weight function: It is assumed that the value is proportional to the fluid density, and it is desired to be constant when the fluid is regarded as incompressible. In this study, the base value for the particle number density is 1.0 because the weight function is normalized in Eq. (8).
With the particle interaction models [Eq. (9)] and the particle number density [Eq. (10)], the governing equations [Eqs. (1) and (2)] are discretized to and respectively, where ̃ is the particle velocity at the previous time step, and the backward Euler method is adopted for the time integration. In this study, the effective radius, r e = 2.5 l 0 , was adopted for Eqs. (11) and (12), where l 0 is the initial particle spacing.
With dividing Eq. (12) by λ, and taking the limit of, λ, κ → ∞, is obtained, where the ratio, γ = κ/λ, is applied. This is a discretized version of Eq. (3), controls the particle number density to be constant, n i = 1, and is equivalent to the discretized version of the null-velocity divergence condition [Eq. (6)]: and the discretized version of the constant density condition [Eq. (7)]: when the initial condition n i | t=0 =1 is given. However, in the region close to the free surface, the initial condition n i = 1 is not satisfied, and even the proportionality between the particle number density n i and density ρ cannot be assumed. In such a region, the null-velocity divergence condition [Eq. (14)] is to be applied instead. Therefore, the control equation [Eq. (13)] is replaced by and this equation is going to be calculated in the full-implicit algorithm shown in the next section.
The free surface boundary in the MPFI method is naturally given by the vacant space from which no force works to the particles. Therefore, no explicit pressure boundary is needed. Besides, to express walls surrounding the fluid particles, a symmetric boundary condition was adopted in this study.

Full-implicit algorithm
The discretized motion equation [Eq. (11)] and the incompressible control equation [Eq. (16)] are solved in a fullimplicit algorithm shown in Fig. 1. First, the velocity u and the pressure Ψ of the particles are calculated at the same time solving Eqs. (11) and (16). In contrast to the previous studies [12,13] where the constant density condition was directly posed and the nonlinear calculation was required, the alternative use of Eq. (16) can avoid the nonlinear calculation because Eqs. (11) and (16) form a linear system with a symmetric coefficient matrix. In this study, the conjugate residual (CR) method [34] was applied for the matrix solver. After the full-implicit calculation, the position x of the particles is updated explicitly as Fig. 1 Full-implicit algorithm in moving particle full-implicit (MPFI) method Although the shear viscosity term was not included in the governing equation in this study, the full-implicit algorithm can be an advantage for the incompressible high-viscous fluid calculation when the viscosity term is incorporated into the matrix equation. In a high-viscous calculation, the diffusion number usually restricts the time step width to be small. In contrast, the implicit velocity calculation may allow the usage of larger time step width because it is basically free from the restriction.

Incompressible control
The incompressible conditions [Eqs. (6) and (7)] are equivalent under the initial condition ρ| t=0 =ρ 0 , and their spatially discretized versions [Eqs. (14) and (15)] are also equivalent to each other under the initial condition n i | t=0 =1. When n i is differentiated with respect to time, is obtained. This is interpreted as a discrete version of the continuity equation [Eq. (4)] because the right-hand side of Eq. (18) is the same as the divergence model in Eq. (9) except for the coefficient, S ∇ . By substituting Eq. (18) into Eq. (14), can be derived. Therefore, Eqs. (14) and (15) are equivalent to each other under the initial condition, n i | t=0 =1, if the time discretization error could be ignored. However, they are different in a practical calculation. When the discretized nullvelocity divergence constraint [Eq. (14)] is adopted, the error accumulates and the density may gradually change. On the other hand, when the constant density condition [Eq. (15)] is adopted, the density is kept strictly constant. The reason for this difference is the time discretization error. Since the time differential, ∂n/∂t, remains in Eq. (18), time discretization error will be introduced in a practical calculation where the finite time step width has to be used. The equation adopted to control the particle number density in this study [Eq. (16)] is a hybrid type, and it can stop the error accumulation to some extent without solving nonlinear equations. To understand how the particle number density is controlled, it is helpful to derive the control equation with respect to the particle number density n i . By substituting Eq. (18) into Eq. (16), is obtained. Even though Eq. (20) is not directly calculated in the program, it underlies this calculation. Since the control equation [Eq. (20)] is not satisfied perfectly because of the time discretization error included in Eq. (18), the particle number density may deviate from the base value. When the fluid is compressed (n i > 1), the deviation is reduced by the upper line on the right-hand side of Eq. (20), and the gradual density change will be suppressed. On the other hand, when the fluid is expanded (n i < 1), the deviation may accumulate and the fluid density may gradually change.
The limitation for the parameter γ can also be known from Eq. (20). Denoting the change of particle number density in a single time step as Δn, the upper line of Eq. (20) can be written as Since the particle number density is explicitly calculated from the particle position, the magnitude of Δn should be smaller than the deviation, n i − 1, to avoid numerical oscillation. Therefore, is to be satisfied at least.

Physical consistency
For the system with conservative force and dissipative force, the extended Lagrangian mechanics for the system with dissipation [11] is useful. The Lagrangian equation is written as where L and D are Lagrangian and Rayleigh dissipative function, respectively. In the framework of the Lagrangian mechanics, the Lagrangian is defined by where T and U are the kinetic energy and potential energy of the system. Moreover, the dissipative function has to be positive to ensure the monotonic decrease in mechanical energy. If the discretized equations in the MPFI method can be fit into the framework of the extended Lagrangian mechanics [Eq. (23)], the physical consistency, i.e., the monotonic decrease in mechanical energy, is assured. Since Eq. (16) is a constraint and it is not easy to relate them to the framework [Eq. (23)], the relation to the equations before taking the limit of λ, κ → ∞ [Eqs. (11) and (12)] is firstly considered 1 3 instead. When the Lagrangian L and the dissipative function D are given as and the Lagrangian equation [Eq. (23)] coincides with the governing equations [Eqs. (11) and (12)], except that the time differential on the left-hand side of Eq. (11) is discretized. Thus, the equations before taking the limit of, λ, κ → ∞ [Eqs. (11) and (12)], can be fit into the framework of the analytical mechanics [11], and they form a physically consistent system with dissipative force. Since the coefficients λ and κ were the arbitrary parameters in deriving the governing equations [Eqs. (11) and (12)] from the Lagrangian [Eq. (25)] and the dissipative function [Eq. (26)], Eq. (12) can be replaced with Eq. (16) by taking the limit of λ, κ → ∞ or by further specifying the parameters as γ = κ/λ = 0. Therefore, it is expected that the discrete system expressed by Eqs. (11) and (16) is also a physically consistent system including dissipation.

Pressure evaluation based on the virial theorem [32]
The pressure Ψ i calculated via Eqs. (11) and (16) suffers large fluctuation. However, in a physically consistent system, a smooth pressure field can be obtained in a physically justifiable way based on the virial theorem [32] which is used in the molecular dynamics. For the smooth pressure field, the virial pressure P is evaluated in a post-process using the interaction force and particle position, which are obtained in the main calculation. The virial theorem [32] gives a relation between pressure, kinetic energy and interaction forces in a statistical manner as where s dim , P, V, K, F i and r i are the spatial dimension, virial pressure, volume of the space, kinetic energy of particles due to fluctuation, force acting on particle i and position of particle i, respectively. The blanket 〈〉 indicates the average with respect to time. The first term on the right-hand side of Eq. (27) is the fluctuation energy of the particles, which can be neglected in this study because neighbor particles mostly move along with each other. The second term on the right-hand side of Eq. (27) relates to the interaction forces among particles. Because the particle interactions in the present method are all defined by pair-wise forces conserving linear momentum, the second term can be rewritten as where F ij and r ij are the interaction force and relative position between particle i and particle j, respectively. The virial theorem evaluates the pressure in the whole system at equilibrium state by averaging the physical quantity with respect to time and space as in Eq. (27). However, pressure at a certain time and at a certain point is to be evaluated in continuum-scale fluid calculations. To evaluate such local values, local averaging is to be taken with assuming the local equilibrium. To take the local average, a test region and a test period are to be given to define a range for averaging. When the range is larger, the obtained value suffers from less fluctuation but has less local information. Therefore, the range for averaging is to be given properly depending on the desired resolution with respect to time and space.
To evaluate the local pressures, the summation in Eq. (28) must be taken over the test region. To obtain equations for the smaller region, the separation of the region is to be considered. A splitting with respect to Eq. (28) is shown in Fig. 2, where one region is divided into two smaller regions. The interaction forces traversing the boundary are halved and added to the summation in each region such that the summation in the original region should equal the sum of those in the two separated regions. This reversibility must be satisfied in dividing the regions. In this manner, the division can be repeated until the single-particle regions are obtained. The equation for the single-particle region with respect to particle i at a certain time step is expressed as where ΔV is the volume of the single-particle region, which is assumed to be constant as where l 0 is the particle distance at the initial state. By replacing the interaction force F ij in Eq. (29) with the pressure gradient force in the governing equation [Eq. (11)], is obtained. With Eq. (31), the virial pressure P i in the single-particle region can be calculated.
Additionally, Eq. (31) can also be used to calibrate the normalization parameter S ∇ . Since the virial pressure P i should be consistent with the pressures Ψ i and Ψ j , the values for P i , Ψ i and Ψ j should be the same when the pressure field is uniform. Therefore, Eq. (31) has to hold when the same value is inserted into P i , Ψ i and Ψ j . This condition determines the normalization parameter S ∇ as In this study, the initial lattice particle distribution was used for this calibration.
The virial pressure P i obtained via Eq. (31) is not smooth because the single-particle region is not large enough to apply the virial theorem. One reasonable choice for the larger test region is the region inside the effective radius. In this study, the region in the radius was adopted. The physical quantity in the larger region can be calculated by summing up the ones in the single-particle regions as where N {i|dij<re} is the number of particles inside the radius. Subsequently, the virial pressure in the region around particle i is calculated as This implies that the virial pressure for the test region can be evaluated by averaging the viral pressure for the singleparticle region. This physically justified averaging process [Eqs. (31) and (34)] enables to obtain smooth pressure field, even when the pressure, Ψ i , directly calculated via the discretized governing equations [Eqs. (11) and (16)] suffers from large fluctuations.
Since the virial theorem [Eq. (27)] is to be applied not only with the spatial averaging but also with the time averaging, it is also possible to take time averaging in a certain period as where p is the period for the averaging. In this study, no time averaging is adopted unless otherwise noted.
For the justification of the pressure evaluated by numerical calculations, the accuracy of the discretization scheme is often discussed by showing that the truncation term goes zero by reducing the particle size and time step size. This feature is usually called "consistency," but here, we call it "scheme consistency" to distinguish from the "physical consistency" shown in Sect. 2.5. The interaction models [Eq. (9)] do not always reproduce the gradient correctly even in a linear field mainly due to the non-uniform distribution of neighbor particles. Therefore, the models may have zeroth-order accuracy at most, and it is difficult to show the "scheme consistency" which stands purely on a mathematical discussion. However, the pressure in the "physically consistent" system can be justifiably evaluated through the virial theorem, without showing the "scheme consistency." This is one of the merits to achieve "physically consistency."

Hydrostatic pressure calculation
The basic property of the MPFI method was studied using a hydrostatic pressure calculation. The initial particle arrangement is shown in Fig. 3. The parameters used in the base case are shown in Table 1. For the comparison, the same geometry was calculated using the MPS method [2] and the SPH method [1], where MPS-SW-MAIN-Ver.2.0 (2010) and DualSPHysics_v4.4 (2019) were used, respectively. In the MPS calculation, the empirical parameters, such as the relaxation factor for PPE calculation and the threshold for surface detection, were tuned so as to avoid particle scattering. The effective radius and the time step width were set the same as the MPFI calculation. In the SPH calculation, the sound speed is set 10.0 m/s which is thought to be large enough for the static pressure calculation where the flow speed is basically zero. The artificial viscosity is set large enough to suppress the oscillation caused by initial unbalance. The cutoff radius, which is the double length of the smoothing length in the SPH method, was set the same as the effective radius in the MPFI calculation, and the smaller time step width, Δt = 0.0001 s, was applied so as to obtain a stable result in the explicit SPH calculation. The pressure distributions at t = 10.0 s obtained by the MPFI, MPS and SPH methods are shown in Figs. 4, 5 and 6, respectively. For the MPFI results, the pressure evaluated using Eq. (34) was shown, whereas for the MPS and SPH results, the pressure outputs from the program were directly shown. The smooth pressure distribution is obtained with all the three methods. The particles stayed at the initial position in the MPFI calculation because the static balance was made by the initial lattice distribution, whereas the particle movement from the initial lattice distribution was observed in the MPS and SPH calculations. The pressure history at the center (A in Fig. 3) is plotted in Fig. 7. The pressure obtained with the MPFI method was almost constant at the theoretical value. On the other hand, the pressure fluctuations were observed in the MPS and SPH calculations. In the MPS calculation, the fluctuation started along with the particle movement and was maintained at a certain level. The pressure on average was a little bit smaller than the theoretical value. Since one of the reasons for the underestimation is the usage of the zeroth-order scheme [2], the higher-order schemes [35] may help to improve the prediction. In the SPH calculation, the pressure periodically oscillated due to the initial unbalance and the oscillation gradually decayed. The decay was due to the artificial viscosity which was set large to suppress the oscillation. The pressure in the SPH calculation was overestimated. The reason for the overestimation was mainly because of the density calculation, which is discussed later.   The density deviations in the MPFI, MPS and SPH calculations are shown in Figs. 8, 9 and 10, where the red color shows the particles with larger density than the base value and the blue color shows the particles with smaller density. In Figs. 8 and 9, the deviation of the particle number density calculated by the MPFI and MPS methods is shown, respectively. In the MPFI calculation, the particle number density was kept constant because the particles moved very little. In the MPS calculation, the deviation was around 2% although the whole system was kept almost incompressible. In contrast, Fig. 10 shows two types of density. One is continuum density, which is calculated by integrating the continuity equation in the SPH method, and is shown on the left side of Fig. 10. The other is summation density, which is the sum of kernel function similar to the particle number density in the MPFI and MPS methods, and it is a dependent variable of the particle position. Since the SPH method did not calculate the summation density, it is calculated afterward and shown on the right side of Fig. 10. These two density evaluations are ideally equivalent, but they were different from each other. The compression at the center (A in Fig. 3) evaluated by the summation density was less than 1%, but the value evaluated by the continuum density was larger than 3%. This contradiction in density estimation is thought to be the reason for the overestimation of the pressure shown in Fig. 7. When only the continuum density used in the SPH calculation increases and the summation density, which is interpreted as the reciprocal of the particle volume, is kept constant, the mass of the particle, which is represented by the product of the density and the volume, will increase. Since the number of particles in the system is constant, the total mass in the system will increase when the mass of each particle increases. Therefore, the difference between the two densities in Fig. 10 implies that the total mass in the system increased in the SPH calculation. This is against the mass conservation law, and it is a possible reason for the pressure overestimation in Fig. 7. The other SPH methods using the continuum density might suffer a similar problem with respect to the contradictory density evaluations.
The history of kinetic energy is shown in Fig. 11. The kinetic energy was almost zero throughout the MPFI calculation because there was almost no particle movement. In the MPS calculation, the kinetic energy increased and was maintained at a certain level. It implies that the particles were kept moving in the MPS calculation. This unintended particle movement is not favorable, not only because it causes pressure fluctuation but also because it might cause system instability like particle scattering. However, the usage of smaller time step size, Δt = 0.001 s, could not suppress this fluctuation, or rather resulted in larger fluctuation as shown in Fig. 11. This fluctuation was against the second law of thermodynamics, because it gave rise to the mechanical energy increase. One possible reason for this instability in the MPS calculation is the formulation of the discretized pressure Poisson's equation (PPE), whose source term diverges in the limit of Δt → 0. The other PPE-based incompressible particle methods [16][17][18][19], which adopted the source term having this property, may also have such undesired dependency on the time step width, Δt. In the SPH calculation, after the kinetic energy arose at the beginning of the calculation, it gradually decayed with the periodic oscillation as the energy in the system was dissipated by the artificial viscosity.
The calculation time of the hydrostatic pressure calculation using the MPFI, MPS and SPH methods is 175 s, 67 s and 76 s, respectively. The calculation time was measured by conducting a single thread calculation with Intel Core i7-8656U 1.9 GHz. The MPFI method took a longer time than the MPS and SPH methods because it solves the larger matrix system. The size of matrix equation was three times larger than that in the MPS method because not only the pressure but also x and y elements of the velocity vector were the unknowns to be calculated. Following the theory of computational complexity of matrix solver, the calculation time could be more than three times longer than that of the MPS method, but it was shorter. This is because the particles were almost stopped in this case, and not so many solver iterations were needed. It is thought that more calculation time is needed in the case where particles move, and in fact, it took 12 min 20 s for the case with initial fluctuation, which is going to be shown later.
Although the calculation time is much longer with the MPFI method, it could calculate the static balance of the incompressible hydrostatic pressure problem much better than the MPS and SPH methods. The MPS method suffered the unintended fluctuation, and the SPH method faced the contradictory density evaluation. They were against fundamental physics, and such an unphysical feature may cause system instability or incorrect evaluation. In contrast, the MPFI method did not show such unphysical behavior and the static balance was calculated very well. This tendency was the same when various particle sizes l 0 = 0.04, 0.02, 0.01 m and time step widths Δt = 0.008, 0.004, 0.002, 0.001 s were applied. In the MPFI calculation, the particles will not move once the balance is reached, and even when it deviates from the balance, the system will return to the balance. Therefore, it is advantageous in calculating static balance problems.
However, the initial lattice distribution seldom appears in a problem where the particles dynamically move around, and it is not certain whether the system can approach the static balance when the lattice distribution is destroyed. Therefore, the static pressure problem is also calculated with a given random initial velocity to see the response against fluctuation. The x and y elements of the initial velocity are set following the Gaussian distribution, N(0,0.06), where the variance, 0.06, was chosen so as not to affect the macroscopic behavior. The calculation is conducted with varying the particle size l 0 = 0.04, 0.02, 0.01 m and the time step width Δt = 0.008, 0.004, 0.002, 0.001 s, respectively. The history of the pressure at the center (point A) in the base case is shown in Fig. 12. The pressure fluctuated due to the particle motion, and the fluctuation got smaller as the system energy was dissipated. The kinetic energy of the system is shown in Figs. 13 and 14. They imply that the MPFI method has an energy dissipating feature, and the dissipation is larger when the particle size is smaller and the time step width is larger. In the limit of Δt → 0, the density control equation [Eq. (16)] is equivalent to the constant density constraint [Eq. (15)]. However, in practical calculations, the time step width has to be finite and the particle number density deviates due to the time discretization error. When the time step width is larger, the density deviation will be also larger due to the larger time discretization error, and it causes the larger energy dissipation. This is why the energy dissipation was larger when the larger time step width was adopted as shown in Fig. 14.

Droplet extension calculation
The droplet extension is one of the benchmarks for an incompressible free surface problem, which is adopted by Monaghan [1] and Ellero et al. [12]. The initial radius of the droplet is 1.0 m, and the initial velocity was given by where x and y are the coordinates whose origin is set at the center of the droplet. The calculation parameters used in the droplet extension calculations are shown in Table 2. The snapshots obtained by the calculation with the particle spacing, l 0 = 0.02 m, are shown in Fig. 15. The deformation of the droplet is calculated with keeping the particle number density almost constant. Theoretically, the length of the short axis, a, and long axis, b, follows  where ω is the initial value of ab. Equation (37) is solved numerically, and the solution was compared with the MPFI results in Fig. 16. The results of the MPFI method agreed well with the solution, and the sensitivity with respect to particle size was small. The history of the total mechanical (37) energy is shown in Fig. 17. The small energy decrease was observed in this calculation.

Standing wave calculation
The standing wave calculation is one of the benchmarks to check the property with respect to wave propagation, which was adopted by Suzuki et al. [13]. The initial particle distribution is shown in Fig. 18. The wavelength λ wave was set the same as the pool width, and the velocity of all the particles was set zero at the initial state. The initial surface displacement was given by where A = 0.1 h and h = 1.0 m were adopted in this study.
The parameters used in the calculation are shown in Table 3. The snapshots obtained from the MPFI calculation with l 0 = 0.02 m are shown in Fig. 19. The elevation of the surface at the center of the pool is plotted in Fig. 20, where the origin in the vertical direction is set at the initial average elevation of the fluid surface. In Fig. 20, the calculation results are compared with the first-and second-order approximations given by   Compared to the theoretical approximations, the wave decay was larger in the MPFI calculation. This is because the MPFI method has an energy dissipating feature. Suzuki et al. [13] reported that the wave decay occurred even when the mechanical energy was conserved because of the local randomness of the velocity distribution. In the MPFI calculation, the local randomness was dissipated and the wave decay was accelerated compared with the calculation conducted by Suzuki et al. [13]. The mechanical energy in the MPFI calculations is shown in Fig. 21. The total mechanical energy monotonically decreased.

Dam break calculation
One of the major advantages of particle methods is that they can easily treat the complex motions of free surface flows including coalescence and breakup. To confirm the ability of the MPFI method in this aspect, the dam break calculation was conducted, and the results were compared with the experiment [33]. The dam break calculation with the same geometry was also conducted in the previous studies, e.g., Asai et al. [24]. The initial particle arrangement is shown in Fig. 22. The initial velocities of the particles are set zero. The calculation parameters are shown in Table 4.
The fluid motion obtained by the MPFI calculation is shown in Fig. 23. The fluid motion in the experiment, which  was provided in the literature [33], was well simulated by the calculation mostly with keeping the density error within 4%. The pressure distribution is shown in Fig. 24. The smooth pressure field can be obtained via Eq. (34) even in the flow with dynamic motion. The comparison with respect to dam front position x is shown in Fig. 25. In the figure, the nondimensional position, x/h, and time, t(g/h) 1/2 , are used, where h is the representative length shown in Fig. 22. The results mostly agreed with the experiment [33].
The pressure history at the wall (A in Fig. 22) is plotted in Fig. 26 and compared with the experimental results [33]. For the pressure evaluation on the wall, the test period was set in applying the virial theorem, and Eq. (35) is used instead of Eq. (34). In specific, the 0.04-s time averaging was applied to remove the fluctuation due to the particle movement. Compared to the experimental result [33], the first pressure peak, which appears when the fluid hit the wall, was underestimated in the calculations, but the second peak, which appears when the fluid fell down again, was well reproduced. The possible reason for the first underestimation is the sparse particle distribution because of the particle reflection at the wall. Just after the fluid hit the wall, some particles detached from the main lump in the calculation, but the fluid in the experiment stayed in a lump at that moment. Even with such deviation, the pressure evaluation based on the virial theorem performed mostly well also in the flow including dynamic motion.
The history of the mechanical energy is shown in Fig. 27. The large energy decay was observed when the fluid collides with each other. To reproduce the fluid motion, the energy decay is also to be calculated properly. Therefore, it is thought that the adequate energy decay was calculated by the MPFI method. However, the energy decay in the MPFI method is not intentionally introduced and is only calculated as a result of solving the incompressible control equation. Hence, the agreement in fluid motion here is just accidental from our current knowledge, and the proper evaluation of the energy dissipation is left as an issue for the future work.

Conclusion
A new incompressible particle method with physical consistency, the moving particle full-implicit (MPFI) method, was developed. The interaction forces were composed of the conservative force and dissipative force, which could be fit into the framework of the analytical mechanics [11]. It implies that the fundamental laws of physics, such as the second law of thermodynamics, were basically satisfied in the proposed system. Instead of directly posing the constant density condition, the control equation which constraints the velocity divergence was adopted to keep the fluid incompressible. The control equation enables to avoid nonlinear equation as in the previous studies [12,13] with introducing a certain energy dissipation. The equation was solved in a full-implicit algorithm, where the velocity and pressure were calculated implicitly at the same time. The implicit calculation only needed to solve the linear symmetric matrix equation once in each time step. Furthermore, a new pressure evaluation procedure based on the viral theorem [32] was proposed for the physically consistent system.
Using the MPFI method, static pressure, standing wave, droplet extension and dam break calculations were conducted.  In the static pressure calculation, the results were compared with those of the MPS and SPH methods. All the three methods could obtain smooth pressure fields. However, in the MPS calculation, unintended fluctuation against the second law of thermodynamics was observed, whereas in the SPH calculation, the density evaluation against the mass conservation law was observed and it resulted in the overestimation of the pressure. In contrast, the MPFI could calculate the hydrostatic pressure balance very well although the calculation time was much larger than the other two methods. In addition, the energy dissipating property of the MPFI method depending on the particle size and the time step width was studied.
In the droplet extension calculation, the change in the length of the short and long axes agreed well with the theoretical solution, and the results depended little on the particle size. It implies that the MPFI method could well capture the incompressible free surface flow in this problem.
In the standing wave calculation, the wave was compared with the theoretical solution. In the MPFI calculation, the large wave decay was observed compared to the theoretical approximation because of the energy dissipative property of the MPFI method.
In the dam break calculation, the results were compared with the experiment [33]. The prediction of the dam front position mostly agreed with the experimental results. The pressure at the wall was also compared with the experiment. It is confirmed that the pressure evaluation based on the virial theorem [32] was valid even under the flow with dynamic motion.    Fig. 22) in the dam break calculations using the particle size of l 0 = 0.04, 0.02, 0.01 m Overall, the MPFI method, satisfying the fundamental laws of physics, had an advantage in calculating hydrostatic balance, and the reasonable results could be obtained in other problems. Compared to the MPS and SPH methods, the MPFI method is thought to have an advantage in a problem where strict incompressibility and static balance are required, and the physical consistency may help conducting simulation without tuning artificial relaxation parameters to cope with unphysical behaviors. Moreover, it may have superiority in the high-viscous fluid calculation because the implicit velocity calculation is free from the restriction with respect to the diffusion number.