An Improved Dynamic Boundary Condition in SPH Method

Smoothed particle hydrodynamics (SPH) was initially proposed by Monaghan and Gingold [1], and by Lucy [2] in 1977. SPH has progressed tremendously due to intensive theoretical work and computational improvements, particularly since the 1990s. This numerical method has been shown to be robust and applicable to a wide variety of fields. The main advantages of SPH come from its pure Lagrangian nature. The method can efficiently deal with the existence of large voids without special treatment. Also, the method can handle large deformations and extreme fragmentation without additional complicated or rather unphysical approaches. From a numerical point of view, its Lagrangian nature brings simplicity and parallelism to its algorithm. In addition, multiphase media can be easily described if the relationships between the phases are determined. In fluid mechanics, fluids are described as discrete particles, and the governing equations are integrated into the Lagrangian form for these particles. The physical quantities of each particle can be obtained by interpolating the neighboring particles' relevant values through a kernel function. In some way, the role of kernel function is like a weight function used to obtain a weighted average of the values of the particles in the support domain. Computationally, the interpolation enables discrete information to be smoothed and continuous in the domain. SPH methods have been successfully used in astrophysical applications and hydrodynamic problems, such as free surface flow [3, 4], gravity current [5]. SPH method can also be applied to many other scientific and engineering applications, such as diesel engines, hydraulic engineering, and nanofluid flow. Boundary condition is one of the key issues to solve these problems which is also very important for other particle methods such as MPS, and considerable efforts have been devoted to it. Different types of boundary conditions can be used, such as ghost particles [6], repulsive particles [3], LUST [7], and so on. Among these boundary conditions, dynamic boundary condition (DBC) was first developed by Crespo, Gómez and Dalrymple [8] and has been widely used. In DBC, boundary particles share the same properties with fluid particles, and they follow the same equations of state and continuity when computing the density variation and the acceleration of fluid particles. The differences are that the boundary particles are not allowed to move, or they can only move according to some external input and the density of boundary particles are kept constant. In this way, all the variables associated with the particles can be calculated in the same manner with a considerable saving of computation time. Although DBC has been proved feasible in many cases, we found it would bring some error due to the asymmetry of the repulsive force field. When the fluid particles collide with the boundary, they may not be reflected specularly. Although this error may be tiny, it should be avoided because the interaction between the fluid and the boundary should be controllable and precise since how the error occurs and propagates is uncertain. In addition, in some simulations, we found a few fluid particles broke through the boundary. This phenomenon is unphysical, and it possibly indicated that the repulsive force was not strong enough. We noticed that local uniform stencil boundary condition (LUST) is presented to solve the two problems. The repulsive force and the symmetry have been improved to some extent in LUST. However, for some exceptional cases, such as waterdrop impact, these problems still exist. This paper presents an improved dynamic boundary condition (IDBC) in SPH methods with three improvements. Firstly, the boundary interaction is enhanced to decrease the number of boundary particles and prevent fluid particles from breaking through the boundary. Secondly, the shape of the support domain attached to the boundary is modified to keep the force field close to the boundary being symmetrical. Finally, a correction factor is introduced to guarantee the unity. Meanwhile, the original kernel is retained, and the algorithm keeps concise.


Introduction
Smoothed particle hydrodynamics (SPH) was initially proposed by Monaghan and Gingold [1], and by Lucy [2] in 1977. SPH has progressed tremendously due to intensive theoretical work and computational improvements, particularly since the 1990s. This numerical method has been shown to be robust and applicable to a wide variety of fields. The main advantages of SPH come from its pure Lagrangian nature. The method can efficiently deal with the existence of large voids without special treatment. Also, the method can handle large deformations and extreme fragmentation without additional complicated or rather unphysical approaches. From a numerical point of view, its Lagrangian nature brings simplicity and parallelism to its algorithm. In addition, multiphase media can be easily described if the relationships between the phases are determined.
In fluid mechanics, fluids are described as discrete particles, and the governing equations are integrated into the Lagrangian form for these particles. The physical quantities of each particle can be obtained by interpolating the neighboring particles' relevant values through a kernel function. In some way, the role of kernel function is like a weight function used to obtain a weighted average of the values of the particles in the support domain. Computationally, the interpolation enables discrete information to be smoothed and continuous in the domain.
SPH methods have been successfully used in astrophysical applications and hydrodynamic problems, such as free surface flow [3,4], gravity current [5]. SPH method can also be applied to many other scientific and engineering applications, such as diesel engines, hydraulic engineering, and nanofluid flow. Boundary condition is one of the key issues to solve these problems which is also very important for other particle methods such as MPS, and considerable efforts have been devoted to it. Different types of boundary conditions can be used, such as ghost particles [6], repulsive particles [3], LUST [7], and so on. Among these boundary conditions, dynamic boundary condition (DBC) was first developed by Crespo, Gómez and Dalrymple [8] and has been widely used. In DBC, boundary particles share the same properties with fluid particles, and they follow the same equations of state and continuity when computing the density variation and the acceleration of fluid particles. The differences are that the boundary particles are not allowed to move, or they can only move according to some external input and the density of boundary particles are kept constant. In this way, all the variables associated with the particles can be calculated in the same manner with a considerable saving of computation time.
Although DBC has been proved feasible in many cases, we found it would bring some error due to the asymmetry of the repulsive force field. When the fluid particles collide with the boundary, they may not be reflected specularly. Although this error may be tiny, it should be avoided because the interaction between the fluid and the boundary should be controllable and precise since how the error occurs and propagates is uncertain. In addition, in some simulations, we found a few fluid particles broke through the boundary. This phenomenon is unphysical, and it possibly indicated that the repulsive force was not strong enough. We noticed that local uniform stencil boundary condition (LUST) is presented to solve the two problems. The repulsive force and the symmetry have been improved to some extent in LUST. However, for some exceptional cases, such as waterdrop impact, these problems still exist.
This paper presents an improved dynamic boundary condition (IDBC) in SPH methods with three improvements. Firstly, the boundary interaction is enhanced to decrease the number of boundary particles and prevent fluid particles from breaking through the boundary. Secondly, the shape of the support domain attached to the boundary is modified to keep the force field close to the boundary being symmetrical. Finally, a correction factor is introduced to guarantee the unity. Meanwhile, the original kernel is retained, and the algorithm keeps concise.

Method description
The main features of the SPH method have been amply described in the references [9][10][11]. The fundamental idea is to consider that a function f(r) can be approximated by: (1) where: r is position; W is kernel function and h is smoothing length that defines the size of the compact support. This approximation, in discrete notation, leads to: where: a and b are particles; Wab is the kernel function between particle a and b; mk and ρk are the mass and density of particle k respectively. A cubic spline kernel developed by Monaghan and Latanzio [12] is used in this study, as follows: where: q r h = with r being the distance between particles; αD equals to 2 10 7 h  in two-dimensional (2D) conditions and 3 1 h  in three-dimensional (3D). Due to the particular choice of the cubic spline kernel, with the first derivative versus q being zero, the tensile correction method proposed by Monaghan [13] is used to prevent particle clumping. This kernel and tensile correction have been thoroughly proved to be feasible.
The momentum conservation equation in Lagrangian form is: (4) where: ν is velocity; P is pressure; ρ is density; g is gravitational acceleration and Γ refers to dissipative terms. SPH methods describe the momentum equation with different approaches through different formulations of the dissipative terms. The artificial viscosity scheme, proposed by Monaghan [10], is a conventional method within fluid simulation using SPH due primarily to its simplicity. In SPH notation, Eq. (5) can be written as: (5) where the artificial viscosity term ab  is given by: ; α and β are constants. α is related to shear and bulk viscosity. β handles high Mach number shocks and is roughly equivalent to the Von Neumann-Richtmyer viscosity used in finite-difference methods. Following Monaghan [10], β is considered to be zero.
The continuity equation can be replaced either by the interpolant: (8) Eq. (7) comes from Eq. (2) in SPH interpolation. Eq. (8) is the continuity equation in SPH form, in which the density change is computed by solving the conservation of mass [9]. In this study, Eq. (8) is used to calculate the fluid density in order to avoid the artificial density decrease near fluid interfaces.
Following Monaghan [5], Batchelor [14] and Cole [15], pressure can be calculated by the following expression: is the speed of sound at the reference density. Particles are moved using the XSPH method according to Monaghan [10,16]: (10) where: 05 .
With this method, the particle is moved depending on not only its own velocity but also the average velocity of the neighboring particles. The XSPH variant has proven useful in the simulation of nearly incompressible fluids such as water, where it keeps the particles orderly in the absence of viscosity.
The time step is controlled by the Courant condition, the force terms, and the viscous diffusion term [16]. A variable time step Δt is calculated according to Monaghan and Kos [4]: (11) in which: Δtf is based on the force per unit mass f, actually it is equivalent to the magnitude of particle acceleration aa and Δtcν combines the Courant and the viscous time-step controls.
For numerical simulations, the boundary condition is one of the key issues. In SPH method, boundary conditions need to be treated specially. When a particle approaches the boundary, the number of neighbors decreases because there is nothing on the other side of the boundary. The lack of neighbors could result in unrealistic effects as some variables are highly dependent on the summation of the variables attached to the neighbors. To solve this problem, many boundary conditions are developed [3,4,7,17]. One of the successful boundary conditions is dynamic boundary condition [8]. In this method, all the particles, no matter which types they are, follow the same equations of continuity and state, but the positions of the boundary particles do not change unless external forces are imposed. One attractive advantage of this assumption is that the implication is simplified, and computational time can be considerably saved. This boundary condition, first presented by Dalrymple and Knio [18], has been tested [8] and widely used [19,20].

Improved dynamic boundary condition
Dynamic boundary condition has been fully described by Crespo, Gómez, and Dalrymple [8]. The basic idea is to calculate the pressure of the boundary particles with the equation of state. The equation of state can be obtained from the first term of the Taylor expansion of Eq. (9) : where: c is the speed of sound, which is assumed as a constant and ρa0 is the initial density of particle a.
To simplify the description of the problem, only two particles are considered here, i.e., one boundary particle and one fluid particle. In the absence of viscosity and gravity, the momentum equation for the fluid particle (Particle A) becomes: (13) Combining Eq. (13) with Eq. (12) gives: (14) When there are only two particles, Eq. (7) can be written as: As we know, Assuming all the particles have the same mass m, we have: Thus Eq. (14) becomes: (18) With the artificial viscosity as described in Eq. (6) , tensile correction, gravity and more particles, we obtain: (19) where the tensile correction term ten f is defined according to Monaghan [13]: 3.1. Improvements on particle interaction Although the dynamic boundary condition has been widely used and tested, there are still some problems. In some situations, some particles can break through the boundary in an unphysical way. A possible method to solve this problem is to enhance the forces between fluid particles and boundary particles. In Eq. (12), the pressure is expressed using the first-order term of the Taylor expansion of the equation of state. Since the exponent 7  = in Eq. (9), the equation of state can be expanded into seventh order, as follows: where: In another way, Eq. (21) can be considered as a binomial expansion if we express it in the following form: With the same assumption, Eq. (19) can be expressed as: ( ) Considering the particular case of a cubic kernel, the differences of the repulsive forces among different orders of the equation are shown in Fig. 1.
As shown in Fig. 1, the higher-order expansion does enhance the repulsive forces, but the difference between the sixth order and the seventh order is not significant. One reason is that the coefficient of the seventh order is oneseventh in Eq. (23). Another reason is that the base θ is less than 1 and the exponent increases which make the repulsive force decreased dramatically in the seventh order. Another interesting phenomenon is that with the higher-order, the curve's peak becomes higher and moves closer to the boundary particle. It means that an energetic particle needs more energy to climb over the peak and this particle must get close enough to the boundary. Therefore, the repulsive force of the boundary can be enhanced to prevent fluid particles from breaking through the solid boundary, meanwhile it exerts little impact on all the fluid particles. Only the particles close enough to boundary particles will be repulsed significantly. The total potential energy created by the repulsive force can be represented by integrating the dimensionless force over the dimensionless distance from 0 to 2 as shown in Fig. 1. The dimensionless potential energy in different orders is shown in Fig. 2.
With the full expansion, the potential energy is increased by more than seven times. It is an important index to evaluate how IDBC works. Using the full expanding form makes all the fluid and boundary particles follow the same equation of state. However, considering the computation efficiency, we need to determine how many orders should be kept according to each order's strength. Fig. 2 The distribution of the dimensionless potential energy in each order

Improvement on support domain
Another problem is that the field of repulsive forces is not uniform and it is also not specular. According to DBC, the boundary described as a group of particles cannot move, but the forces are also discretized. The field of repulsive force is circular because the support domain of the particles is a circle. The superposition of the force fields of the boundary particles results in the edge of the whole force field being unsmoothed. When fluid particles reach the boundary, they may not be reflected specularly. For some cases, such as droplet impact, the result will be asymmetrical [21]. For an ideal boundary, it is unphysical and unacceptable. Furthermore, an uneven force field may lead to some errors out of control.
We set a given fluid particle at (x0, y0) and some neighbor boundary particles distribute from (x0-Δx, yB) to (x0 +Δx, yB). The density ρB of each boundary particle is the same, and the pressure PB is calculated using the equation of state. We can express the repulsive force on the fluid particle as Eq. When Eq. (25) is expressed in SPH form, the summation of the repulsive forces could not be independent of x if the distribution of the particles is not symmetrical. Thus, we introduced IDBC to make it only related to the distance in the normal direction.
Therefore, we change the form of the support domain of the boundary particles in IDBC. In theory, a boundary particle represents a group of corresponding boundary elements. The contour line of the repulsive force should be parallel to the boundary so that the shape of the support domain is changed from a circle to a rectangle, as shown in Fig. 3. In the compilation, when calculating the distances between fluid particles and boundary particles, the projection in the normal direction of the boundary is only considered. Then the fluid particles at the same distance to the boundary are repulsed equally. The length in the normal direction is 4h, and the length in the tangential direction is two times the distance between the particles, as shown in Fig. 3.
Fluid particles will search neighbor boundary particles in the rectangle support domain. where: n is the normal vector at the boundary. They are used in the governing equations between fluid particles and boundary particles. Eq. (26) is applied to determine the density variation and compute the repulsive force for fluid particles.
As is known, the coefficient of the kernel function is related to the support domain. This coefficient αD is used to keep the integration of the kernel function W(r, h) within the support domain S equal to 1 as expressed in Eq. We select six particles in a row as the boundary and calculate the repulsive force generated by them in a program which is specially developed with C++. The dimensionless forces field in DBC and IDBC with different orders are shown in Fig. 4. The expansion with higher order enhances the repulsive force obviously. This enhancement will prevent fluid particles from breaking through the boundary. The equipotential plane is smoothed by the new shape of the support domain, especially near the boundary. Furthermore, the peak of the equipotential plane is transformed into a ridge. This improvement can make the reflection specular. The form of the kernel function is kept in the vertical direction, ensuring zeroth or first-order consistency at the boundary's edge.

Model performance test
There are some open-source codes for SPH, such as SPHysics [22] and DualSPHysics [23]. However, they are too sophisticated for our simple test cases. For convenience, we developed a program based on C++ and OpenMP to test DBC and IDBC. LUST was tested using DualSPHysics.

Test case 1: Single-particle moving inside a box
The movement of a single particle inside a box is the case used to test the performances of different boundary conditions. With this case, we can make it clear that how the boundary particles repulse the single fluid particles. DBC has been proved to be useful and acceptable in an individual case, as shown in Fig. 5, in which the fluid particle hits the boundary vertically. the boundary particles, and the hollow one is the fluid particle As we know, not all the particles get close to the boundary vertically or symmetrically. In DBC, the field of repulsive force is not uniform. If a particle rushes to the boundary obliquely, the field passed through by the particle is not conservative in DBC, as shown in Fig. 6 (up). It can be seen that the field passed through by the particle is asymmetric in DBC, so is the repulsive force. As a result, in DBC the fluid particle is not reflected specularly, and the energy is not conservative.
With the modified shape of the support domain and the correction factor, this problem is solved in IDBC. As shown in Fig. 6 (down), the modified field in IDBC is symmetric. A test case is designed to compare the differences between DBC and IDBC, as shown in Fig. 7. In a square box with the length of every side being 60 m , a particle moves from (x0, y0) = (30,36) m with an initial velocity ν = (20,-20) m/s. Viscosity and gravity are neglected. The CFL number is 0.05. To make the differences distinct in a short time, the speed of sound c close to the real speed in water is set to 1460 m/s. Although this high value of sound speed is not usually used in simulations, it is high enough to make the differences distinguishable without other unphysical assumptions as the error accumulates over time. Fig. 6 The particle rushes to the boundary obliquely. The big circles and big rectangles represent the boundary particles' support domain in DBC and IDBC, respectively. Solid dots are boundary particles, and the hollow one is the fluid particle Fig. 7 Sketch of the initial setting in the test case for checking the performance on the symmetry of different boundary conditions. The dashed line is the expected route of the fluid particle In this case, the fluid particle rushes to the boundary at 45 degrees in the square box. If the fluid particle is reflected specularly, this particle's trajectory should be a closed rectangle, as shown in Fig. 7. If the boundary generates any error, it will be amplified over time and illustrated by the fluid particle's trajectory. The error will bring unphysical fluctuation to the variables, such as kinetic energy. The results are shown in Fig. 8.
From Fig. 8, it is seen that IDBC is more accurate than DBC for this case. In this extreme test, the fluid particle cannot return the initial position after only one loop in DBC, and the energy fluctuates unphysically. However, IDBC reflects the fluid particle symmetrically and decreases the error generated by the boundary. In the LUST test, the fluid particle was also reflected symmetrically but earlier than that in DBC and IDBC. This phenomenon indicates that the repulsive force generated by LUST is stronger than DBC and IDBC when the particle is at a distance larger than 1.0h to the boundary.
In general, the smoothing radius is chosen to include enough neighbors to limit the impacts of the misalignments between the particles. Therefore, the single particle test case is useful to test the performance of the improved boundary treatment method in detail. a b Fig. 8 The simulated results of the movement of a fluid particle in a box with DBC and IDBC. a) the trajectory of the fluid particle; b) the kinetic energy of the fluid particle

Test case 2: Waterdrop Impact
Waterdrop is an ideal example to test IDBC with a group of fluid particles. This case describes a small waterdrop smashes to the boundary vertically driven by gravity. For this simple test case, if there is no other error, the result should be symmetrical. The waterdrop with a radius of 1.0 mm drops from the height of 2.0 mm with an initial velocity ν = (0,-2) m/s, the artificial viscosity coefficient α = 0.2 for fluid particles and α = 0. between fluid particles and boundary particles. The CFL number is 0.05. The speed of sound is 50 m/s, which is selected to guarantee that the system is weakly-compressible and the computational efficiency is acceptable simultaneously.
To distinguish the differences between DBC, LUST, and IDBC, we set the boundary particles in two ways, as shown in Fig. 9. One is symmetrical and another one is asymmetrical. The separation between the particles, Δx, is 0.2 mm. The displacement between the particles of the two kinds of boundaries is 0.3 Δx.
where: N is the number of fluid particles. This coefficient describes how symmetrical these results are and how the error generated by the boundary conditions affects the accuracy. The results are shown in Fig. 11.
As shown in Fig. 11, all three methods work well with the symmetrical deployment of boundary particles. However, when the boundary particles are deployed asymmetrical (in terms of the waterdrop), the differences between DBC, LUST, and IDBC become obvious. In theory, no matter the boundary particles are symmetrical or not, the results should not be changed too much with the same boundary condition. In IDBC, the result with the symmetrical boundary is close to that with asymmetrical boundary. However, in DBC and LUST, obvious differences appear. These results indicate that with the asymmetrical boundary, the error generated by DBC or LUST is larger than that by IDBC. To test the strength of these three boundary conditions, we also tried to increase the velocity of the waterdrop. With LUST, some particles begin to break through the boundary at 16 m/s, and all the particles escape at 50 m/s. At 65 m/s particle begin to break through in DBC while IDBC still works well that all the fluid particles are reflected into the fluid domain symmetrically. We should notice that in LUST three layers of boundary particles are deployed, and LUST also introduced several fictitious particles for each fluid particle. Although DBC and LUST work perfectly for many other cases, the performance of IDBC is better in symmetry and repulse force.

Conclusions
Improved dynamic boundary condition (IDBC) is developed to solve two problems in DBC. One problem is that the particle may break through the boundary in some situations, such as thin boundary and high-velocity cases.
Another problem is the error generated by the unphysical rough boundary. Three improvements were made in IDBC to solve these problems: a) expansion with a higher order, b) a new shape of the support domain, and c) a correction factor introduced for the kernel function. With these improvements, the potential energy produced by the repulsive force is increased over seven times, and the force field is smoothed.
The validity of the improvements has been firstly tested with a simple case that describes a single particle's movement in a box. In this case, the particle rushes to the boundary obliquely. The performance of DBC is disappointing that the particle departures the expected trajectory in a few steps. This error is decreased remarkably in IDBC.
The validity of the improvements has been further tested with a waterdrop case. This case is tested with different deployments of boundary particles and different boundary conditions. Either the deployment of the boundary particles is symmetrical or asymmetrical, IDBC keeps the results almost the same and stronger than LUST or DBC.