Abstract
Discrete multibody interaction and contact problems and the multiphase interactions such as the sand particles airflow interactions by Aeolian sand transport in the desert are modeled by using the different kernel smoothing lengths in SPH method. Each particle defines a particular kernel smoothing length such as larger smoothing length which is used to calculate continuous homogenous body. Some special smoothing lengths are used to approximate interaction between the discrete particles or objects in contact problems and in different field coupling problem. By introducing the Single Particle Model (SPM) and the Multiparticle Model (MPM), the velocity exchanging phenomena are discussed by using different elastic modules. Some characteristics of the SPM and MPM are evaluated. The results show that the new SPH method can effectively solve different discrete multibody correct contact and multiphase mutual interference problems. Finally, the new SPH numerical computation and simulation process are verified.
1. Introduction
Smoothed Particle Hydrodynamics (SPH) is a mesh-free Lagrangian method developed in recent years. Originally, SPH was developed for astrophysics and cosmological applications in three-dimensional open space [1–4]. It has become a powerful and flexible technique to perform multidimensional hydrodynamic simulation of large number of complex nonlinear phenomena. The resolution of the method can be easily adjusted with respect to variables such as density, velocity, and acceleration. It is a robust and conceptually simple method, although it suffers from unsatisfactory performance due to lack of consistency. Because of several benefits over traditional grid-based techniques, the SPH method has been increasingly used to model fluid mechanics [5–8], even though the accuracy is significantly lower compared to the sophisticated grid-based techniques. The SPH method is also used in many solid mechanics applications due to its advantage of dealing with larger local distortion than grid-based methods [9–11]. The basic idea of the SPH method is that the state of a system (the problem) can be discretized into a set of particles. Each particle has attributes such as mass, velocity, and thermal pressure. Also particles move with the fluid flow. During computation interactions of particles are determined by an interpolation of nearby particle distribution using a special kernel function, which has a given length-scale so-called smoothing length .
Numerical modeling for discrete multibody interactions and multifield coupling using SPH method still have some challenges. The choice of the influence radius in kernel is very important as it directly affects the accuracy and efficiency of the SPH calculation. If the radius selected is too small for the homogeneous continuous body, the particles inside of influence domain will not be enough. This leads to reduction of the calculation accuracy. On the other hand, if the radius selected is too large for discontinuous discrete multibody, the premature and overdomain evaluation will occur. Also, it cannot reflect the deformation of the locality, which also affects the calculation accuracy. In particular, to evaluate impacts of discontinuous discrete bodies, dynamic contact and multifield coupling problems may also be judged by kernel estimation by changing smoothing length, so there is no need to define special fine contact area or grids. Therefore, the problem becomes how to define a suitable smoothing length in kernel function for interaction of discontinuous discrete bodies, continuous homogenous bodies, and multifield coupling. Most existing adaptive SPH calculations are based on the nearest neighbor approach [12–14], in which broader kernels are used in regions of lower density. Determining the smoothing length according to particle’s density may not be able to effectively avoid interactions of no contacting particles from different bodies, because in the contact problems density differences of particles may be small. Most of those methods calculate each particle’s smoothing length in every time step, which may influence the efficiency of the calculation.
The windblown sand problem in desert regions is one of the typical multiphase interaction problems. In desert regions, the movement of sand caused by wind has been listed as one of the major environmental problems in the world today. The most effective way to study wind-sand multiphase flow is numerical simulation and experimental measurement. Experimental measurement in blown sand physics has been a main research method, and numerical simulation methods can study mechanism and characteristics of wind-sand multiphase flow from a microcosmic perspective. Wind-sand flow is a complex, nonlinear, self-organized air-sand two-phase flow [15, 16]. Grains of sand in the air have different trajectories because of their different size, different forces, and discrete tiny multibodies. The movement form of sand particles in the wind can be roughly summarized into three types: saltation, suspension, and creep. The root of all these movements lies in air driving sand after takeoff. So the master sand particles interaction mechanism is particularly important [17–19]. For the airflow field, there are also very mature mesh-based numerical computation methods.
However, in the mesh-based methods such as Finite Difference Method (FDM), Finite Element Method (FEM), and Boundary Element Method (BEM), because of mesh tangling and distortion it is hard to track gear transmission process, the trajectory tracking of sand particles, and so forth, especially when calculating large deformation and crack propagation. In this case computationally extensive remeshing is needed in calculation process, which affects accuracy. Coauthors Imin and Geni [20] presented different kernel smoothing length SPH method and used the method in gear mesh impact problem. Based on this research, in this paper we further validate the method by including more test models. Validation of method is performed by exploring its ability to reproduce an elastic impact through the collision of two rigid particles and an inelastic impact by making two rigid clusters of particles collide. Then in this paper, in order to solve problems such as the discrete sand particles interactions by Aeolian grain transport in the desert, each particle defines a different kernel smoothing length. For example, larger smoothing lengths are used to calculate continuous homogenous body, and some special smoothing lengths are used to calculate interactions between the discrete particles contact and the different field coupling. As a result it can effectively solve different bodies and field mutual interference problems. Finally, numerical calculation and simulation are used to verify the effectiveness of the proposed method. Results show that the method can be extended to deal with multiphase interactions problems.
2. Fundamentals of the SPH Method
In the SPH method the first step is to approximate a function and its gradient using integrals of kernel function based on interpolation method. Continuous partial differential equation is transformed into integral equations. In the second step, continuous forms of integral equations are discretized into discrete equations using the particle approximation method [6].
For any function , the value at a point can be approximated aswhere denotes kernel approximation of the function and are position vector and is a kernel function, the distance , and the smoothing length , which defines the influence area of kernel function.
Using integration by parts, Gauss theorem, and property of kernel function, gradient can be approximated asThe function and its gradient and integral equations (1) and (2) can be discretized using the SPH method as follows:where and are particle indices, and , respectively, denote mass and density of the particles , is the number of particles in the influence domain of particle , and , , . In the above equation, the sum is calculated only in influence area of instead of the entire computational area. The influence area is determined by the radius of kernel function, as shown in Figure 1 in a 2D example.
3. SPH Discretization
The governing equations for dynamic fluid can be expressed by Navier-Stokes equations (the equations for the rates of change of density, velocity, energy, and position) as follows.
The governing equations for dynamic fluid can be expressed by Navier-Stokes equations (the equations for the rates of change of density, velocity, energy and position) as follows. Conservation of mass is as follows: Conservation of momentum is as follows: Conservation of energy is as follows: Position equation is as follows:where and are used to denote the coordinate directions, is the velocity, is the density, is the total stress tensor, is the external force, such as gravity, is the internal energy, and is the position vector.
In the above equations the total stress tensor consists of two parts: the isotropic pressure and the viscous stress and is expressed asIn Newtonian fluids, the viscous shear stress should be proportional to the shear strain rate denoted by through the dynamic viscosity and written aswhere .
So the energy equation can be rewritten as
In solid mechanics, the constitutive model, in general, permits the stress to be a function of strain and strain rate. For the anisotropic shear stress, if the displacements are assumed to be small, the stress rate is proportional to the strain rate through the shear modulus :where .
By substituting the SPH approximations for a function and its derivative using (3), a set of SPH particle discrete Navier-Stokes equations can be obtained as Discretized equations of energy equation for fluid and solid can be written asNote that by using different numerical approximation it is possible to get different forms of SPH equations for the same partial differential equation.
4. Multibody and Multiphase Coupling
In SPH method, the interaction and the contact area are determined by kernel estimation, so there is no need to define special contact areas or particles. As a result, the choice of the influence radius of kernel is very important because it affects the accuracy and efficiency of the calculation directly. When calculating the continuous multibody interaction, such as gear meshing and ball bearing transmission, if the selected radius of kernel is too small, there are not enough particles in the influence domain. This causes reduction of accuracy. On the other hand, if the radius selected is too large, approximation will not reflect locality of the deformation; this also affects accuracy. In addition, when calculating the discrete tiny granular multibody interactions, such as sand particles interactions, if the selected radius of kernel is too large, there are many particles in the influence domain; interaction will occur between two particles before substantive contact. The similar issue occurs in multiphase coupling problems such as the lubrication with inclusion and ball bearing interactions or airflow and sand-flow interactions in the desert. In particular, when a pair of continuous bodies contacts each other, interaction of particles that belong to different bodies will happen near the contact surface and will be taken into account in computation. In this case, if the influence radiuses of all particles are the same, the interaction will occur early between two bodies due to the particles in the kernel radius. This causes that the two bodies always cannot directly contact each other on the surface; hence there remains a gap during the process.
The more appropriate method to solve this problem is to introduce different kernel radiuses. So in this paper different kernel radiuses are defined by choosing different smoothing lengths to compute multibody and multiphase interaction problems. The basic idea of this method is to choose large smoothed length for evaluating the continuous body and to choose small one for evaluating the discrete body as shown in Figures 2(a) and 2(b).
(a) Continuous body
(b) Discrete body
The most frequently used kernel function in the SPH method is the cubic B-spline kernel [21], so in this study the same kernel function is used. The typical form of the cubic B-spline kernel is shown in the following:where equals , , and in one-, two-, and three-dimensional cases, respectively, is defined as , and is the distance between the two points and . Usually the smoothing length is set to a constant value in the computation domain. In this study different smoothing lengths are introduced to evaluate different phenomena of multibody and multiphase interactions. It is clear from (14) that the kernel function depends on the smoothing length ; however, the relative distance is not changing with in the evaluation domain . This is a great advantage in order to satisfy the consistency condition [21] for kernel function with different smoothing length .
Figures 3(a) and 3(b) show the images of cubic B-spline kernel and derivatives when choosing different smoothing lengths of , , , and , respectively. It is obviously seen from Figures 3(a) and 3(b) that adjusting the smoothing length in the corresponding domain can reduce or increase the influence of kernel function ; therefore, it maintains consistency conditions of kernel function in all cases. A larger value of kernel radius will create averaged interactions of each particle in the body, so the calculating area will be bonded into stronger continuous body as shown in Figure 2(a). The small value of kernel radius will create very local interactions so the calculating area will be separated into discrete pieces more easily as shown in Figure 2(b). Meanwhile the consistency conditions are also satisfied.
(a) Kernel function
(b) Derivative of kernel function
The different kernel radiuses are defined by choosing different smoothing lengths to calculate continuous multibody interaction problem as shown in Figure 4(a) [20] and multiphase interaction problem as shown in Figure 4(b). The following expression is used to correctly evaluate two different bodies’ interactions:The diagonal terms , are large kernel radiuses used to evaluate the interaction of particles inside the same body and represent strong continuous bodies. The coefficient is determined by corresponding kernel function used in calculation and the smoothing length is equal to 1.2 times of the distance between two particles [22]. Whereas the off-diagonal terms , are smaller kernel radiuses used to calculate and evaluate the interaction between different bodies at contact surfaces when two bodies are not continuous and separated from each other. Generally we choose an appropriate smoothing length which keeps kernel radius equal to 1.0 times of the distance between two particles.
(a) Multibodies interactions
(b) Multiphase interactions
For two different phases the following expression is used to evaluate of airflow field as fluid phase and sand-flow field as granular phase interaction. In this study to compute multiphase flow problem such as sand-driving flow as shown in Figure 4(b), three different kernel radiuses are used as defined in the following:The term is an air-to-air particles smoothing length. Due to the airflow continuity a large kernel radius is used to evaluate the interaction of air-to-air particles and represents strong continuous airflow field. In this study, the smoothing length is set to 1.5 times of the distance between two particles. The terms and represent the sand-to-air and air-to-sand particles smoothing lengths, respectively. Because the airflow field continuously involves the sand particles, the smoothing length in this study is 1.25 times of the distance between two particles. Lastly, is a sand-to-sand particle smoothing length. Considering the sand particles are discrete with each other, a smaller kernel radius is used to evaluate the interaction of sand particles with each other and represents discrete tiny sand particle field. Therefore, the smoothing length in this study is set to 0.55 times of the distance between two particles. The coefficient is a constant value determined by corresponding kernel function.
5. Validation of New Method
For validation of this new method, different kernel smoothing lengths with different elastic modules and different particle models are used. The two kinds of influence radius are explored, one is original model by setting in all cases as shown in (17a) and another is a new model with and the as shown in (17b):
The cubic B-spline for kernel function is used in approximation with . The material is assumed to be an elastic hard material and its properties are set as follows: Poisson ratio and density kg/m3 and the elastic modulus is set as GPa. for elastic impact model and GPa. for rigid or perfectly elastic impact model. Then two main special particle models, Single Particle Model (SPM) and Multiparticle Model (MPM), are obtained as shown in Figures 5(a) and 5(b) to evaluate multibody elastic impact phenomena. The initial conditions of SPM and MPM are identical such that object A with an initial velocity 10 m/s impacts object B in quiescent state as shown in Figures 5(a) and 5(b). The MPM consists of multiparticles. It is also constructed by square models. The MPM square models consist of 1331 particles.
(a) SPM
(b) MPM square model
Figure 6 shows the results of the velocity exchange in the SPM and MPM by using original SPH method with the same smoothing length. It is shown that the velocity is exchanged nearly half between objects A and B and they stick together and move with the same velocity in both SPM and MPM square models when GPa. Additionally, in the SPM the velocity exchange process is very close to the one in the MPM square model. So from perspective of velocity exchanges on impact behavior, the single particle shape is close to the multiparticle square shape and represents a perfect inelastic impact. However, in SPM and MPM two objects A and B earlier interact with each other; this may cause some error in interaction for discrete bodies.
(a) SPM
(b) MPM square model
Figure 7 shows the velocity exchange results of SPM and MPM square models by using different smoothing lengths when GPa. It is shown that the velocity is exchanged nearly half between objects A and B in the SPM as shown in Figure 6(a), but in the MPM square model there are some differences with Figure 6(b). The velocity of A is smaller than B, and they did not stick together. This represents an elastic collision. The other benefit of this new method by using different smoothing lengths is that the exact interaction between two objects A and B is improved in both SPM and MPM for discrete bodies.
(a) SPM
(b) MPM square model
It is obviously seen from Figures 6 and 7 that by using different smoothing lengths the exact interaction between two bodies is improved and realized. Moreover, the SPM represents perfect inelastic impact, while the MPM square model represents weak elastic impact. However, in real particles such as sand or inclusion it always shows strong elasticity when impacting each other. If the MPM models are used in SPH method, the numerical simulation requires too much time and memory due to the increasing of the particle number. So the challenge is to represent the perfect elastic phenomenon of tiny hard particle by using the single particle. For this purpose, the influence of elasticity such as elastic module is explored.
Figures 8(a) and 8(b) show the results of impact velocity exchanging processes by using the original SPH method with the same smoothing length and the new SPH method with different smoothing lengths, respectively. The results show that the velocity exchanges between two objects A and B are increased with increasing elastic module in the SPM in both the original and the new SPH methods. The new SPH method gives more perfect elastic collision behavior and represents the correct impact process. The similar behavior is obtained in the MPM model as shown in Figures 9(a) and 9(b). Based on these numerical results it is clear that increasing elasticity module can increase the perfect elastic impact characteristics by using the single particle, and the multibody and multiphase coupling such as gear teeth correct meshing and contact, tiny sand-particles impact, and air-sand two-phase interaction are improved by using the different kernel smoothing length.
(a) Original SPH method
(b) New SPH method
(a) Original SPH method
(b) New SPH method
6. Application in Wind-Sand Multiphase Flow Numerical Simulation
Coauthors Imin and Geni [20] used different kernel smoothing length SPH method in the gear mesh impact problem. Based on these researches, the method was firstly applied in the numerical simulation of gas-solid two-phase flow in this paper. Different smoothing length was adopted in the different phase of wind-blow sand flow, so the interaction between the two phases can be accurately simulated. By setting the smoothing length and the distance between the two particles to judge whether a collision between the sand has taken place and solving speed vector after impact, the change of sand motion after collision can be accurately reflected without rebuilding the collision model.
The two kinds of models for wind-sand multiphase flow numerical simulation are generated as the flat sand bed surface model and the small dune models as shown in Figures 10 and 12. To evaluate two different phases of wind flow field and sand flow field, (16) is used and three different kernel radiuses are chosen to represent air-to-air, air-to-sand, and sand-to-sand multiphase interactions as mentioned above.
(a) Step = 20
(b) Step = 5400
As shown in Figure 10(a), the upper layer of sand begins to produce a lift-off velocity and trends to jump under the action of airflow on the stillness of the sand bed. Meanwhile, the top of the sand remains at the same velocity and direction. Due to the increase of the time step, the sand leaps from the bed surface as the velocity at the top of sand gradually increases; then the velocity changes while the direction basically remains unchanged. At this time step, the grains of sand in some areas at the lower layer of sand bed produce velocity because of the impact by airflow and other sand. As shown in Figure 10(b), when the time step proceeds to 5400, most of the sand in the sand bed has velocity, which means the sand bed begins to loosen or accumulate. At this time step, the top layers sand begins to have a different saltation length and height caused by the top of sand bed. The velocity and direction of the sand grains have a very large difference. As shown inside the box of Figure 10(b), some grains of sand in the sand bed have a trend to jump whereas some sand particles change the direction down to the sand bed as shown in the circled area of Figure 10(b). These particles can cause sand bed to be more compact. One of the direct causes of sand ripple is the loosening or tightening of the sand bed.
Figure 11(a) shows the changes of velocity and position of sand particles at the top of the sand bed in the process of takeoff. Initially, the sand particles at the top of the sand bed are standstill, as shown in Figure 10(a). Then the sand particles begin to generate velocity by the acting force of the airflow; the magnitude and direction of the velocity have little changes at this moment. Only after it reaches certain velocity, the sand particles jump from the sand bed surface.
(a) Step = 0~12600
(b) Step = 18000
(c) Step = 21900
While the airflow acts on the sand particles, the sand also brings certain reactive force to the airflow. This reactive force causes the airflow velocity to change near the grains of sand. Changes of airflow velocity caused large differences in velocity and direction of sand particles at the same layer; see Figure 10(b). The velocity and displacement of sand at the same layer would change gradually, and sand particles move at different velocities and produce different trajectories as shown in Figure 11(a). At this time step, the originally flat sand bed surface morphology such as the sand accumulation and erosion substantially changes; this shows in the form of sand ripple on the bed surface.
Figures 11(b) and 11(c) show the magnitude and direction of velocity of sand particles at the top of the sand bed in step = 18000 and 21900, respectively. At steps 18000 and 21900, the most of jumping sand have already started to come down to the bed surface.
Initial Condition. Sand particles are quiescent at the beginning; the initial airflow speed is defined as follows [23]:
Boundary Condition. SD and ST are fixed boundary conditions; SL and SR are periodic boundary conditions. SI is the boundary between airflow and sand. A row of ghost particles are set on the bottom boundary SD and the top boundary ST. These particles have a repulsive force to other particles to prevent particle penetration boundary.
Figures 13(a) and 13(b) show the velocity distribution of airflow field at the time steps 0.2 ms and 70 ms, respectively. It is shown that the airflow velocity changes when the airflow is close to the sand dune area and airflow speed increases gradually along the slope of windward side. After the airflow throughout the sand dune leeward side, the speed decreases and the direction gradually return to the initial state.
(a) ms
(b) ms
Figure 14 shows that the sand particles’ speed increases with the interaction of airflow. At the beginning, the sand particles on the top of sand dune windward site get larger speed and the speed causes jumping movement. In the process of sands saltation, the sands move along the airflow direction meanwhile sinking motion occurs because of the effect of sand’s self-gravity. Eventually, the sand saltation system reaches equilibrium.
(a) ms
(b) ms
7. Conclusion
In this paper two different models are introduced to evaluate the velocity exchanges during the impact process. It is shown that in the SPM the velocity exchange during the impact process is very close to the one in the MPM square model, so from perspective of velocity exchange on impact the single particle shape is close to multiparticle square shape and it represents the perfect inelastic impact phenomenon. However, the two objects A and B in the SPM and MPM influenced each other at the very early stage of the process; this may cause some error in interaction for discrete bodies.
The different kernel radiuses are defined by choosing different smoothing lengths to approximate continuous and discrete multibody interactions as well as multiphase interactions. The new SPH kernel smoothing length is represented and modeled for evaluating the multibody correct contact and multiphase interaction problems.
By using the new SPH kernel smoothing length, the results show that the velocity is exchanged nearly to half in the SPM, but in the MPM square model there are some differences. In this case the velocity of object A is smaller than object B. Additionally; the two objects (A and B) do not stick together. This represents an elastic collision. The other benefit of this new method by using the different smoothing length is that the exact interaction between two objects A and B is improved and realized in both SPM and MPM for discrete bodies.
The velocity exchanges between two objects A and B are increased with increasing elastic module in SPM and MPM in both original and new SPH methods. The new SPH method gives more perfect elastic collision behavior and represents the correct impact process. So the SPM represents a perfect inelastic while the MPM square model represents a week elastic characteristics. In real particles such as sand or inclusion it always shows a strong elasticity when impacting each other. If the MPM is used in the SPH method, the numerical simulation has high time and space complexity due to the increase of the particle number. It is clear that by increasing the elasticity of the SPM it can increase the perfect elastic impact characteristics. This will bring great convenience to solve the multibody and multiphase coupling such as tiny sand-particles impact and air-sand two-phase interaction with relatively lower time and space complexity.
The new SPH method is applied to the evaluation of sand-air multiphase interactions. The two kinds of models for wind-sand multiphase flow numerical simulation are applied as the flat sand bed surface model and the small dune models with the new SPH method. The results give a good agreement with real sand-air flow phenomenon and represent air-to-air, air-to-sand, and sand-to-sand multiphase interactions.
Conflict of Interests
The authors declare that there is no conflict of interests regarding the publication of this paper.
Acknowledgments
The authors wish to acknowledge the support of the China natural science foundation (no. 51075346) and National Key Basic Research and Development Program (973 Program no. 2011CB706600).