A Dual-Support Smoothed Particle Hydrodynamics for Weakly Compressible Fluid Inspired By the Dual-Horizon Peridynamics

A dual-support smoothed particle hydrodynamics (DS-SPH) that allows variable smoothing lengths while satisfying the conservations of linear momentum, angular momentum and energy is developed. The present DS-SPH is inspired by the dual-support, a concept introduced from dual-horizon peridynamics from the authors and applied here to SPH so that the unbalanced interactions between the particles with different smoothing lengths can be correctly considered and computed. Conventionally, the SPH formulation employs either the influence domain or the support domain. The concept of dual-support identifies that the influence domain and the support domain involves the duality and should be simultaneously in the SPH formulation when variable smoothing lengths are used. The DS-SPH formulation can be implemented into conventional SPH codes with minimal changes and also without compromising the computational efficiency. A number of numerical examples involving weakly compressible. fluid are presented to demonstrate the capability of the method.

Peridynamics [Silling (2000); Silling, Epton, Weckner et al. (2007); Nikravesh and Gerstle (2018); Zhao, Tang and Xue (2018); Shafiei (2018)] is a formulation of solid mechanics that is oriented toward deformations with discontinuities, especially fractures. Though derived from different physical background, Lagrangian smoothed particle hydrodynamics shows great similarity with the state-based peridynamics [Ganzenmüller, Hiermaier and May (2015)]. Peridynamics is regarded as a physics-based method for solid mechanics, and the attempt to solve fluid problem has not been explored. On the other hand, the conventional SPH is restricted by constant smoothing length for all particles, though several numerical techniques have been proposed [Nelson and Papaloizou (1994); Monaghan (2002); Springel and Hernquist (2002)]. Variable smoothing length in SPH can significantly improve the numerical efficiency and allows for adaptive analysis. Dual-horizon peridynamics allows variable horizon domain for each point by defining two horizon domains for each point, where the horizon describes the direct force interaction and dual-horizon represents the reaction force interaction. Based on the similarity between SPH and peridynamics, it is possible to apply the dual support (inspired by the concept of "dual-horizon" in peridynamics) in the field of SPH to solve the variable smoothing length problem.
SPH is a Lagrangian method based on the kernel approximation, where the partial differential equations are transformed into integral form with the so-called kernel interpolation technique [Liu and Liu (2003)]. The computational domain is represented by a set of particles which carry the physical properties, i.e., mass, density, velocity, position, pressure and internal energy. The particles move and properties change with the time due to the interactions between neighboring particles.
One constraint on the kernel approximation is that the smoothing lengths are required to be constant for all particles. However, a computationally efficient SPH implementation requires locally refined regions where variable smoothing lengths can be used. In order to introduce variable smoothing lengths, several methods have been proposed e.g., the averaged kernel method [Hernquist and Katz (1989)] or correction methods with ∇h [Nelson and Papaloizou (1994); Springel and Hernquist (2002); Vacondio, Rogers and Stansby (2012) ;Vacondio, Rogers, Stansby et al. (2013)]. The averaged kernel method uses the averaged smoothing lengths or averaged kernel functions, where the conservation laws can not be satisfied. In the correction methods, the gradient of the smoothing length must be calculated to determine the optimal smoothing length, and a modified coefficient is introduced in all SPH formulations to preserve the conservations laws (more details will be presented in §2). The optimal smoothing length is calculated iteratively, which makes the implementation of SPH more complicated.
The difficulty of varying support domain for SPH arises from that the single support of one particle fails to account for all force interactions with other particles. In contrast with the correction methods aforementioned, we introduce the concept of "dual-horizon" in dual-horizon peridynamics [Ren, Zhuang, Cai et al. (2016); Ren, Zhuang and Rabczuk (2017)] into SPH for varying the smoothing length for each particle. Meanwhile, the dual-support SPH in fluid may be viewed as the fluid version of dual-horizon peridynamics based on the similarity of SPH and peridynamics in solid.
In this paper, we develop a SPH formulation which naturally includes variable smoothing lengths. The formulation named dual-support smoothed particle hydrodynamics (DS-SPH) is simple, satisfies the conservation of linear momentum, angular momentum and energy when variable smoothing lengths are employed. The paper is outlined as following. In §2, the theoretical background of SPH is reviewed. In §3, the concepts of support and dual-support are introduced, based on which we present a general method to reformulate the SPH equations. In §4, we convert the traditional SPH formulation into the dual-support SPH. In §5, the conservations of momentum, angular momentum and energy are proved. In §6, the implementation of DS-SPH is discussed. Three numerical examples are tested to demonstrate the performance of the DS-SPH method in §7.

Kernel function
Let i denote the point in the global domain Ω; r i being the current coordinates vector of point i, H i is the support domain associated with i with a radius of h i . The kernel function W i is defined as where r ij = |r ij |, r ij = r i − r j . It can be seen that if There are a wide variety of kernel functions including the Gaussian kernel, cubic-spline, quintic spline or Wendland kernel [Wendland (1995)], to name only a few. For more discussions on these kernel functions, readers are referred to Liu et al. [Liu and Liu (2003); Dehnen and Aly (2012)]. As the support in SPH fluid is often defined by the current configuration of the particles, the kernel function is termed as the Eulerian kernel. Aside from the Eulerian kernel, there exists another type of kernel function called Lagrangian kernel whose support domain is based on the material coordinates in the initial configuration. The Lagrangian kernel is capable of overcoming the tensile instability when large deformation is encountered in the context of SPH for modeling solid [Rabczuk, Belytschko and Xiao (2004)]. However, the Eulerian kernel is more suitable for fluid simulation compared with the Lagrangian kernel [Rabczuk, Belytschko and Xiao (2004)].

Kernel interpolation
Kernel interpolation theory starts with the identity where f is an arbitrary scalar variable, δ is the Dirac delta function and V j is the volume associated with point j. This integral is approximated by replacing the delta function with a smoothing or kernel function W with finite width h i , e.g., where W has the property Integration by parts and neglecting all surface terms, the derivative of f (r i ) is derived as Hence the derivative of a function is transferred to the kernel function by kernel interpolation theory.

Completeness of kernel function
The completeness is related to the capability of the kernel approximation to reproduce a polynomial function of certain degree exactly. Zero-order completeness is necessary to represent the rigid body modes while first-order completeness is a necessary condition for constant strain state. According to Belytschko et al. [Belytschko, Krongauz, Dolbow et al. (1998)], first-order completeness is defined by First-order completeness for the derivatives of the kernel approximation for any particle i in 1, 2 and 3 dimensions should satisfy Here, I is the identity matrix. SPH in its continuous form (or integral form) fulfills zero-and first-order completeness and hence the 'standard' kernel function as given in Eq. (5) is sufficient. However, not even zero-order completeness is guaranteed for the discrete SPH form as shown by Nguyen et al. [Nguyen, Rabczuk, Bordas et al. (2008)]. Therefore, many correction methods have been developed such as the symmetrization proposed by Monaghan [Monaghan (1988)], Randles & Libersky correction [Randles and Libersky (1996)], Johnson & Beissel correction [Johnson and Beissel (1996)] and Krongauz-Belytschko corrected derivatives [Krongauz and Belytschko (1997)]. In the current paper, two of those approaches are briefly summarized.
The gradient correction of kernel function [Bonet and Lok (1999)] in discrete form satisfying Eq. (8) is where ∆V j = m j /ρ j . The gradient correction guarantees the conservation of angular momentum. The zero-order (Shepard filter) is given bỹ Such correction is often adopted to reduce the pressure oscillations for particles near the boundaries or close to free-surfaces when density summation is used instead of the continuity equation. Based on Eq. (10), the kernel gradient correction (or corrected derivatives method) is given bỹ It is not difficult to verify that both Eq. (9) and Eq. (11) satisfy Eq. (8).
The gradient correction Eq. (9) is much simpler than the mixed correction Eq. (11). In the current paper, the gradient correction is employed. The formulations with zero-order correction or gradient kernel correction are straightforward in SPH. It requires only the replacement of the kernel functions with their associated corrected kernels, i.e., W ij →W ij and ∇ i W ij →∇ i W ij or ∇ i W ij →∇ iWij . Subsequently, we derive all equations based on the standard SPH formulation but the equations for the associated correction methods are straight forward.

Variable smoothing length in SPH
There are several methods to deal with the variable smoothing length issue. The first method is the averaged kernel function [Hernquist and Katz (1989)] given by The averaged kernel between paired particles with different smoothing lengths guarantees the anti-symmetrical pair-wised forces and thus preserves the symmetry of the particle interactions [Liu and Liu (2003)]. The second method is based on the correction term ∇h [Nelson and Papaloizou (1994); Springel and Hernquist (2002); Vacondio, Rogers, Stansby et al. (2013)] in the equations of motion. The key idea is to relate the local number density of particles with the smoothing length and to keep the variable inside the smoothing sphere constant [Springel and Hernquist (2002)], i.e., The equation of motion in Springel et al. [Springel and Hernquist (2002)] is where N is the number of all particles, p i is the pressure and f i is defined as The quantities ∂ρ i ∂h i can be computed along with the densities themselves; more details are found in Springel et al. [Springel and Hernquist (2002)].
Monaghan [Monaghan (2002)] proposed a similar formulation for the smoothing length. In his work, the momentum equation is where Ω i is given by The coefficient Ω i is also introduced in the continuity and energy equations. Both expressions in Eq. (15) and Eq. (17) are related to the term ∇h i . The basic idea is to find the proper value h i which satisfies Eq. (13). In order to find the desirable smoothing length, Eq. (13) is transformed into a set of two simultaneous equations which is computed at the location of particle i, where η is a parameter specifying the smoothing length in units of the mean particle spacing (m/ρ) 1/d , d is the number of dimensions. These two equations can be solved simultaneously using standard root-finding methods such as Newton-Raphson or Bisection [Price (2012)]. However, this inevitably increases the computational cost and make the SPH implementation more complex.

Support domain and dual-support domain
In this section, the key concepts of support and dual-support are presented and the duality of dual-support is proved. The new concept provides great flexibility to convert the traditional constant support SPH to dual-support SPH allowing for variable smoothing length for each particle.

Support and dual-support
The conventional SPH formulations based on variable smoothing lengths consider the unbalanced interaction with averaged kernel functions, correction terms ∇h. The basis of all these methods is a single support domain. However, the single support domain cannot elegantly resolve the unbalanced interaction with different support radii. In SPH with variable smoothing lengths, one common situation as shown in Fig. 1 The unilateral interaction violates the Newton's third law. Therefore, a single variable support is not sufficient to define the interactions between particles and the new concept of support and dual-support is introduced subsequently. A similar concept is the "dual-horizon" in peridynamics [Ren, Zhuang, Cai et al. (2016)].

Support
The support H i for point i is defined as a domain related to i. When the domain is centered at i with a radius of h i , the support H i can be given as In traditional SPH, the support of a particle is a domain where the kernel function is positive (see p60 in Liu et al. [Liu and Liu (2003)]). The definition in Eq. (19) implies the positivity of kernel function since the kernel function is calculated in the support, and the kernel function is zero outside of the support domain. Therefore, the support in Eq. (19) is the same as that in conventional SPH. Dual-support of i is defined as a union of the points whose support domain includes i, denoted by In the notation of dual-support H i , the superscript prime indicates "dual", and the subscript i denotes the object particle. When the support is defined as circle or sphere centered at that point, i's dual-support can be expressed as For any point i, the shape of H i is arbitrary and depends on the sizes and shapes of supports as well as the locations of the particles. For example, as shown in Fig. 2, the dual-support with respect to point 0 contains particles {1, 2, 3, 4}, whose supports are denoted by thin solid circles. Particles {5, 6} are not included in the dual-support of point 0 since their supports do not include point 0. For case with constant smoothing length, as j ∈ H i ⇔ i ∈ H j , H i is equal to H i and then support and dual-support degenerate to the constant support in traditional SPH.
It is worth mentioning that the dual-support domain with respect to the support domain is similar to the relation between the influence domain and support domain in meshfree methods [Hernquist and Katz (1989); Liu and Liu (2003)]. The dual-support domain actually represents the influence domain. However, to the authors' knowledge, the duality of the influence domain and support domain was not identified in the literature. In the conventional SPH formulation, either the support domain or the influence domain is employed. In our study, the influence domain and the support domain should be employed simultaneously in the SPH formulation when variable smoothing lengths are used. For the simplicity of derivation, we use the dual-support domain to represent influence domain.

The duality of dual-support and support domain
Let F(i, j) be any expression depend on two points i, j. The duality of dual-support-which we have proven in Ren et al. ]-states that the double integral of the term F(i, j) in dual-support domain can be converted to the double integral of the term F(j, i) in support domain: The key idea lies in the fact that the term F(i, j) can be both interpreted in H i and H j . For the sake of completeness, the proof of duality of dual-support is provided in Appendix A.

The dual formulation based on support and dual-support
Let ρ be any scalar, f be any scalar field. The SPH formulation of ∇f i ρ i can be obtained as In order to derive the DS-SPH formulation, we simply exploit the duality, Eq. (23), and replace the integral over the support of the first term in Eq. (25) with the integral over the dual-support: This obviously yields the DS-SPH formulation: Note that for points with identical smoothing length, (25) forms the key step to convert traditional constant support SPH to dual-support SPH. Let φ be any scalar and A be any vector. With the operation defined in Eq. (25), the following dual support formulations can be devised The second order derivative can be obtained in a similar way. One expression recommended by Brookshaw et al. [Brookshaw (1985); Morris, Fox and Zhu (1997)] for second derivative is which requires only first spatial derivatives, where v ij = v i − v j , η is a small number introduced to avoid a zero denominator during computations and is set to 0.1h. Eq. (29) is the combination of SPH formulation of first derivative and first order finite difference.
4 Dual-support smoothed particle hydrodynamics In this paper, we apply the dual-support formulation to fluids. The application of DS-SPH in solid is refereed to Dai et al. [Dai, Ren, Zhuang et al. (2017)]. The Lagrangian form of the fluid dynamics equations is given as where d dt is the material time derivative, σ the Cauchy stress; q = −κ∇T the heat flux density vector, κ the thermal conductivity, T the temperature,ṡ the source term; b the body force density, the internal energy density, e = v 2 2 + the specific energy per unit mass. Eq. (33) and Eq. (34) are equivalent. Two constitutive models for the fluid are where p is the pressure; µ is the dynamic viscosity, and λ is the second coefficient of viscosity. As a detailed derivation of the SPH equations is provided elsewhere, see e.g., [Liu and Liu (2003)], we omit it in this manuscript. The continuity equation in SPH form is given by [Liu and Liu (2003)] In the Lagrangian formulation, the conservation of mass is naturally satisfied, hence we are not going to recast the continuity equation into the dual-support formulation. The continuity equation is often replaced by the summation density approach

SPH formulation with dual-support
It can be shown [Randles and Libersky (1996)] that the SPH formulation with constant smoothing length for the equation of motion is given by and the energy equation by where the term on heat conduction is similarly expressed by Eq. (30).
Eq. (40) can be simplified via the momentum equations Eq. (39), which leads to the energy equation based on the internal energy In order to allow variable smoothing lengths, applying the procedure in Eq. (25)-Eq. (39), Eq. (40) and Eq. (41), we obtain the following dual-support formulations for the equation of motion and energy equations,

Two special cases
For fluid with material constitution given by Eqs. (35), (42) and (44) can be further simplified into, respectively where Π ij is the artificial viscosity, which is necessary to model strong shocks and prevent particles from interpenetration [Benz (1990)] in the absence of physical viscosity and is given by ; c i is the sound speed associated to particle i, α and β are constants that are all typically set around 1.0 [Monaghan (1988)].
Another viscosity in fluid is the physical viscosity. When the artificial viscosity in Eqs. (45) and (46) is replaced with the physical viscosity using Eq. (30), SPH can be reformulated as,

Conservation of the basic laws
In this section, the conservation laws, i.e., the linear momentum conservation, angular momentum conservation and energy conservation, in DS-SPH are discussed. Due to the Lagrangian description of motion, the total mass is conserved naturally as pointed out previously.

Conservation of linear momentum
The conservation of momentum requires which represent the continuous form and discrete form, respectively. The variation of momentum is derived as dP dt = In the fourth step, the duality of dual-support is used. Therefore, the conservation of momentum is satisfied.

Conservation of angular momentum
The conservation of angular momentum requires In the proof of angular momentum, we replace { 1

The variation of angular momentum is derived as
In the fourth step, the duality of dual-support is used. Temporally, for the sake of simplicity, we abbreviate H i → H, ∆V j → ∆V , r ij → r, W ij → W and σ i → σ . Let capital letters I, J, K, L be the dimensional index.
where IJL is the permutation symbol and r I refers to the I-th component of r,∇W is the function with gradient correction, and the symmetry of the Cauchy stress tensor σ and the completeness of∇ i W ij are used. Therefore, Eq. (52) is zero, and the angular momentum is conserved when the linear completeness is satisfied.

Conservation of energy
The conservation of energy requires dE dt = The variation of specific energy is derived as In the third step, the duality of dual-support is used. So the specific energy is conserved.

Implementation of DS-SPH
At each step, the values at step t + ∆t are calculated based on the known variables as r Many numerical schemes can be used to integrate the SPH formulation, e.g., the Leap-frog prediction-correction scheme and the Verlet-velocity scheme [Verlet (1967)].
When the energy equation is considered, the Leapfrog prediction-correction scheme is preferred. The Leapfrog prediction-correction scheme comprises three steps as below 1. Prediction step DS-SPH requires the summation over the support and dual-support for each particle. However, the forces from the dual-support domain are obtained from other particles' support domains. For any particle j ∈ H i , the force f ij for pair ij is calculated and added to particle i. Since j ∈ H i ⇔ i ∈ H j , −f ij is the force from dual-support H j of particle j, thus it can be added to particle j. Hence, when summing over all other particles, the forces from H i are automatically calculated. In other words, the force from one particle's support domain is reusable for the other particle's dual-support domain.
Let us consider Eq. (45). In the implementation, we loop over all particles j ∈ H i and add the force 2 )∇ i W ij to particle i. Since j ∈ H i ⇔ i ∈ H j , we can simply add −f ij to particle j. In contrast, traditional SPH requires computing −m j ( p i Hence, DS-SPH is more efficient compared with the traditional SPH. In contrast with the kernel average method which depends on two particles' kernel functions, any particle in DS-SPH only concentrates on its own kernel function and support domain, which is independent with the other particles. Eq. (18) indicates that the smoothing length is inverse proportional to the density. Let ∆x i denote the 'size' of particle i. Without loss of generality , we assume that the shape of the particle is a cube in 3D or square in 2D. Then ∆x i = m i /ρ i , where d is the number of dimensions. Since each particle's mass is fixed, ∆ x i decreases when the density ρ i i ncreases. At m(≈ 50) time steps, we update the support radius based on the new estimated particle size by where n ≈ 1.5 ∼ 3. It can be seen that the smoothing length decreases with increasing particle density; the low(high) density, which indicates relatively sparse (dense) neighbors, causes the support domain to expand (shrink) adaptively.

1D heat conduction
Consider a 1D bar of with length L =50 cm. The bar is discretized with a particle spacing of ∆x =1 cm or ∆x =0.5 cm; the heat diffusion coefficient is α = 1.0 × 10 −4 m 2 s −1 . The left half is assigned an internal energy of e 0 l = 1 J/m, the right half e 0 r = 2 J/m. The total energy in the 1D bar is E total =0.75 J. There is neither potential nor kinetic energy present in this simulation. The only non-zero contribution is the total internal energy. The energy profile given by is compared to the analytic solution at 4.0 s. Three cases as shown in Tab. 1 are considered. The first case is modeled by conventional SPH with constant smoothing length. The second case is simulated by conventional SPH with variable smoothing lengths but without additional treatment. The third case is implemented with our dual-support SPH. In order to test the influence of the transition of smoothing length, two particle spacings ∆x =1 cm, ∆x =0.5 cm in Case II and Case III were employed , more specific, the particle spacing in the interval of x = 0.3L − 0.5L is ∆x =0.5 cm. The L 2 error in internal energy is given by The numerical results agree well with theoretical solution, as shown in Tab. 2. The error of the conventional SPH formulation with varying smoothing length is roughly twice as high as the error with the dual-support SPH version.

Water droplet flow
In this section the flow simulation of an elliptical water droplet is presented. The main purpose of this example is to demonstrate that the total linear and angular momentum are conserved for the dual-support SPH in the absence of external forces. This example is tested by two formulations, one is the Eq. (39) with only one varying support for each particle (namely, single support SPH), the other is Eq. (42) (dual-support SPH). No treatment of varying support radii is employed for the first formulation. The gradient correction of kernel function (cubic-spline kernel) is employed in two formulations. The continuity equation Eq. (31) is calculated directly by the divergence of velocity. The velocity gradient is calculated by The geometry of the bubble is a circle of 1 m radius without external forces but with initial velocity field of (−100x, 100y) m/s. The bubble is discretized with 1979 particles, whose distribution is shown in Fig. 3. The cell's volume represents the particle's volume. In order to test the influence of variable smoothing lengths on the conservation laws, the upright part of the bubble is discretized with small particles. The initial smoothing length is selected as two times of the particle size, i.e., h i = 2∆x i . The particle size is estimated by assuming the shape of square.
The total stress is calculated by Eq. (36). The pressure is calculated by the following equation of state for water [Monaghan (1994)], where γ is a constant; γ = 7 in our simulations; ρ 0 is the reference density, c 0 = 1400m/s is the sound speed at the reference density. p 0 = c 2 0 ρ 0 /γ is the artificial bulk modulus [Morris, Fox and Zhu (1997)]. In this example, the artificial bulk modulus is p 0 = 285.714 MPa, viscosity coefficients λ = 0 and µ = 0.5 kg m −1 s −1 , initial density ρ 0 = 10 3 kg/m 3 . Within the simulations, the shape of the bubble should remain elliptical, the value of ab and w is the initial value of ab. Fig. 4 shows the geometry of the water droplet simulated by single support SPH and the dual-support SPH. The result by DS-SPH agree well with the theoretical shape denoted by solid lines, while the result by traditional SPH using only one support domain is affected by the particle distribution. The reason is that the traditional SPH formulation is only suitable for constant smoothing length; when encountering the variable smoothing length, the accumulated unbalanced forces drive the final result away from the theoretical solution.
In addition, when tracking the total linear momentum (Fig. 5), and total angular momentum ( Fig. 6), the results by single support SPH change significantly with time, whereas that by dual-support SPH fulfill all conservation laws. The non-zero initial total momentum and angular momentum are due to the unsymmetrical discretization. The variation of linear momentum and angular momentum for single support SPH is due to the variation of the smoothing length for each particle at every 50 steps. Figure 3: The particle discretization of water droplet. Each cell represents one particle

2D dam break over dry bed
The dam breaking experiment, which was described in Koshizuka et al. [Koshizuka and Oka (1996)], is a benchmark problem to test the accuracy of SPH code by Violeau et al. [Violeau and Issa (2007)] and Crespo et al. [Crespo, Gómez-Gesteira and Dalrymple (2007)]. The tank is 4 m long, the initial volume of water is 1 m long and its height 2 m, as shown in Fig. 7. The system is solved with a leapfrog prediction-correction scheme, using a cubic-spline kernel without kernel gradient correction, specular reflection boundary condition by Eq. (61), artificial viscosity, α = 1, β = 1. The density is calculated by Eq. (37). Fluid particles were initially placed on a staggered grid with zero initial velocity. In order to employ a large time increment, the sound speed is set as 100 m/s, which is 10 times larger than the maximum flowing speed. The specular reflection boundary is given by Total linear momentum x, single support SPH y, single support SPH x, dual−support SPH y, dual−support SPH where v denotes the velocity, n is the inward normal direction of the wall at that point. Note that the specular reflection boundary only changes the direction of the particle's velocity when the particle is approaching the boundary, thus the kinetic energy is not affected. Such good property enables us to track the global energy during the simulation. Although the equation of state given by Eq. (59) is irrelevant to the energy state, the energy equation Eq. (46) is considered in order to track the internal energy since the artificial viscosity converts the kinetic energy into the internal energy.
Three cases are run to shown the capabilities of the dual-support SPH. Case I is based on traditional single support SPH with variable smoothing length and Case II adopts the DS-SPH formulation. Case III uses constant smoothing length during the simulation. The initial smoothing length is set as h i = 2∆x i . The smoothing length for each particle in Case I and Case II is updated with Eq. (55) at every 200 time steps. The particle spacing of Case I (5,600 particles) and Case II (5,600 particles) is ∆x = 2.5 × 10 −2 m in Ω 0 , ∆x = 1.25 × 10 −2 m in Ω 1 , whereas only one particle spacing ∆x = 1.25 × 10 −2 in Case III (12,800 particles). There is a sharp change of smoothing length in the interface of Ω 0 and Ω 1 for Case I and Case II. Such transition deteriorates the traditional single support SPH, while dual-support SPH can reduce the adverse effect. The Case III based on traditional constant support SPH is served as the reference.
As shown in Fig. 8, the toe velocity of Case II and Case III agreed well with experimental data, whereas that of Case I was affected by sharp transition of smoothing lengths.
The fluid domain is marked with four colors so that the deformation of the interfaces can be tracked. The deformation of water column for three cases at different time is shown in Fig.  9 and Fig. 10. For case I, the blue zone pushed the red zone and caused the simulation to deviate far away from the reference results given by Case III; the interfaces became irregular and the blue zone spreaded over the bottom. The results of Case II agreed well with that by Case III; the interfaces were continuous and smooth in all steps. Hence the dual-support SPH can reduce the adverse effect of variable smoothing lengths to minimum. When the front toe hit the right wall, for Case II and Case III, the maximal density variation for front toe is smaller than 2%, as shown in Fig. 11 Figure 8: The evolution of dam front toe in experiment [Koshizuka and Oka (1996)  The total energy for three cases were compared in Fig. 12. Fig. 12(a) shows that the total energy of Case I increased over time, which indicates that the spurious force has a great impact on the global energy conservation. In fact, the spurious force is the unbalance force interaction between two particles with different smoothing lengths. For Case II and Case III, the total energy variation arising from the time integration is less than 4% . Therefore, the dual-support SPH formulation conserves the global energy in the absence of external force doing work. The dual-support formulation is a good alternative to varying the smoothing lengths in SPH.

Conclusions
In this paper, we have introduced the dual-support domain, the dual part of support domain, based on which the conventional SPH was reformulated into dual-support SPH.
We identified that the duality of influence domain and support domain in conventional SPH method. This formulation enables the conservations of momentum, angular momentum and energy when variable support domains are employed. In DS-SPH, the change of smoothing length is achieved with ease and the modified coefficients in Eq. (14) and Eq. (16) are eliminated. The implementation of DS-SPH shows that the DS-SPH can reduce the computational cost compared with traditional SPH.
Three numerical examples were presented to validate the dual-support SPH. The first numerical example showed for heat diffusion problem the energy is conserved when variable smoothing lengths are utilized. The third numerical example verified that the conservations of basic laws for DS-SPH are well preserved while the traditional SPH not.
The last example showed traditional SPH is greatly affected by the variable smoothing lengths and DS-SPH can eliminate the adverse effect.
The concept of dual-support facilitates the support-variable SPH formulation. To some extent, the concept of support domain and dual-support domain is similar to Newton's third law, considering the direct force and reaction force for paired particles. The proof of conservation is obtained easily with the aid of duality of dual-support and support domain. Therefore, the dual-support can be also applied in the SPH on solid or SPH on magnetohydrodynamics. The present method is also promising for multi-scale analysis where the models with different length scales can be bridged by using different smooth length settings.