Numerical Simulation of the Generalized Newtonian Free Surface Flows by a Density Reinitialization SPH Method

A periodic density reinitialization smoothed particle hydrodynamics (PDRI-SPH) method is proposed to treat the generalized Newtonian free surface flows, which is based on the concept of Taylor series expansion. Meanwhile, an artificial stress term is also presented and tested, for the purpose of eliminating the unphysical phenomenon of particle clustering in fluid stretching. The free surface phenomena of a Cross model droplet impacting and spreading on an inclined rigid plate at low impacting angles are investigated numerically using the proposed PDRI-SPHmethod. In particular, the effect of the surface inclination and the different regimes of droplet impact, spreading and depositing on an inclined surface, are illustrated; the influence of surface inclination on the tensile instability is also concerned. The numerical results show that the accuracy and the stability of the conventional SPH are all improved by the periodic density reinitialization scheme. All numerical results agree well with the available reference data.


Introduction
The problems of free surface flows for polymers are important in today's industry, such as the structured reactors, surface coating, container filling in the food, and pharmaceutical industries of polymers.All the flows involved almost exhibit nonlinear behavior, for example, the viscoelastic or shearthinning behavior.In these processes, the impacting, spreading, and depositing of liquid droplets on solid surface play a crucial role.
In the early stage of research, many methods based on the Eulerian description of motion are mainly presented to capture the complex free surface of polymers, including particle in cell (PIC) [1], marker and cell (MAC) [2], volume of fluid (VOF) [3], level set [4], and phase-field [5] methods.These methods are based upon grid-based numerical methods such as finite difference methods (FDM) and finite element methods (FEM) that are commonly used to solve the Navier-Stokes equations.However, it is difficult for the simulation of large deformation.
In order to overcome the shortcomings of grid-based methods and effectively handle the problem of large deformation, the various mesh-free methods [6][7][8] or particle methods have been proposed in a Lagrangian framework.Among the various particle methods [9][10][11][12], the smoothed particle hydrodynamics (SPH) method [9,13] is the earliest one and it is also most widely used.The SPH has the following main advantages over grid-based methods.(1) It handles convection dominated flows and large deformation problems without any numerical diffusion.(2) Complex free surfaces are modeled easily and naturally without the need of explicit surface tracking technique.(3) It is easy to program for complex problems compared with grid-based methods.In 1994, it was firstly used to deal with fluid mechanic problems [14].Since then, it has been extensively studied in many areas such as viscous flows [15,16], incompressible fluids [17,18], multiphase flows [19,20], geophysical flows [21,22], viscoelastic flows [23,24], and viscoelastic free surface flows [25].
Unfortunately, the consistency between mass, density, and occupied area cannot be enforced exactly (see [15,26]) when the evolved particle density is obtained by the continuity equation in standard SPH method [23][24][25].In this work, a periodic density reinitialization method based on 2 Mathematical Problems in Engineering the corrective kernel estimate [27] of a Taylor series expansion is proposed to overcome the problem of the consistency between mass, density, and the occupied area.Moreover, we can know that the tensile instability is related to the sign of both the stress and the second derivative of the kernel function as noticed by [28].And then, a changing artificial stress term is presented and tested to remove the unphysical phenomenon of particle clustering in the simulations of a generalized Newtonian droplet impact and spreading on an inclined rigid plate, which is different from the one in [25].
This paper has been directly motivated by the polymer industry where materials tend to be shear-thinning but not necessarily viscoelastic.Due to the fact that the Cross model [29] or some similar model can describe the shearthinning behavior better, then here we choose the Cross model.In general, the phenomena of free surface can be complex.Therefore, the two-dimensional shear-thinning free surface flows of a Cross droplet impact and spreading on an inclined rigid plate are discussed.During these processes, the effect and ability of the mentioned above periodic density reinitialization method and artificial stress term for capturing the complex polymer free surfaces are also analyzed.
The structure of this paper is organized as follows.The governing equations for the Cross model are introduced in Section 2. Section 3 describes the PDRI-SPH discretization of the Navier-Stokes equations, including artificial viscosity, boundary conditions, and temporal discretization of the governing equations.In particular, the density reinitialization method and tensile instability are also discussed in Section 3. In Section 4, the validity of the proposed PDRI-SPH combined with the mentioned artificial stress term is first tested.Subsequently, a numerical example of a Cross model droplet impacting on an inclined dry surface is solved to demonstrate the capability of the PDRI-SPH method in handling generalized Newtonian free surface flows.Some concluding remarks are reported in Section 5.

Governing Equations for the Cross Model
2.1.Governing Equations.In a Lagrangian frame, the generalized Newtonian fluid is governed by the conservation of mass and momentum equations, together with a nonlinear constitutive equation.The isothermal, incompressible fluid is usually described by the following equations: where  denotes the fluid density, V  the th component of the fluid velocity,   the (, )th component of the total stress tensor,   is the th component of the gravitational acceleration, and the / = / + V  ⋅ (/  ) is the material derivative.The spatial coordinates   and time  are the independent variables.
The total stress tensor in ( 2) is commonly made up of the isotropic pressure  and the components of extra stress tensor   : where   = 1 if  =  and   = 0 if  ̸ = .In order to study a shear-thinning polymer material, the relating constitutive equation must be provided.

Cross Model.
In this paper, the surface problem based on the generalized Newtonian fluid is mainly considered.The generalized Newtonian fluid displays more complex fluid characters than the Newtonian fluid, and the constitutive model for describing the generalized Newtonian fluid can be derived from the Newtonian model; that is, the viscosity is variable for the generalized Newtonian fluid.Several constitutive models for describing generalized Newtonian fluid have been proposed, in which the Cross fluid model with four parameters (see ( 6)) [29,30] is usually used to represent the polymer material and describe the shear-thinning behavior during polymer processing.Then the typical constitutive model of Cross fluid with four parameters is employed to study the influences of shear-thinning behavior on the free surface in polymers impact process.It is worth noting that the surface tension is not considered in the following simulations.The measure of fluid droplet is centimeter level in the following simulations, so that the effect of surface tension can be omitted according to the physical knowledge.In addition, the detailed description of Cross model with four parameters can also be found in [29,30].
The extra stress tensor  for the generalized Newtonian fluid based on the Cross model [29,30] is expressed as where ( γ ) = ( γ ) is the dynamic viscosity and the ( γ ) is the kinematic viscosity.The symmetric strain rate tensor d is defined as The kinematic viscosity ( γ ) represents the shear-thinning nature of the fluid; it is defined as where ,  0 ,  ∞ , and  are given positive constants.The γ is the local shear rate defined by and the symbol "tr" denotes the trace of matrix.Here, the positive constants  and  are all chosen equal to 1.0.Monaghan [14] and Morris et al. [15]).Here, the incompressible flows are also treated as weakly compressible flows using the following equation of state [24]:

Equation of
where  is the speed of sound and  0 is reference density.An artificial, lower sound speed is usually used to avoid the instability and extremely small time steps.To keep the density variation of fluid less than 1% of the reference density, the Mach number  ( ≡ /, where  is a typical reference velocity) [9] must be smaller than 0.1.In other words, the sound speed must be ten times higher than maximum fluid velocity.

PDRI-SPH Formulation
3.1.Discretization Schemes of Standard SPH.The SPH method [9,13] is based on the interpolation theory, which is the theory of integral interpolates using a kernel function.Namely, the fluid domain Ω is discretized into a finite number of particles, where all the relevant physical quantities are approximated in terms of the integral representation over neighboring particles.Each particle carries a mass , velocity k, and other physical quantities depending on the problem.Any function (r) defined at the position r = (, ) can be expressed by the following integral: where  represents the kernel function (or smoothing function) and ℎ denotes the smoothing length defining the influence area of .The kernel function  is usually required to meet three properties, namely, the regular condition, Dirac delta function property, and compactly supported condition [9,13].In addition, the smoothing function is also usually chosen as an even function over Ω.
According to (9), the integrating principle by parts and the divergence theorem, the particle discretization scheme of standard SPH for a function (r), and its first derivative at the position r = (, ) of the particle  can be written in the following condensed forms: where   and   are the mass and density of the th particle, and   = (r  ).The   /  represents the occupied volume by the th particle.The   = (|r  − r  |, ℎ),   /r  = −  /r  .In this paper, the cubic spline function is chosen as the smoothing function which is the function about   = |r  − r  | and  =   /ℎ.Then it reads for 2-D as follows: In order to have an accurate interpolation, the smoothing length ℎ should be chosen bigger than the mean interparticle distance.Here, the smoothing length ℎ is given by 1.5 0 with  0 as the initial distance between neighboring particles.The compact support domain size is 2ℎ.
Considering the discrete gradient equation (11) and the following identity: (1/)(  /  ) = (  /)/  + (  / 2 )(/  ), the particle discretization schemes of the governing equations can be obtained at the particle : Introducing the velocity gradient then the particle approximation schemes of constitutive equation for the Cross model can be obtained as

Density Reinitialization Method.
In the standard SPH method, each particle has a fixed mass.If the number of particles is constant, mass conservation is intrinsically satisfied.However, the consistency between mass, density, and the occupied area could not be enforced exactly (see [15,26]) if the evolved particle density is determinated by the evolution equation ( 1) for simulating the weakly compressible flows.Although the density field is periodically reinitialized by applying the following equation for removing this problem in [15]: the particle approximations (17) do not have first order accuracy and  0 ,  1 consistency for boundary regions or irregularly distributed particles (see [31]).Therefore, the above periodic density reinitialization method (see ( 17)) cannot well alleviate the above problem.
In order to well overcome this problem, we use a secondorder accurate particle approximation scheme based on Taylor series expansion (see [27,31]) to periodically reinitialize the density field: where the corrected kernel function  Tay  is given by Here   =   −   ,   =   −   , and   is replaced by   /  .The particle approximations scheme (18a)-(18c) possesses  0 ,  1 consistency for boundary regions or irregularly distributed particles (see [31]).
When the density reinitialization method is applied in above standard SPH method, an inversion matrix of 3 × 3 should be solved for each fluid particle; thus, the computing time is increased slightly.Considering the computational cost and the efficiency of using periodic density reinitialization, we can apply this procedure every fixed (about 10∼40, see Section 4.1) time step in our numerical simulations.In particular, for the purpose of preserving that the corrected particle approximations (18a) at least have  0 consistency on the whole domain, the matrix ∑  (r  )  may be replaced by ∑      if the matrix of (18b) is singular (it occurs occasionally).Now, the above discretization schemes of standard SPH combined with the periodic density reinitialization method may be called the "PDRI-SPH" method.

Artificial Viscosity Model.
According to previous works [15,32], the artificial viscosity is first introduced to enhance the numerical stability and accuracy in the simulations of strong shock problems [32].On the other hand, the artificial viscosity term guarantees the conservation of the angular momentum without external force when it is added into the momentum equation of TSPH schemes.For that reason, the artificial viscosity term is usually also considered and employed in the SPH simulations of viscous or viscoelastic fluid flows problems with large deformation, which can be seen in recent works [18,25,33].Through the simulations   of viscoelastic droplet impact problem in [25] and our numerical simulation experience of using SPH or improved SPH method, we find that it is necessary to employ an artificial viscosity term in the discrete momentum equation (14) for improving the numerical stability (see Section 4.2).Here, the artificial viscosity term is also added to the discrete momentum equation ( 14) of PDRI-SPH method, which is usually chosen as [32,34] where The 0.01ℎ 2 term is included to prevent numerical divergence when two particles get too close to each other.The  Π and  Π are usually chosen approximately equal to 1.In the artificial viscosity, the first term associated with  Π involves shear and bulk viscosity, while the second term associated with  Π is similar to the von-Neumann-Richtmeyer viscosity for resolving shocks and is very important in preventing unrealistic particle penetration.

Artificial Stress Model.
In 1995, the "tensile instability" was first investigated in detail by Swegle et al. [28], which pointed out that the phenomenon of unphysical clustering of particles arises when the standard SPH method is applied to Euler problem.At present, a number of methods have been proposed to remove the tensile instability in elastic dynamics of solid materials.The artificial stress method [35,36] is one of the most successful approaches, which has been successfully extended and applied to non-Newtonian fluid free surface flows [25].In [35,36], the authors think the "tensile instability" is mainly caused by tension (positive stress in tension), so that the adopted artificial stress term [25,35,36] is only related to the positive stress.As noticed by [28], the "tensile instability" is related to the sign of both the stress and the second derivative of the kernel function, which implies that the instability is caused by not only the tension but also the compression (negative stress in compression).Therefore, we use the following artificial stress term by extending the conclusions in [28,36] to eliminate the "tensile instability": where  = (0, ℎ)/( 0 , ℎ),   =   (|r  − r  |, ℎ)/  ( 0 , ℎ).
The components of the artificial stress tensor    are given as where  is a positive parameter (0 <  < 1), and the Introducing the   =    (   +   ), the artificial viscosity term (19a) and the artificial stress term (20a) are added to the discrete momentum equation ( 14) of PDRI-SPH and we can obtain Usually, the particle positions are updated by the following equation:

Boundary Condition Treatment.
In most engineering problems, the physical boundary might be the surface of rigid bodies enclosing fluid or enclosed by fluid, fully or partially.The boundary can be stationary or in motion.We know that the treatment of boundary conditions is very important in the numerical simulation process using SPH method.Several methods for treating rigid wall boundary conditions have been presented in previous work.There are mainly two methods; that is, (1) the solid walls may be simulated  by particles, which exert repulsive force by employing an artificial repulsive force (see [14]) on inner fluid particles to prevent them from penetrating the wall.(2) The wall boundary conditions can also be modeled either by fixed particles [37] or by virtual particles that mirror the physical properties of inner fluid particles.The above methods of boundary treatment have been discussed in 2009 [38], and the literature shows that the virtual particles approach has better stability and affectivity than the artificial repulsive force method.So, the boundary particles in this work do not employ an artificial repulsive force instead of adopting the virtual particles on approaching real particles to prevent fluid particles from penetrating rigid walls.
As shown in Figure 1, two types of virtual particles are used to implement the boundary conditions on a rigid wall.The first type virtual particles are located right on the rigid wall, namely "wall particles." The density of wall particles is not evolved unlike Morris et al. [15].Meanwhile, the nonslip condition is enforced on the solid wall and the positions of wall particles remain fixed in time.If the no-slip condition was not considered in the simulations, the fluid particles may penetrate the wall and the numerical simulations will be terminated.The pressure on the wall particles is calculated according to the following approximation formulation: where  represents the index of a wall particles and  denotes the index of its neighboring fluid particles only.
The second type virtual particles are placed just outside the solid wall and fill a domain with at least a range of depth comparable with the compact support of the kernel used in the computations, which are called "ghost particles" and have fixed density and positions.The velocity and the pressure on the ghost particles are computed in the following way.
(1) For each ghost particle , we assign a corresponding point  just on the rigid wall and point  inside domain, respectively.Meanwhile, these three kinds of points lie in a line which is perpendicular to the wall.
(2) In order to calculate conveniently we can define the normal distances   and   of the points  and  to the rigid wall, respectively.
(3) The pressure   and velocity V  I for the ghost particles can be obtained through the following linear extrapolation: where L represents the vector of variables (, V  ).To specify the values for L  , the interpolation formulation (23a) is applied again.Here, we let the   =  0 .
Moreover, the following total stress-free condition must be satisfied in the computational domain for surface particles: where n denotes a unit normal vector to the surface.In this paper, the surface tensor is neglected and a Dirichlet boundary condition of zero pressure is given to the surface particles.This condition, that is, (24), is satisfied naturally by the PDRI-SPH method.
3.6.Time Integration Scheme.In order to better illustrate the effect of the density reinitialization method, a suitable time integration schemes is chosen necessarily in practice.
Considering that the predictor-corrector scheme possesses second-order accuracy and better stability, we chose the predictor-corrector scheme for solving the system of ordinary differential equations ( 13), (21), and (22).The predictor step consists of an Eulerian explicit evaluation of all quantities for each particle where X  represents the vector of the unknown variables (  , V   ) and Γ  denotes the vector of right-hand sides of ( 13), (21), and (22).In the corrected step, the updated value of X  at the end of each time step is given by To ensure the numerical stability, the time step and space step must satisfy the well-known Courant-Friedrichs-Lewy (CFL) condition.According to [15], we may choose the following stability condition: where   is the hydrodynamical force acting on the particle, and  0 = / 0 is the kinematic viscosity.

Test Examples and Numerical Simulations for the Cross Model
Firstly, the effect and the validity of modified models for simulating generalized Newtonian free surface flows using the PDRI-SPH method are obviously demonstrated through applying the periodic density reinitialization method, artificial viscosity model, and artificial stress model.Subsequently, the capacity of PDRI-SPH method for solving a Cross model droplet impact, spreading and depositing on an inclined rigid plate, is shown in Section 4.3.

Effect of the Periodic Density Reinitialization Scheme.
In order to show the effect of PDRI-SPH method and compare with the conventional SPH method, the two-dimensional benchmark problem of the stretching of an initially circular water drop is simulated using PDRI-SPH and SPH, respectively, without enforcing the artificial viscosity, artificial stress, and rigid boundary condition.This example has been used in the literature [14,33] of using SPH method, and its corresponding analytical solution can be obtained from [14,33].
All the physical quantities are the same as those in [33, Figure 1], the reference density  0 = 1000 kgm −3 , the viscosity  = 0.001 kgm −1 s −1 , and the speed of sound  = 1400 ms −1 .The initial geometry of the water drop is a circle of radius  = 1 m with its center located at the origin ( = 0,  = 0).There is no external forces but initial velocity field The number of fluid particles is 1961 and corresponding to the initial distance  0 = 0.02 m, and the time step  = 10 −5 s.During the stretching process of water drop, the water drop remains elliptical shape and the value of ã ⋅ b (ã is the semiminor axis and b is the semimajor axis) remains constant.We let "num" denote the interval time step of the periodic density reinitialization of PDRI-SPH, and the PDRI-SPH method becomes the SPH method if num = ∞.Figure 2 shows the positions of 1961 particles calculated by PDRI-SPH method with different interval time steps num = 10, 20, 40, 60, ∞ at the time 0.01 s.The distributed particles of using PDRI-SPH (num ̸ = ∞) are all more uniform and the outer surfaces are all far smoother than the ones of SPH method (num = ∞).From Figure 2, we can observe that the better results belong to num = 10 and 20.In fact, the more uniformly distributed the particles are, the better the numerical accuracy is (see [31]).In other words, the accuracy of numerical results using SPH can be improved by periodic density reinitialization with appropriate "num." For further exhibiting the merit of the PDRI-SPH method, the pressure field distribution (/10 7 ) and numerical accuracy obtained using PDRI-SPH are shown in Figures 3 and 4, respectively.We can know that the problem of pressure oscillations of using SPH can be effectively reduced by PDRI-SPH (see Figure 3).The pressure distribution has certain defect around the boundary region in Figure 3, due to the reduced particles on the boundary.Figure 4 demonstrates that the PDRI-SPH has better accuracy than the standard SPH.
In a word, the effect of the density reinitialization method used in the standard SPH method is obvious.Through the results of Figures 2-4, and considering the computational cost and the effect of PDRI-SPH with different "num, " we choose the interval time step num = 20 in all the following numerical simulations.

Validity of the Artificial Viscosity and Artificial Stress Models.
In this subsection, the example of a droplet impact and spreading on a horizontal or an inclined plate is considered.The initial state of a droplet impact on an inclined surface is shown in Figure 5(a).When a droplet impacts on the inclined rigid plate, the shape of the droplet distorts and spreads symmetrically ( = 0 ∘ ) or asymmetrically ( ̸ = 0 ∘ ) relative to the point of impact, which is shown in Figure 5(b).We define the positive value of the elongation  back and the negative value of the elongation  front (see Figure 5(b)), in which asymmetry increases with time.The front edge of the droplet spreads forward, while the back edge spreads backward or slips forward.All the physical parameter values for a Newtonian droplet are chosen as follows.Its initial diameter and velocity are  = 0.02 m and  0 = −1 ms −1 , respectively.The total viscosity is  = 4 Pa ⋅ s, the reference density is  0 = 10 3 kgm −3 , the speed of sound is  = 12 ms −1 , and the gravitational force acts downwards with   = −9.81ms −2 : in this subsection, 1961 fluid particles, 251 wall particles, and 753 ghost particles.The initial spacing  0 = 0.0004 m and the time-step  = 1 × 10 −5 .The height of dropping is  = 0.04 m from the center of drop to the center () of inclined rigid wall (see Figure 5(a)).For the Cross model droplet,  0 = 4 × 10 −3 m 2 s −1 ,  ∞ = 4 × 10 −4 m 2 s −1 , and the other parameter values are the same as the case of Newtonian droplet.
Although the artificial viscosity term (19a) is adopted for simulating the Newtonian drop case with considering a horizontal plate in [25], the role of the artificial viscosity has not been obviously demonstrated in [25].Here, the effect of artificial viscosity is illustrated in Figure 6.And then, we can obtain three advantages of using the artificial viscosity in the example of droplet impact.(1) The particles are more uniformly distributed than those without using it.(2) The numerical accuracy and stability are improved.(3) The phenomenon of unphysical clustering becomes weakened.Note that the artificial stress term and the density reinitialization method are all not considered in Figure 6.In all subsequent simulations, the artificial viscosity ( II = 1,  II = 2) is adopted.
Figures 7 and 8 show the effect of the artificial stress for simulating a Newtonian/Cross model droplet impact on horizontal/inclined ( = 30 ∘ ) rigid plate using PDRI-SPH method with different artificial stress parameters.It can be seen that the droplet fractures unrealistically for the problem of droplet impact without the artificial stress term ( = 0), and the simulations may be eventually diverged.The phenomenon of fracture is observable for the Newtonian droplet impact on horizontal rigid plate with artificial stress  = 0.2, but it is much severer when the Newtonian droplet impacts on inclined plate.For the Cross model droplet, the unphysical fracture is obvious no matter how the droplet impacts on horizontal or inclined rigid plate even if the artificial viscosity is adopted.Observing Figures 7 and 8, we can get the following.(1) The problem of tensile instability occurs more evidently for the Cross droplet than the Newtonian case when the droplet impacts on horizontal rigid plate.(2) A droplet impacts on rigid plate fracture more likely at low impact angles  ̸ = 0 ∘ than  = 0 ∘ .In fact, the tensile instability is also related to the ratio of the kinematic viscosity  0 ,  ∞ for a Cross droplet impact on plate.Here, we can find that the fracture is avoided completely by increasing the value of  up to 0.8 for the Newtonian/Cross model droplet impact on inclined plate.In other words, it is necessary that the artificial stress parameter  is chosen appropriately for simulating a droplet impact on inclined rigid plate at low impact angles using PDRI-SPH method.

Numerical Simulations Based on Cross Model Using PDRI-SPH.
In this subsection, we mainly focus on the PDRI-SPH/SPH method combined with the artificial viscosity  ( II = 1,  II = 2), artificial stress ( = 0.8), and boundary condition for simulating a Cross droplet impact on inclined rigid plate at different low impact angles  = 15 ∘ , 30 ∘ , 45 ∘ , respectively.For the purpose of comparison, the Newtonian droplet case is also considered.The number of fluid particles is set to 7845, corresponding to the initial spacing  0 = 0.0002 m, 501 wall particles, and 1503 ghost particles.The time-step is 5 × 10 −6 .The other physical parameters values are the same as those in Section 4.2.
The effect of the proposed periodic density reinitialization method is obviously shown by predicting the pressure distribution for the problem of droplet impact in Figure 9.At the short time of droplet impact, the phenomenon of pressure oscillations occurs for the SPH method combined with the above improved models.The pressure oscillations grow near the rigid plate varying with time and later progressively destroy the whole pressure field, resulting in making its physical interpretation and possible practical use difficult.However, the pressure filed maintains a much smoother character obtained using the proposed PDRI-SPH method than the SPH, especially on the boundary regions.We can also know that the wall particles and ghost particles contribute to the evolution of the density of the fluid particles; pressures on both fluid and virtual particles increase when fluid particles are near the rigid wall.The presented boundary treatment is strong enough to prevent fluid particles from penetrating the rigid wall without employing an additional artificial repulsive force.
Figure 10 shows the comparison of the numerical results obtained using SPH or PDRI-SPH method for the width of a Newtonian and Cross droplet varying with dimensionless time.The Newtonian/Cross droplet spreads symmetrically along the wall after impact with time.The PDRI-SPH results much closer to the results in [39] than the SPH results in the numerical simulations of Newtonian droplet impact.There are also certain differences between the numerical results of using PDRI-SPH and those of using SPH for solving the width of a Cross droplet impact on horizontal rigid plate with time in Figure 10(b).Considering the analysis of Section 4.1 and the results in Figure 9, it is not difficult to believe that the results obtained using the PDRI-SPH method are more reliable than those using SPH method.From Figure 10, we also can observe that the width of a Cross droplet becomes much larger than the corresponding Newtonian droplet case Mathematical Problems in Engineering  at the short time (about dimensionless time  ≤ 3) of droplet impact, due to the shear-thinning behavior of the Cross model droplet (see [29,30]).
In order to further demonstrate the feasibility and the credibility of the proposed method to simulate the impact problem, the numerical convergence of the PDRI-SPH results with different smoothing length ℎ and particles number   (along the -axis direction) is shown in Figure 11.From Figures 10 and 11, we can get that (a) the proposed method is convergent to simulate the Newtonian or Cross droplet impact on horizontal rigid plate under different smoothing length; (b) the results of Newtonian droplet for ℎ = 1.5 0 are more accurate than those for ℎ = 1.2 0 by observing Figures 10(a We can observe the shapes of a Newtonian drop spreading over inclined rigid surfaces ( = 15 ∘ , 30 ∘ , 45 ∘ ) at different dimensionless times from Figure 12.It can be seen that the first phase of impact involving the initial deformation of the droplet for all the cases of  = 15 ∘ , 30 ∘ , 45 ∘ is similar to that of impact angle  = 0 ∘ ; namely, the front edge spreads forward and the back edge spreads backward.Subsequently, the front edge spreads forward and the back edge slips forward, which can be clearly observed in Figure 14(a).The phenomenon of precipitation appears versus time, which becomes more evident with an increase of impact angle.The PDRI-SPH results are similar to the results of a water droplet impact on inclined surface obtained using VOF in [40,41].
From Figure 13, we can get that the shapes of the Cross drop spreading over inclined rigid surfaces ( = 15 ∘ , 30 ∘ , 45 ∘ ) at different dimensionless time are different from the case of Newtonian (see Figure 12) under the same total viscosity.There are some obvious differences between the Newtonian drop and the Cross model droplet case, which are shown in Figures 12-14.The speed of a Cross droplet spreading over inclined rigid plate is faster than its Newtonian counterpart.Moreover, an indentation is formed only for the Cross droplet with time at lower impact angle ( = 15 ∘ ) because of the shear-thinning of the Cross fluid.
To further exemplify the reliability of the proposed method for simulating the droplet spreading over inclined rigid surface, Figure 15 shows the numerical convergence of evolution of a Newtonian droplet and a Cross droplet spreading on an inclined plate with impact angle 30 ∘ and smoothing length ℎ = 1.5 0 .Observing Figure 15, the numerical results for   = 81 are very close to those for   = 101, which demonstrates that the proposed method possesses preferable numerical convergence for simulating the impact problem.In short, it is feasible and reliable to simulate the impact problem of a Newtonian or Cross droplet spreading on an inclined plate using the proposed PDRI-SPH method.

Conclusions
Aiming at the deficiency of standard SPH method, a periodic density reinitialization method which is called the PDRI-SPH method is proposed to preserve the consistency between mass, density, and the occupied area.In order to verify the validity and ability of the proposed PDRI-SPH, the benchmark problem of drop stretching is simulated by PDRI-SPH.Due to the density reinitialization, the PDRI-SPH has better accuracy than the SPH, and the distributed pressure field is much smoother likewise.Meanwhile, an artificial stress is successfully presented and tested to simulate a Cross droplet impact onto an inclined rigid plate of using PDRI-SPH.The effect of the proposed PDRI-SPH combined with an artificial viscosity, an artificial stress, and boundary condition treatment is further shown in the physical problem of impact of droplet onto inclined rigid plate and compared with the standard SPH.All the numerical results declare that the proposed PDRI-SPH has some merits comparing with the SPH, and it is a powerful tool to simulate the complex free surface for generalized Newtonian fluid.It is expected to be widely used and further improved to solve complex free surface flows in future.

Figure 2 :
Figure 2: The position of particles calculated by PDRI-SPH method at the time 0.01 s.

Figure 5 :
Figure 5: The initial state of a droplet impact (a) and side view of a droplet on an inclined rigid plate (b).

Figure 6 :
Figure 6: The particles distribution predicted by standard SPH method for a Newtonian droplet impact on horizontal rigid plate at dimensionless time  = 5.2.Without artificial viscosity (first row), artificial viscosity (second row).

8 Figure 7 :
Figure 7: The shape of a Newtonian droplet impact on horizontal (left) or inclined ( = 30 ∘ , right) rigid plate obtained by PDRI-SPH method with  II = 1,  II = 2 and different artificial stress parameter , at dimensionless time  = 6.2.

8 Figure 8 :
Figure 8: The shape of a Cross droplet impact on horizontal (left) or inclined ( = 30 ∘ , right) rigid plate obtained by PDRI-SPH method with  II = 1,  II = 2 and different artificial stress parameter , at dimensionless time  = 4.3.

Figure 9 :
Figure 9: The distribution of pressure contour for a Newtonian/Cross drop impact on horizontal rigid plate obtained using PDRI-SPH (left column) or SPH (right column) at different dimensionless time.

Figure 10 :
Figure 10: Comparison of the numerical results obtained using SPH or PDRI-SPH method for the width of a Newtonian (a) and Cross (b) droplet impact on horizontal rigid plate varying with dimensionless time.
) and 11(a), which implies that it is credible to adopt the ℎ = 1.5 0 in the simulations of Section 4.3; (c) the credibility of the PDRI-SPH for simulating the impact problem based on the Cross fluid is further verified by Figure 11(b).