Recent Advances in Solution Methods for Nonlinear Evolution Equations, Fluid Flow, and Heat and Mass TransferView this Special Issue
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 SPH method, 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.
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 , MAC , VOF , and LS  methods. Moreover, mesoscopic simulation methods have been proposed including MPS , LBM , and DPD . 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  and Gingold and Monaghan  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 , non-Newtonian fluids flows [11–14], incompressible fluids [15–17], multiphase flow [18–20], free surface flows [21, 22], large deformation , and the dynamic response of elastic-plastic materials [24–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.  used the finite difference method to simulate unsteady free surface flows. Using MAC, the problem of impacting drops was simulated by Harlow and Shannon [27, 28]. Xie et al.  treated the problem of flow phenomenon after droplet with different diameters and velocities impacting the fluid films by using MPS. Prosperetti and Oguz  predicted the changing of free surface by BEM when droplet impacted the fluid films. Fang et al. , Rafiee et al. , and Ma et al.  simulated a single droplet impacting a rigid plate by SPH. Various modified algorithms for SPH were employed to improve stability and accuracy in these literatures. However, these methods for simulating the impacting droplet are based on the fixed boundary without considering the deformation of plate. Regarding the plate as a rigid solid surface is a feasible solution for the free surface problem when the hardness of solid is much higher than fluid. In fact, the absolutely rigid solid does not exist. The stability of system will be significantly impaired because of microdeformation of solid or fluid. Taking elastohydrodynamic lubrication, for example, the microdeformation of the wheel gear or bearing can lead to many fatal consequences in practice. Therefore, it is necessary to consider the interaction between solid and fluid to analyze the effect of deformation on the whole system.
The purpose of this paper is to extend and test the ability of the SPH method for the solid-liquid interaction problem. In Section 2, the governing equations are described in detail and the SPH formulation is modified to cope with the special issues. The modifications and improvements include the correction algorithm of density on solid-liquid interface, artificial viscosity, and artificial stress. The process of a droplet impacting an elastic plate with different elastic modulus was simulated and some results are displayed in Section 3. The paper ends with some concluding remarks in Section 4.
2. 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  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–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 . Therefore, the cubic spline function is chosen as follows:in two dimensions, the normalization factor is .
2.3. 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 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.
2.4. 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 asThe most common form of the artificial viscosity suggested by Monaghan [10, 39] is where , and denotes the sound speed. The parameter 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 , which means that smoothing of velocity will only take place if the particle spacing is less than . 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  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 and the -term should be retained to prevent unphysical particle penetration . 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.
2.4.2. 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  and Gray et al. . 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 asThe detailed formulation of the artificial stress is as follows:where , , and is the initial distance of particles. The components of the artificial stress tensor are given by [31, 32]whereIn these equations, is a parameter and it is set to 0.2 for the Newtonian fluid [31, 32].
2.5. 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.
2.5.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  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 fluid-solid 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 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.
2.5.2. The Solution of Momentum Equation
The key of momentum equation (10) is to calculate the stress tensor. 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 . ConsiderThe constitutive equation for the additional elastic stress of Newtonian fluid is as follows:The equation of state is as follows:where is the dynamic viscosity of fluid, 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 (, 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.  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. Considerwhere and 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).
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:
3. 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.
3.1. 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 m and m·s−1. The density and dynamic viscosity of droplet are set to kg·m−3 and 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 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 spacing for this model is set to cm corresponding to 10203 particles in all (7698 of which are fluid particles). For the fluid-solid coupling problems, the process of a droplet impacting an elastic plate is simulated. The elastic plate with thickness of 0.01 m is placed at m and the length of plate along -direction is m. The solid material parameters are chosen to be kg·m−3, , and the elastic modulus E ranges from to . The time step of is chosen to ensure numerical stability in simulation.
3.2. The Values of and in the Artificial Viscosity
To check on the reasonable values of and for the SPH method studied here, the particle positions of droplet with different values of and at the dimensionless time are shown in Figure 2. Here, the superscripts and are for solid andfluid, respectively. It can be seen that the -term of artificial viscosity has a significant influence on the computed solution. Without the -term () for solid, no matter what the value of for fluid is, the droplet deforms unrealistically and some fluid particles penetrate the plate (shown in Figures 2(a) and 2(b)), which consequently causes the simulations to diverge. Besides, in the case with the -term () for solid and without the -term () for fluid, there are some observable voids in the droplet and the particle positions of the droplet become disorderly (shown in Figure 2(c)). It can be seen from Figure 2(d) that the void phenomenon is avoided completely by retaining the -term () for solid and fluid. The models without the -term () have also been conducted, and the results are 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 and is suitable for the model in this paper.
3.3. 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 (), the droplet retains its rough shape and the deformation of the plate is not obvious. The droplet touches the plate at the nondimensional time . In the second phase (), 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 (), 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.
3.4. 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.
(a) The width of droplet
(b) The width of plate
(a) The normal stress of droplet
(b) The normal stress of plate
(a) The velocity of droplet
(b) The velocity of plate
3.4.1. 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.  and with a meshless method by Fang et al. . 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 fluid-solid interface boundary. In contrast to the work by Fang et al. , 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 ). 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.
3.4.2. 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 . 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.
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.  and the other SPH method of Fang et al. .
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.
Conflict of Interests
The authors declare that there is no conflict of interests regarding the publication of this paper.
This work was supported by the National Basic Research Program of China (973 Program, no. 2011CB706602), China National Natural Science Foundation (no. 11072209), and Open Projects of State Key Laboratory for Strength and Vibration of Mechanical Structures (Xi’an Jiaotong University).
J. E. Welch, F. H. Harlow, J. P. Shannon, and B. J. Daly, The MAC Method: A Computing Technique for Solving Viscous, Incompressible, Transient Fluid Flow Problems Involving Free Surfaces, Los Alamos Scientific Laboratory, Los Alamos, NM, USA, 1966.
J.-Y. Deng, J. Yang, L.-Q. Wang, H.-X. Du, and N. Wang, “Simulate of liquid-solid two-phase flow by smoothed particle hydrodynamics method,” Journal of Zhejiang University (Engineering Science), vol. 41, no. 5, pp. 853–858, 2007.View at: Google Scholar
L.-Q. Ma, M.-B. Liu, J.-Z. Chang, T.-X. Su, and H.-T. Liu, “Numerical simulation of droplet impact onto liquid films with smoothed particle hydrodynamics,” Acta Physica Sinica, vol. 61, no. 24, Article ID 244701, 2012.View at: Google Scholar
G. R. Liu and M. B. Liu, Smoothed Particle Hydrodynamics: A Meshfree Particle Method, World Scientific, River Edge, NJ, USA, 2003.
I. P. Natanson, Theory of Functions of a Real Variable, Ungar, New York, NY, USA, 1960.
Q. Hongfu and M. Lijun, “Computational simulation of hypervelocity penetration using adaptive SPH method,” Transaction of Tianjin University, vol. 12, pp. 75–78, 2006.View at: Google Scholar