Applying surface tension as pressure boundary condition in free surface flow analysis by moving particle simulation method

A model that introduces surface tension as a pressure boundary condition, named the surface tension as pressure (STP) model, was developed for free surface flow analyses by the moving particle simulation (MPS) method. The STP model assigns to surface particles the liquid pressure of Laplace’s formula. The model is an alternative to previous models that apply surface tension as volume force such as the continuum surface force model. Problems that appeared when using the volume force models, such as the dependencies of calculation results on particle resolution and pressure gradient accuracy, were solved by using the STP model. Calculations predicted the theoretical values of the internal pressure of a 3D spherical droplet and the oscillation period of a 2D elliptic droplet over a wide range of surface tension coefficients and droplet sizes with errors less than 10%. Since the STP model is easy to implement, does not increase computation cost from previous models, and does not require surface reconstruction or additional marker particles, the model is suitable for practical and large-scale free surface flow problems that involve violent deformation of the liquid surface such as liquid atomization.

the liquid mass conservation is guaranteed, and simulations are executable at low calculation cost. Therefore, introducing surface tension into free surface flow analyses by the MPS method both accurately and simply is an important issue.
From the viewpoint of how surface tension is applied to the momentum equation, we sorted surface tension models for computational fluid dynamics into two groups: the volume force form and the boundary condition form (Table 1). A majority of studies using particle methods apply surface tension in the volume force form. On the other hand, this paper proposes a model that is in the boundary condition form, which is a minor approach.
The continuum surface force (CSF) model is in the volume force form and is one of the most popular surface tension models in both Euler grid methods and particle methods. Theoretically, for inviscid fluids and when the surface tension coefficient is constant, the pressure difference across the gasliquid interface should satisfy Laplace's formula [8]: Here, p l , p g are the pressures of the liquid and gas phases, σ is the surface tension coefficient, and κ is the interface curvature (positive when the center of curvature is on the liquid side). In Euler grid methods, applying a pressure boundary condition to an interface requires interface reconstruction. The reconstruction is difficult especially when the interface is topologically complex. Therefore, Brackbill et al. [9] decided to distribute the following volume force to grid cells close to the gas-liquid interface.
f SV = −σ κδn (2) δ is a delta function, andn is a unit surface normal vector which points to the gas side. In this way, Laplace's formula can be indirectly satisfied, just by adding the force f SV to the source term of a flow solver. The CSF model was introduced to the SPH method by Morris [10] and to the MPS method by Nomura et al. [11]. The force f SV was simply applied to particles in the vicinity of the interface in a two-phase analysis. Today, the CSF model is commonly used in free surface flow analyses as well [2,12].
Other surface tension models used in particle methods can be sorted into two groups [13]. The cohesive pressure approach considers that the negative pressure yielded from the equation of state functions as surface tension [13,14]. The potential force approach applies a pairwise force that mocks intermolecular force [15,16]. While the theoretical basis are different, all three approaches share the idea of applying surface tension as volume force.
The previous models have been validated through calculations of benchmark problems such as the internal pressure of a stationary droplet [10,13,15] and the oscillation period of an elliptic/square droplet [10,11,13,15,16]. However, we found that the results might differ according to calculation conditions such as particle resolution. Moreover, an unnatural pressure distribution appeared near the surface, as pointed out by other researchers [17]. As we will discuss in Sect. 4, we suspect that the problems stem from applying surface tension as volume force. Let us recall here that the CSF model was developed for Euler grid methods in order to circumvent the complex interface reconstruction procedure. In free surface flow analyses using particle methods, applying a pressure boundary condition to surface particles is easy. A surface tension model does not necessarily have to be in the volume force form.
The boundary condition form refers to applying the pressure jump condition denoted by Laplace's formula directly to the gas-liquid interface. The approach is already a common measure in an Euler grid method called the ghost fluid method (GFM) [18]. Liu et al. [19] enabled solving the pressure Poisson equation (PPE) for a discontinuous pressure field, and Desjardins et al. [20] applied the method to a gasliquid two-phase flow calculation. Sussman et al. first used the CSF model [21] but switched to using an approach similar to the GFM [22]. If a flow solver can treat the interface as a discontinuity, applying surface tension as a boundary condition is the preferred choice today.
There are a few models for particle methods that are in the boundary condition form. Fürstenau et al. [23] (labeled FST in Table 1) basically used the CSF model, a volume force model, in their incompressible SPH calculations. Meanwhile, when solving the PPE, the boundary condition at the surface was determined by Laplace's formula. The aim was to give a well-posed Dirichlet boundary condition to the PPE, and this successfully improved numerical stability. However, as we will discuss in Sect. 2.3, we anticipate that applying volume force becomes unnecessary if another pressure gradient model is chosen. The moving surface mesh method developed for the MPS method by Matsunaga et al. [17,24] applies surface tension solely in the boundary condition form. Marker particles are placed outside liquid particles. At each marker particle, surface tension is first evaluated as force, but the modeling arrives at a boundary condition for the PPE.
Calculations of benchmark problems were all highly accurate and stable. A possible drawback is that handling the marker particles requires additional computation time and memory. In some industrial applications such as liquid atomization, the liquid surface area is extremely large and grows rapidly, and breakup and coalescence of droplets occur violently. A simpler surface tension model would be beneficial for such practical simulations.
In this paper, we propose the surface tension as pressure (STP) model for free surface flow analyses by the MPS method. The model directly assigns the liquid side pressure of Laplace's formula to surface particles. No volume force is applied. The STP model is simple and requires low cost. Details of the STP model and how we implemented it into MPS calculations will be discussed in Sect. 2. After briefly describing in Sect. 3 the numerical methods and conditions, we will present results and discussion of two test problems, a stationary droplet and an oscillating droplet, in Sect. 4. A conclusion will be given in Sect. 5.

Basic idea
The basic idea of the STP model is presented here using Fig. 1. The figure is a one-dimensional schematic of the pressure distribution across a gas-liquid interface. Following the CSF model, we assumed that the viscosity of the gas and liquid phases is negligible, and the surface tension coefficient is constant. Therefore, the pressure difference between the gas and liquid phases should satisfy Laplace's formula, Eq. (1). Ideally (Fig. 1a); the width of the gas-liquid interface is infinitesimally small, so the pressure should change discontinuously across the interface. Volume force form CSF model [9] CSF model [10,11] (FST a : prediction step) Cohesive pressure approach [13,14] Potential force approach [15,16] Boundary condition form Applying pressure boundary condition in the GFM [19,22] FST a : PPE step Moving surface mesh method [17,24] STP model of the present study

(a) (b) (c)
In conventional calculations using a volume force model (Fig. 1b), a surface particle's pressure is fixed to the gas atmospheric pressure, which is usually zero. This is a common method in both the original (semi-implicit) MPS method [1] and the explicit MPS method [25]. By doing so, the surface particle is assumed to be on the gas side of the gas-liquid transition region in terms of pressure, although it is a liquid particle. The volume force f SV pushes the surface particle to the liquid side. Due to the increase in particle density, the pressure of internal particles rises, and the surface particle settles at a position where the surface tension force and the pressure gradient force balance. The resulting pressure difference between the surface and internal particles is expected to equal Laplace's formula. However, this is not necessarily guaranteed. The pressure gradient model's accuracy can affect the evaluated pressure gradient and thus the internal pressure. In addition, the steep pressure change at the surface can be a source of numerical errors when evaluating the pressure gradient.
On the other hand, the STP model directly gives the liquid pressure of Laplace's formula, p l = p g + σ κ, to a surface particle (Fig. 1c). If that pressure is higher than pressures of internal particles, the pressure gradient force −∇ p will pull the surface particle toward the internal particles. The increased particle density leads to an increase in internal pressure. If the surface pressure is lower than the internal pressure, the surface pressure is pushed outwards, and the internal pressure decreases. At the time when the calculation reaches a steady state, the pressure gradient of each particle approaches zero, so the pressures of internal particles will be equal to that of the surface particle. As a result, the pressure difference between the gas and liquid phases satisfies Laplace's formula. The pressure gradient within the liquid particles will not be a differential across a steep pressure change, so problems caused by numerical errors are less likely to occur. The STP model is a natural model that gives the pressure of the liquid phase to a liquid particle.

Implementation to MPS method
The implementation of the STP model to the MPS method is simple. At each time step, first, surface particles are detected. Next, the surface normal vector and the surface curvature are evaluated at each surface particle. Finally, the liquid pressure in Eq. (1), is given to each surface particle i. For free surface flow analyses, the gas pressure p g is the atmospheric pressure ( p ∞ = 0) throughout the region. In this study, the STP model was introduced to the explicit MPS method. Details of the surface pressure evaluation adopted in this study are described in the following paragraphs. Particles that satisfy Eqs. (4) and (5) were detected to be surface particles: Here, n i is the particle density evaluated by and n 0 is the initial particle density. n i is the mean direction to neighbor particles estimated by r i j is the relative position between particles i and j, r i j is the length of that vector (| r i j |), and w(r ) is a weight function. Constants β 1 and β 2 were set to 0.97 and 0.2, respectively. Equation (4) is the standard surface particle criterion of the MPS method [1], and Eq. (5) is an additional criterion similar to that of Khayyer et al. [26]. The surface normal vector was simply evaluated using n i : The surface curvature κ i was evaluated as the divergence of the surface normal vector. We modified the MPS method's standard divergence model to evaluate the divergence over the liquid surface: The summation of the right-hand side is executed over neighbor surface particles ( j ∈ surface). Therefore, the summation should be normalized by the surface particle density: The dimension number in Eq. (9) should be d = 2 in a threedimensional calculation because the divergence is evaluated on a two-dimensional liquid surface. For the same reason, d is 1 in a two-dimensional calculation. Improving the curvature evaluation has been an important issue in studies using the CSF model because the error in the curvature is a common source of numerical instabilities. Some studies used the SPH/MPS methods' standard divergence models [10,27], while others proposed improved methods [12,[28][29][30]. We used Eqs. (8)(9)(10) to simplify the calculations. Figure 2 shows a result of the surface particle detection, the normal vector evaluation, and the curvature evaluation in a typical stationary droplet calculation. Details of the MPS calculation will be introduced in Sect. 3. Surface particles were successfully detected, and the evaluated curvatures were close to the theoretical value. The relative error between the average curvature over all surface particles and the actual curvature of a sphere surface was within 3% throughout the calculation time.

Selection of pressure gradient model
The pressure gradient force at the surface plays an important role in adjusting the pressure of internal particles, as discussed in Sect. 2.1. A surface particle should be pulled toward internal particles if the surface pressure is higher than the internal pressure, and vice versa. However, some pressure gradient models do not meet this requirement. Therefore, a pressure gradient model should be properly selected.
There are several pressure gradient models for the MPS method. The MPS standard model [1] is written in the following formula: Commonly, in the MPS method, p i is replaced to the minimum pressure inside the effective radius p min i [1]: Other examples are the CMPS method developed by Khayyer and Gotoh [31] and the model suggested by Oochi et al. [25]. Equation (11) is purely a weighted average of pressure gradients between particle pairs. Meanwhile, Eq. (13) and the other models implicitly contain an additional term so that the pressure gradient force between particle pairs always become repulsive. Using repulsive force models is advantageous in improving numerical stability [32]. A similar discussion exists in the SPH method. According to Price [33], a pressure gradient evaluation similar to Eq. (11) (Eq. (82) of Ref. [33]) is accurate, while using an evaluation that yields repulsive force between particles (Eq. (31) of Ref. [33]) is advantageous because the particles will be automatically rearranged to a regular distribution.
The pressure gradient force always being repulsive is a problem when using the STP model. This is because a surface particle will never be pulled toward internal particles. Moreover, the surface particle will be applied excess force outwards the liquid bulk, which leads to an expansion of particles. Therefore, a model that applies attractive force toward low pressure particles, for example Eq. (11), must be selected for pressure gradient estimation at surface particles.

Velocity correction for surface particles
An artificial velocity correction was introduced to surface particles. If a particle's velocity is pointed outwards and too large, the velocity's normal component is corrected to that of neighbor particles. The aim is to maintain particle density at the surface by preventing surface particles from spreading out.
Such a correction is unnecessary in conventional calculations because the volume forces constantly pull surface particles to the liquid bulk (Fig. 3a). In contrast, in calculations using the STP model, the force that pulls a particle to the liquid bulk is the pressure gradient term (Fig. 3b). When a particle has started to leave the liquid bulk, the pressure inside the liquid bulk becomes low due to the decreased particle density, and the surface pressure is increased by the STP model according to the increased curvature. As a result, the pressure gradient force acts to the liquid bulk direction. The problem is that the pressure gradient force becomes weaker as the distance between particles increases, mainly because of the weight function. If a surface particle gains large velocity in the outwards direction by some reason, the pressure gradient force can be insufficient to pull back the particle.
As a countermeasure to that situation, the velocity correction was introduced. The correction prevents surface particles from traveling too far so that the pressure gradient force works properly. The correction itself does not function as surface tension because it does not actively pull back particles and is not dependent on surface curvature. It just helps keeping particles together so that the STP model works.
The velocity correction was applied to surface particles whose particle density was lower than the initial particle density at the surface n 0,half , and whose relative velocity to neighbor particles was pointing outwards, n 0,half is the density evaluated at surface particles in an initial particle arrangement (Fig. 4a) and was 0.71n 0 for the calculations in this paper. The mean velocity of neighbor particles u ave i was evaluated by The velocities of particles that satisfy conditions of Eqs. (14) and (15) were corrected at each time step by the following formula, so that their normal components equal the normal component of u ave i (Fig. 4b).
In the first calculation case of this paper, a stationary spherical droplet, the velocity correction was applied to 9% of surface particles at each timestep in average. Mean velocity of neighbor particles Surface particle velocity

Corrected velocity
Initial particle density at the surface Initial particle density : surface particle

Overview
The motions of stationary and oscillating droplets were solved by the explicit MPS method developed by Oochi et al. [25]. The governing equation is the incompressible Navier-Stokes equation: μ is the viscosity coefficient. Figure 5 is a flow chart of the MPS calculations of this study. Details of the numerical method and conditions will be presented in the following subsections.
Calculations using the CSF model and the STP model were compared. When using the CSF model, surface tension was introduced by applying volume force to surface particles in the velocity/position prediction step. As pressures of surface particles, the gas atmospheric pressure ( p g = p ∞ = 0) was given. In calculations using the STP model, no volume force was applied, and pressures evaluated by the STP model were given to surface particles. Therefore, surface tension was considered in the velocity/position correction step.

Surface tension evaluation by the CSF model
The volume force of the CSF model is evaluated by Eq. (2) and was applied to surface particles. Detection of surface particles and evaluation of normal vector and curvature were conducted by the methods described in Sect. 2. Therefore, the intensity of surface tension, σ κ, was evaluated in the same way as in the STP calculations, while the form of applying surface tension, i.e., the volume force form or the boundary condition form, was different.
The CSF model was applied to a single layer of particles. The delta function δ was approximated as 1/l 0 , where l 0 is the initial particle spacing. By doing so, a line integral of the delta function across the liquid surface equals one. Applying surface tension force to a single layer of particles is seen in several studies [11,27,29], while in other studies, the force is spread to several particle layers because the delta function is evaluated by a smoothed color function [10]. We tested applying force to single and several particle layers and confirmed that the main problem of the volume force models, that the results depend on particle resolution and pressure gradient accuracy, appear in both cases. Results of the former case are presented in this paper.

Pressure evaluation for internal particles
Pressure of internal particles was evaluated by the formula of the explicit MPS method: The atmospheric pressure was set to zero. c s is a virtual sound speed that does not necessarily have to be the real sound speed of the fluid. There is a restriction that the virtual Mach number | u|/c s should not exceed 0.2 [34]. The virtual sound speed was set to 10 m/s, which meets the restriction.
In the MPS method, pressure of internal particles near the surface tends to be unnaturally low because the particles' densities n i are lower than n 0 . To remove such unphysical phenomena, pressure below the atmospheric pressure (negative pressure) is modified to the atmospheric pressure (zero) [25,35]. However, the internal pressure being negative at around a concave surface, where the surface pressure is negative because κ < 0, is natural. Therefore, the pressure modification was altered as follows. If there appeared a concave surface somewhere in the calculation region, the lower limit of pressure was set to the minimum pressure of all surface particles. The new pressure modification is expressed as follows: Elsewhere in this paper, the modified pressure p mod i is simply referred to as p i to avoid equations being lengthy.

Pressure gradient evaluation
Two sets of pressure gradient models with different accuracies were tested. The reason is to test our hypothesis that calculation results depend on the accuracy of pressure gradient evaluation when using a volume force model (Sect. 2.1). As discussed in Sect. 2.3, a pressure gradient model that applies attractive force toward a particle with lower pressure should be selected for surface particles. For internal particles, selecting a model that enhances numerical stability is a better choice. Therefore, the following models were chosen as the first set and is called the MPS standard model in this paper: The MPS standard model is said to be zeroth-order accurate when particles are randomly distributed [24]. The second set of models is the gradient correction (GC) model by Khayyer and Gotoh [36]: where is a correction matrix that guaranties first-order accuracy. Figure 6 shows a comparison between pressure gradients evaluated by the MPS standard model and the GC model. Pressure distribution was artificially determined so that the gradient is one in the horizontal direction. The MPS standard model underestimated gradient at surface particles. This is because the weight function is standardized by the initial particle density, n 0 , which is larger than the actual summation of the weight function. On the other hand, the GC model was accurate at all particles.

Other information
Physical properties of water, ρ = 1000 kg/m 3 , μ = 1.0 × 10 −3 Pa · s, and σ = 7.2 × 10 −2 N/m, were used except when specified. The viscous term was evaluated by the MPS method's standard model [1]. The following weight functions were used: The first and second weight functions are the ones recommended by Yamada et al. [37]. The particle densities n i , n 0 , n surf i , n 0,half , the curvature κ i , and the mean velocity of neighbor particles u ave i were evaluated using the first one, while the pressure gradient was evaluated using the second one. The third one was used for Eq. (7) so that the surface normal vector would not be too sensitive to the nearest neighbor particle. The effective radius r e was set to 3.1l 0 in all the summations. The dynamic stabilization (DS) scheme of Tsuruta et al. [38] was used to stabilize calculations. The DS scheme applies the minimum force needed to dissolve overlapping of particles. The velocity correction for surface

Internal pressure of stationary droplet
First, a three-dimensional spherical droplet was calculated to investigate the pressure distribution within the liquid phase. Theoretically, the pressure should be uniform at a steady state, equal to the Laplace pressure, where R 0 is the droplet radius. The droplet was initially placed in a spherical shape. Figure 7 shows pressure distribution inside the calculated droplets. The droplet diameter D 0 (= 2R 0 ) was 5 mm, the particle resolution was D 0 /l 0 = 32, and pressure gradient was evaluated by the GC model. Pressure of each particle was time-averaged over 100 timesteps in order to remove the pressure oscillation typically seen in explicit MPS calculations [34] and make the figures readable. When surface tension was applied as volume force by the CSF model (Fig. 7a), pressure of internal particles far from the surface was typically lower than the theoretical Laplace pressure (labeled A in the figure). An unnatural pressure distribution was observed near the liquid surface, as we mentioned in Sect. 1. A pressure spike was seen at particles of the second layer from the surface (B), while the pressure of the third layer was lower than neighboring particles (C). As another problem, there was an unnatural spatial gap between the second and the third layers (D). In contrast, when surface tension was applied as a pressure boundary condition by the STP model (Fig. 7b), the pressure distribution inside the droplet was closer to uniform, and the value was near the theoretical Laplace pressure (E). Several particles near the surface showed high or low pressure (F), although the number of those particles was far less than in the CSF case. The particles were evenly distributed.
The mechanism for the pressure spike and the spatial gap seen in the CSF case can be explained as follows. Suppose that particles are initially spaced at even intervals (Fig. 8a  left). When the surface tension force f SV is applied to a surface particle, labeled No. 1 in the figure, the particle moves toward particle No. 2. Consequently, the particle density rises at around particle No. 2, and its pressure rises (Fig. 8a center). Particle No. 1 stops at a position where the surface tension force and the pressure gradient term −∇ p balance. Next, particle No. 3 moves toward the inside of the liquid, because of the pressure difference between particle No. 2 and inner particles. The other particles also move around until the pressure gradient at their position becomes zero (Fig. 8a right). Particle No. 2 is pushed outwards by the pressure difference between particles No. 1 and 3, but is pushed back by the dynamic stabilization (DS) model to prevent particle No. 2 from getting too close to particle No. 1. At this moment, the pressure of particle No. 3 does not reach the pressure of particle No. 2. This is because the pressure gradient at particle No. 3 becomes zero when the average pressure of the outer side (particles No. 1 and 2) and the inner side (particles No. 4, 5, etc.) balance. As a result, only the pressure of particle No. 2 becomes relatively high. Since particle No. 2 is continuously pushed toward particle No. 1 by the pressure gradient force, a gap appears between particles No. 2 and 3. From the discussion above, it is clear that the appearance of the pressure spike and the spatial gap between particles is inevitable as long as surface tension is applied as volume force.
On the other hand, when the STP model is used, the pressure jump due to surface tension is directly given to surface particles (Fig. 8b left). The pressure jump equals the Laplace pressure when the droplet is spherical. After a transient state (Fig. 8b center), the calculation reaches a steady state when the pressure gradient of all particles approaches zero (Fig. 8b  right). Consequently, the pressure of all particles becomes close to the Laplace pressure. The pressure fluctuation seen at the surface (labeled F in Fig. 7) is a result of particle oscillation at the surface, which we will discuss later using Fig. 10. Figure 9 shows the internal pressure of a droplet at several calculation conditions. The internal pressure was evaluated as an average of pressure over all internal particles through timesteps 800 to 1000. The droplet diameter was 5 mm. When the CSF model was used, the internal pressure varied depending on the particle resolution and the accuracy of the pressure gradient model. The internal pressure was high in a low particle resolution condition probably because the number of particles in the second layer from the surface, which possess high pressure, was relatively large. The internal pressure was higher when the MPS standard pressure gradient model was used. This was because the MPS standard model underestimates pressure gradient at the surface, as shown in Fig. 6.
Since the pressure gradient force that acts on a surface particle to the outward direction was underestimated, the surface particle was pushed further to the inward direction, resulting in a high internal pressure. When the GC model was used, the internal pressure was half the amount of the theoretical Laplace pressure. The results demonstrate that even if an accurate pressure gradient model is used, the resulting pressure jump at the liquid surface does not necessarily satisfy Laplace's formula. Meanwhile, when the STP model was used, the dependency of the internal pressure on the calculation condition was relatively small. This was because the process of the internal pressure approaching the Laplace pressure, which is explained using Fig. 8c, functions regardless of the particle resolution and whether the pressure gradient is underestimated or not. The relative error to the theoretical Laplace pressure was within 10%. The results demonstrated that using the STP model is more accurate on reproducing the Laplace pressure compared to the CSF model. All internal pressures of the STP cases in Fig. 9 were lower than the Laplace pressure. The oscillation of surface particles is probably the cause. Figure 10a shows a schematic of the oscillation. If a surface particle is far from internal particles, the internal pressure is low, and the surface particle is pulled to the inward direction. Oppositely, when the surface particle is close, the high pressure of the internal particles pushes the surface particle outwards. As a result, every surface particle continuously oscillates around its neutral position, with the pressure gradient force acting as a restoring force. The problem is that the restoring force is weak when the surface particle is distant from the internal particles. Figure 10b is an example of the restoring force as a function of the surface particle's position. Here, an ideal situation was set to evaluate the restoring force. Particles were evenly placed at an interval of unit length and one surface particle was moved in the z direction from its neutral position z = z 0 . The pressure gradient was evaluated by the MPS standard gradient model, and the effective radius was 2.1. Since the restoring force is small when z > z 0 , a surface particle will stay in the z > z 0 side for a longer time during the oscillation. Consequently, the time-averaged internal pressure will be lower than the neutral value. This is probably the reason why the internal pressure, which was the time-averaged pressure of internal particles, was underestimated in the STP cases. Figure 11 shows the internal pressure when the surface tension coefficient σ and the droplet diameter D 0 = 2R 0 were changed. The STP model was used, the particle resolution was D 0 /l 0 = 32, and the pressure gradient was evaluated by the GC model. The calculated pressure clearly followed the trend described by Eq. (26). In all cases the relative error to the Laplace pressure was within 10%. Therefore, the STP model can be applied to static problems of liquid other than water and droplets with various sizes. The applicability of the model to dynamic problems will be discussed in the next subsection.

Oscillation period of elliptic droplet
As a second demonstration of the STP model, the oscillation of a two-dimensional elliptic droplet was calculated. The oscillation period was theoretically derived by Lamb [39], which is written as The droplet was initially placed as an ellipse whose ratio of the major and minor axes was 2. Figure 12 shows visualized images of an oscillating droplet, and Fig. 13 shows the time change of the droplet's width. The mean diameter D 0 (= 2R 0 ) was 5 mm, the particle resolution was D 0 /l 0 = 32, and pressure gradient was evaluated by the GC model. Qualitatively, the droplet's oscillation was reproduced in both the CSF and STP cases, while the oscillation period and the amplitude differed among the two cases. Figure 14 shows oscillation periods at several calculation conditions. The oscillation period T was evaluated as an average of three wavelengths, as indicated in Fig. 13. The oscillation period was more accurately reproduced when the STP model was used. In the CSF cases, the oscillation period had a large dependence on the pressure gradient accuracy. When using the GC model, the oscillation period was more  Fig. 14 Oscillation period of a droplet initially placed in an elliptic shape at different calculation conditions than 40% longer than the theoretical value. This suggests that the surface tension was underestimated. It agrees with the underestimation of internal pressure of a stationary droplet in the same condition. On the other hand, in the STP cases, the relative error with the theoretical value was within 10%, independent of the pressure gradient accuracy. Figure 15 shows the dependency of the oscillating period on the surface tension coefficient and the droplet mean diameter. The STP model was used, the particle resolution was D 0 /l 0 = 32, and the pressure gradient was evaluated by the GC model. The results followed the theoretical correlation, Eq. (27), within an error of 4%. This shows that the STP model is applicable to dynamic problems over a wide range of surface tension coefficients and droplet sizes.
According to Fig. 13, the oscillation amplitude was attenuated by time, and the attenuation rate differed among calculation models. This paragraph will discuss possible causes of this tendency. Viscosity and artificial viscosity can be causes of the attenuation, as pointed out by previous researches [15]. The velocity correction explained in Sect. 2.4 might be another cause because the velocity correction acts as a limitation to the motion of surface particles. However, they are not causes of the difference among calculation models because the viscosity and the velocity correction were equally applied to all calculation cases. The attenuation rate was larger when using the MPS standard pressure gradient model and the STP model. The tendency suggests that the difference in attenuation depends on the pressure gradient force applied to surface particles. In a free surface flow, the liquid surface deforms when the pressure distribution inside the liquid phase changes. Suppose that the internal pressure was increased by p from a previous state. Figure 16 is a schematic of this situation. The change in force acting on a surface particle, f , is the difference that occurs to the pressure gradient force. For example, when using the MPS standard gradient model, the change in force is expressed as follows: As we have discussed earlier, the MPS standard model underestimates pressure gradient at surface particles. Similarly, the force f might have been underestimated, which resulted in preventing the liquid surface from deforming. The force f was probably greater in the CSF case compared to the STP case because the distance with the nearest particle (r 12 in Fig. 16) was shorter. As a result, the oscillation was more enhanced in the CSF case. Fig. 15 Oscillation period of an elliptic droplet when the surface tension coefficient (a) and the droplet mean diameter (b) were changed. The droplet mean diameter was 5 mm in Fig. (a), and the surface tension coefficient was 7.2 × 10 −2 N/m in Fig. (b). Solid lines indicate the theoretical correlation, Eq. (27)

Conclusion
This paper proposed the surface tension as pressure (STP) model for free surface flow analyses using the MPS method. The pressure boundary condition denoted by Laplace's formula is directly applied to the liquid surface. The model is an alternative to conventional models that apply volume force to surface particles. Two basic problems were calculated to test the model and compare it with the CSF model.
First, a three-dimensional stationary droplet was calculated. When using the CSF model, an unnatural pressure spike and a gap between particles was observed just under the surface. The internal pressure varied between cases of different particle resolution and pressure gradient accuracy. Meanwhile, when the STP model was used, the pressure inside a droplet was close to uniform, and its value was less dependent on the calculation condition. The theoretical relationship of the internal pressure with the surface tension coefficient and the droplet diameter was successfully reproduced. The results showed that the STP model functions well in static problems.
The second evaluation was calculating the oscillation of an elliptic droplet in two-dimensions. Again, results of the STP cases were less dependent on particle resolution and pressure gradient accuracy compared to the CSF results. The STP results followed a theoretical formula over a wide range of the surface tension coefficient and the droplet diameter. The results demonstrated that the STP model is applicable to dynamic problems as well.
Several shortcomings of the STP model were presented. A velocity correction was necessary to maintain particle density at the surface. The internal pressure of a stationary droplet was several percent lower than the theoretical value in all cases, which indicates that the actual surface tension applied to the liquid phase was slightly underestimated. The influence of these drawbacks should be assessed when using the model.
Nevertheless, the STP model was better than the CSF model in terms of independencies on numerical condition. The independency of calculation results on particle resolution is an advantage when simulating phenomena that contain droplets with various sizes, in other words, various particle resolutions. The GC model can be selected to improve the overall accuracy of a simulation. As a simple, low cost, easy to implement, and accurate surface tension model, we expect that the STP model will contribute to extending particle methods to simulations of practical and large-scale free surface flow problems. unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecomm ons.org/licenses/by/4.0/.