Simulation of Droplet Impacting on Elastic Solid with the SPH Method

The phenomenon of droplet impacting on solid surfaces widely exists in both nature and engineering systems. However, one concern is that the microdeformation of solid surface is difficult to be observed and measured during the process of impacting. Since the microdeformation can directly affect the stability of the whole system, especially for the high-rate rotating components, it is necessary to study this phenomenon. Aiming at this problem, a new numerical simulation algorithm based on the Smoothed Particle Hydrodynamics (SPH) method is brought forward to solve fluid-solid coupling and complex free surface problems in the paper. In order to test and analyze the feasibility and effectiveness of the improved SPHmethod, the process of a droplet impacting on an elastic plate was simulated. The numerical results show that the improved SPH method is able to present more detailed information about the microdeformation of solid surface. The influence of the elastic modulus of solid on the impacting process was also discussed.


Introduction
The phenomenon of droplet impacting on solid surfaces widely exists in nature.Simulation of this kind of problems has always been a difficult and important research area in the computational fluid dynamics (CFD).Two basic numerical methods used nowadays are developed from Euler's theory and Lagrange's theory, respectively.Over the past few decades, the grid-based methods based on Euler's theory, such as the finite element method (FEM), the finite volume method (FVM), and the finite difference method (FDM), have been widely applied in the simulation of the flow.However, dealing with large deformations of the moving free surfaces and interaction of the multiphase flows are difficult to achieve for the grid-based methods.Some methods which are good at capturing the free surface and regenerating the grid have been proposed, such as PIC [1], MAC [2], VOF [3], and LS [4] methods.Moreover, mesoscopic simulation methods have been proposed including MPS [5], LBM [6], and DPD [7].Although these methods can deal with the unsteady flows with complex free surfaces, they also lead to high computational complexity and large amount of data processing.
Smoothed Particle Hydrodynamics (SPH) method which is a meshless method based on Lagrange's theory was proposed by Lucy [8] and Gingold and Monaghan [9] in the 1970s.The basic idea of the SPH method is to discretize the computed domain into a large number of particles instead of the fixed grids.These particles can freely move according to the governing conservation equations but be restricted by the material properties.The SPH quantities can be obtained as weighted averages from the particle values over the region of interest.Therefore, the SPH method is able to handle the complex free surfaces and the large deformation problems easily.Recently, the SPH method has been adopted for various problems, including Newtonian fluid flows [10], non-Newtonian fluids flows [11][12][13][14], incompressible fluids [15][16][17], multiphase flow [18][19][20], free surface flows [21,22], large deformation [23], and the dynamic response of elastic-plastic materials [24][25][26].
Droplet impacting on a rigid plate is a typical instance of free surface flows and attracts much attention.In recent years, different methods have been used to simulate this problem and some useful conclusions were obtained.For example, Tomé et al. [26] used the finite difference method to simulate

Introduction of the SPH Method
2.1.Basic Principles of the SPH Method.In the SPH method, the computed domain is discretized into several continuous particles with material properties.These physical quantities are obtained by integral representation of function [34] as follows: where ⟨()⟩ is the approximation of (),  is the vector position, the smoothing length ℎ defines the influence area of the kernel function , and  satisfies the following conditions: ( The integrated form of () can be discretized as a summation over all the particles in the influence domain as follows: When the kernel function is differentiable, the approximations derivative of () based on (3) is derived as 2.2.Kernel Function.The kernel function is a key element in the SPH method to ensure the accuracy of the algorithm.Many possible forms of the kernel have been analyzed and compared in these literatures [34][35][36][37].In order to balance the computational accuracy and efficiency, the cubic spline kernel is adopted as one of the most widely used kernel functions in SPH method [10].Therefore, the cubic spline function is chosen as follows: 0,  ≥ 2; (5) in two dimensions, the normalization factor is   = 15/7ℎ 2 .

Governing Equations.
Under the tremendous impact, the solid material deforms obviously and solid particles characterized by SPH method move like fluid.Therefore, the governing equations for high strain hydrodynamics with material strength were proposed in these literatures [24,38] to simulate the impacting and penetrating problems.Those equations including the mass and momentum conservation equations have exactly the same form as those for transient compressible fluid flow.In a Lagrangian frame, the governing equations are written as [34] where  denotes the density,  is the mass of particle,  is time, and   and   are the velocity and the acceleration due to external forces of the th component, respectively.  is the (, )th component of the stress tensor.Calculation of   for solid and fluid is different and more details are shown in Section 2.5.2.

Artificial Viscosity and Artificial Stress
2.4.1.Artificial Viscosity.In SPH method, the artificial viscosity is employed to allow the algorithm to be capable of modeling shock waves or simply to stabilize a numerical scheme.Considering the artificial viscosity Π  in (7), the momentum equation could be obtained as The most common form of the artificial viscosity suggested by Monaghan [10,39] is where   = (  +   )/2,   = (  +   )/2,   = ℎ]    /( 2  +  2 ),   =   −   , and c denotes the sound speed.The parameter  2 for   prevents singularities and it should be small enough to prevent severe smoothing of the viscous term in the high density regions.Normally, this is achieved by taking  2 = 0.01ℎ 2 , which means that smoothing of velocity will only take place if the particle spacing is less than 0.1ℎ [10].Moreover, the  Π -term produces shear and bulk viscosity, while the  Π -term is similar to the Von Neumann-Richtmyer viscosity to handle high Mach number shocks.The values of  Π and  Π are not critical, but they have significant influence on the computed results.Monaghan [10] suggested their values be around 1.0 and 2.0, respectively.But the  Πterm may result in larger shear viscosity when modeling flows with physical viscosity.It has been proposed that this term should be removed by setting  Π = 0 and the  Π -term should be retained to prevent unphysical particle penetration [16].For the rigid boundary case without the  Π -term of artificial viscosity, there was fracturing phenomenon in the Newtonian impacting droplet, but the computation can continue [31,32].However, unlike the rigid-plate case, for considering the deformation of elastic-plate and fluid-solid interaction under impact in this paper, it is found necessary to retain the  Πterm to avoid the divergence and penetration in simulation.The detailed results given in Section 3.2 show that values 1.0 and 2.0 are suitable for  Π and  Π , respectively.

Artificial Stress.
In SPH method, the tension force can introduce instability for the solid deformation, while the tensile stress can also become unstable for the free surface.There are many methods to improve the stability of tension and shear.The most commonly used one is the "artificial stress" proposed by Monaghan [40] and Gray et al. [41].The basic idea of artificial stress is to introduce a small repulsive force between a pair of neighboring particles to prevent them from getting too close when they are in a state of tensile stress.To do this, the momentum equation is modified as The detailed formulation of the artificial stress   is as follows: where  = (0, ℎ)/(Δ, ℎ),   =   (|  −   |, ℎ)/  (Δ, ℎ), and Δ is the initial distance of particles.The components of the artificial stress tensor   are given by [31,32]   =   cos 2  +   sin 2 , where In these equations,  is a parameter and it is set to 0.2 for the Newtonian fluid [31,32].

Fluid-Solid
Coupling.The processing of two-phase coupling is a popular interest of research in computational fluid dynamics.Many techniques have been proposed to accurately and effectively simulate the interface.However, most of them only analyzed the interaction of solid on fluid, such as flow around a cylinder, or the interaction of fluid on solid like feathers in the air movement.Obviously, the interaction between two phases should be analyzed to ensure the accuracy of the simulation results.In this work, an improved algorithm of fluid-solid coupling based on the SPH method is presented.The specific algorithm is as follows and the flow chart for the numerical simulation is shown in Figure 1.

The Solution of Density.
In the SPH method, density is a critical variable in solution of the governing equations.
(1) The Density of Fluid or Solid.In general, density can be directly calculated by (3).However, integral on truncated boundary for the free surface case may result in large errors.A method which has been proposed by Monaghan [39] can be used to solve this problem.The key idea of the method is that the reference density is equal to the real density of the object and then the increase of density is calculated by (6).Finally, density of the computed domain is obtained at the end of each time step.
(2) The Density of the Fluid-Solid Interface.For the fluidsolid coupling problem, the calculation accuracy of density on the interface is very important.If particles of different materials in the domain are not regarded as neighboring particles, the change of density cannot reflect the actual status.On the contrary, particles of different materials as The flow chart for the numerical simulation (the particle  is a certain particle and the particle  is the particles in the influence domain of particle ).
neighboring particles contributing to calculation can cause severe distortion of density.A correction algorithm for the calculation of density is proposed as In this algorithm, particles of different materials in the influence domain are regarded as neighboring particles, and the numerical results prove that the above algorithm is correct and practical.Moreover, (15) can be written as (6) when the particles in the influence domain are the same material with the same density (  =   ).Therefore, ( 15) is adopted for the solution of density..The key of momentum equation ( 10) is to calculate the stress tensor.

The Solution of Momentum Equation
In this paper, the calculations of solid and fluid are totally different because the droplet is treated as a Newtonian fluid while the plate is treated as an elastic solid surface.Besides, the deformation degree of plate and the effect of deformation on the droplet should be considered.Therefore, the solution of momentum equation for solid and fluid is as follows.
(1) The Extra Stress Tensor of Fluid.The extra stress tensor   is calculated based on the isotropic pressure  and the additional elastic stress   .Consider The constitutive equation for the additional elastic stress   of Newtonian fluid is as follows: Mathematical Problems in Engineering 5 The equation of state is as follows: where  is the dynamic viscosity of fluid,  0 is the reference density, and  is the speed of sound at the reference density.For compressible fluid flow, the density variation is proportional to the square of the Mach number  2 ( 2 ≡ /, where  is a typical reference velocity).In this paper, we control the density variation to within 1%.On the fluid-solid interface, solid particles in the influence domain of the fluid particle should be contributed to calculate the additional elastic stress of the fluid particle according to (17).
(2) The Extra Stress Tensor of Solid.The solid particle as elastic solid accords with the relation of stress-strain: where  is the elastic modulus and is also called Young's modulus; ] is Poisson's ratio.On the fluid-solid interface, fluid particles in the influence domain of the solid particle should be also contributed to calculate the strain of the solid particle according to (20).
(3) The Fluid-Solid Interface.The treatment of boundary condition is one of the difficulties in the simulation based on SPH method.In this paper, in order to avoid errors caused by integral on truncated boundary condition, different material particles in the influence domain of the particle  should be contributed to the solution for the particle .In other words, solid particles in the influence domain of the fluid particle should be regarded as fluid particles to solve the governing equations of fluid, while the governing equations of solid should be calculated in the same way.It means that the force exerted by solid particle on fluid particle has to react in the opposite direction on the solid particle.Therefore, the interaction between solid and fluid has been considered in the solution for the governing equations on the interface.However, nonphysical penetration and particles doping may occur on the fluid-solid interface during the process of impacting.The algorithm of interface proposed by Liu et al. [42] is quite effective in simulation of high-velocity impact and underwater explosion.There is a pair of balancing forces on two particles when the two particles are approaching each other.And this force has the same direction as the centerline of particle.Consider where  1 and  2 are equal to 12 and 4, respectively. depends on the specific application, and it is usually set to about 10 times the maximum velocity of fluid field.
Considering the balancing forces, the momentum equation on the fluid-solid interface could be obtained based on (10). Consider 2.6.Time Integration.In this work, Leap-Frog (LF) method is chosen for solving the time-dependent ordinary differential equations of the SPH method.In order to prevent numerical divergence, the time-step Δ should satisfy the following restricted conditions [34,38]: where  is the hydrodynamics force (per unit of mass) acting on the particle .The speed of sound for solid is even more complicated and difficult to calculate than fluid.And the normal formula for calculation is as follows:

Numerical Results
The two-dimensional Newtonian droplet impacting on a rigid plate as a typical instance for free surface flows has been investigated by other researchers [26,31,32].Based on the related works, the plate in this paper is considered to be an elastic solid surface and the model is simulated by present method to analyze the complex free surface of droplet and the deformation of solid.

Description of Model.
In order to compare the results from different methods, the parameters of droplet given in these literatures [26,31,32] are adopted in this section.The initial diameter and velocity of droplet are  = 0.02 m and  = −1 m⋅s −1 .The density and dynamic viscosity of droplet are set to   = 1000 kg⋅m −3 and  = 4 Pa⋅s.At the initial time, the droplet is in a stress-free state and the center of the droplet is at position (0.0 m, 0.04 m).Then under the influence of gravitational forces taken as  = 9.81 m⋅s −2 , the droplet begins to fall along the -direction.The computations and results under different initial particle spacing were obtained and analyzed in literatures [31,32].To ensure the accuracy and efficiency of the simulation, the reasonable initial particle  shown to be diverge.It may be concluded that, for the droplet impacting the elastic-plate case, both terms of the artificial viscosity should be retained to keep particles in order which improves the numerical stability.Therefore, the conventional choice of  Π = 1 and  Π = 2 is suitable for the model in this paper.

The Process of Impacting.
In this section, results from two-dimensional simulations of a droplet impacting on an elastic plate are shown to demonstrate the performance of the proposed SPH method for the fluid-solid coupling problems.Figure 3 shows the results at different time during the impacting.The stress tensor and velocity field of the process are shown in Figures 4 and 5, respectively.The shape of the droplet and the deformation of the plate at different dimensionless time are presented in Figures 3(a) and 3(b), respectively.It can be seen that the process is divided into three phases.In the first phase ( ⩽ 1.75), the droplet retains its rough shape and the deformation of the plate is not obvious.The droplet touches the plate at the nondimensional time  = 1.33.In the second phase (1.75 <  ⩽ 3.0), the lengths of the droplet and plate in -direction increase rapidly.The shape of the droplet changes obviously and the impacting surface of plate becomes slightly concave due to the elasticity of the solid.In the third phase ( ⩾ 3.0), the deformation of the droplet and plate changes with time slowly.
The stress tensor given here is   and   because the changing trend of   is the same as   .It is understood from Figure 4 that the normal stress tensor   and the shear stress tensor   increase rapidly at the initial stage of impacting.Moreover, the shear stress tensor at the edge of droplet is relatively large and symmetric.In the second phase, the stress tensor of the droplet and plate decreases rapidly when the droplet is deformed significantly in -direction.Finally, the stress tensor decreases gradually and goes to zero when the shape of droplet has little change.With reference to Figure 5, it can be seen that the fluid velocity in -direction   of the contact point decreases sharply but the velocity in -direction   increases rapidly when the droplet touches the plate.And the velocity field is symmetrically disposed about a vertical axis.In the second phase, the contact area increases along with the deformation of droplet.The fluid velocity   of the contact area goes to zero, while the fluid velocity   at the edge of droplet is relatively large and symmetric.Moreover, the elastic plate stretches and deforms under the force of the impacting.In the third phase, the velocities in -direction and -direction decrease rapidly and finally go to zero.

Analysis of Results.
In order to analyze the influencing factors of the deformation degree, the process of the droplet impacting on the plate is simulated.The plate is regarded as a rigid solid and an elastic solid with different elastic modulus, respectively.The lengths of droplet and plate along -direction are shown in Figure 6.The stress tensor at the central point of droplet (the initial coordinate is (0.0 m, 0.04 m)) and the stress tensor at the central point of the plate upper surface (the initial coordinate is (0.0 m, 0.0 m)) are obtained and shown in Figure 7. Figure 8 shows the velocity of these points along -direction during the impacting process.

Comparison with the Results of Other Methods.
The deformation of droplet in the process of impacting on a rigid wall has been predicted with grid-based method by Tomé et al. [26] and with a meshless method by Fang et al. [31].The former used the FDM method, while the latter used the SPH method which is roughly the same as the method adopted in this work.However, the main difference is the treatment of the boundary condition.This problem is a critical factor in SPH method because it determines the algorithm on the fluid-solid interface.Fang et al. used two types of virtual particles to model the rigid boundary conditions.The virtual particles of the first type were located right on the solid boundary, while the virtual particles of second type were called "image particles" to solve the particle deficiency problem near the boundary.And both of the two types of virtual particles behaved entirely similar to the real fluid particles.They had the same material properties and contribute to the usual SPH expressions for velocity, pressure, and stress gradients.Therefore, the solid plate was not the real solid and the two-phase problem had been converted into the one-phase problem.In this paper, a new algorithm based on SPH method is provided for the fluidsolid interface boundary.In contrast to the work by Fang et al. [31], the particles located on the solid plate in this work carry properties of solid material, such as density, the elastic modulus, and Poisson's ratio.Here, in order to resemble the no-slip condition on the rigid-plate boundary, velocities of the plate particles are set to zero at the end of each time-step.The computed width of the droplet is compared with those of literatures [26,31] in Figure 6(a) to verify the feasibility of present method.It was found that the results are in good agreement with those of the other two methods at the initial stage of impacting.And the differences between the results obtained by different methods can be associated with the difference in treatment of the fixed boundary.In the SPH proposed by Fang, the velocities and stresses of the boundary virtual particles were obtained according to the distance between real particles and virtual particles (for further details, see [31]).However, the velocities of the boundary grids in the FDM and the velocities of the boundary particles in the improved method are set to zero.And the stresses of the boundary are calculated from the governing equations.Therefore, the differences between the results of the FDM and the improved SPH shown in Figure 6(a) are small.

The Effect of the Elastic Modulus.
As shown in Figure 6, the width of droplet and plate in -direction increases during impacting, which means that droplet and plate begin to deform.It can be seen from Figure 6(a) that the elastic modulus of solid affects the deformation degree of droplet.The elastic modulus has slight influence on the deformation of droplet at the initial stage of impacting, but the effect begins to emerge in the second and third phase.As a whole, the width of droplet decreases with the increase of the elastic modulus.Obviously, the width of droplet is the smallest when the plate is regarded as a rigid solid, because the solid with the larger elastic modulus tends to maintain the rigid character.In Figure 6(b), the length of plate increases rapidly because the elastic plate can absorb some of the impact.In the first and second phase, the length of plate increases with the increase of the elastic modulus.However, the smaller the elastic modulus is, the faster the speed rate of the length grows.And finally the plate with the smallest elastic modulus has the greatest deformation.The result of simulation was in accordance with the facts that the solid with the smaller elastic modulus is more flexible and its load bearing capacity is lower.Figure 7 shows the normal stress tensor   of the central point for droplet and the central point on the plate upper surface during the impacting process.It is shown in Figure 7(b) that the normal stress tensor of the plate upper surface reaches the maximum when the droplet touches plate at the nondimensional time  = 1.33.In addition, the maximum of the normal stress tensor for plate decreases with the decrease of the elastic modulus.The reason is that the plate with the smaller elastic modulus is easier to relieve impact efficiently and its load bearing capacity is lower.Meanwhile, the central point of droplet does not touch the plate, so the normal stress tensor   of this point increases during impacting and reaches the maximum value when this point of droplet impacts the plate (shown in Figure 7(a)).And then the normal stress tensor decreases intensely and has a period of fluctuation because the elastic plate acts as cushioning.Eventually, the normal stress tensor tended to be stable and can be ignored.
The velocities of droplet and plate along -direction during whole deformation process are obtained and shown in Figure 8.As shown in Figure 8(a), the velocity of droplet decreases and goes to zero when the deformation of droplet and plate tends to be stable.It is shown in Figure 8(b) that the central point of plate on the upper surface begins to move downward under impact when the droplet touches the plate.And the velocity of this point increases with the decrease of elastic modulus because the plate with the smaller elastic modulus is easier to elastically deform.Finally, the velocity of plate comes down to zero after a period of fluctuation.

Conclusion
This paper presents an improved SPH method to deal with the fluid-solid coupling problems.In this work, the SPH method is modified by adding a number of features, such as the correction algorithm of density and artificial viscosity.And the results for a droplet impacting the fixed boundary obtained by this method are in good agreement with those from the FDM of Tomé et al. [26] and the other SPH method of Fang et al. [31].The correction algorithm of density proposed in this paper is used to avoid the distortion of density on the interface.Moreover, an artificial viscosity is also added to the SPH momentum equation to reduce the tensile stress instability.And it is found necessary for the impact problem to retain both terms of artificial viscosity to achieve a better stability and more accurate free surface.The results show that the proposed method is able to accurately capture both the shape of droplet and the microdeformation of elastic solid surface under the force of impacting.The influence of the elastic  modulus on the whole process is also discussed.Obviously, the deformation tendency is inversely proportional to the rigidity of the material, which is coincident with the practical law.
Additionally, further more detailed studies are needed for the flow phenomenon of the impacting droplet with different diameters and initial velocities, which will be performed more intensively in future work.

Figure 1 :
Figure1: The flow chart for the numerical simulation (the particle  is a certain particle and the particle  is the particles in the influence domain of particle ).

Figure 2 :
Figure 2: Particle positions of the droplet with different values of  Π and  Π at the dimensionless time  = 2.0.
+ 08 E = 5E + 08 E = 1E + 09 SPH with the fixed wall at present FDM with the fixed wall (Tom SPH with the fixed wall (Fang et al.  = 1,  = 2) SPH with the fixed wall (Fang et al.  = 0,  = 2) e et al.) + 08 E = 1E + 09 SPH with the fixed wall at present (b) The width of plate

Figure 6 :
Figure 6: The width of droplet and plate during the impacting process.
The normal stress   of plate

Figure 7 :
Figure 7: The normal stress tensor   of droplet and plate during the impacting process.

1 ) 1 )
+ 09 SPH with the fixed wall at present Velocity V y of drop (m•s −(a) The velocity   of droplet + 08 E = 1E + 09 SPH with the fixed wall at present Velocity V y of plate (m•s −(b) The velocity   of plate

Figure 8 :
Figure 8: The velocity   of droplet and plate during the impacting process.