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

Authors: Jinlian Ren; Weigang Lu; Tao Jiang
Journal: Mathematical Problems in Engineering
 (2015)

---

## Abstract

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-SPH method. 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.


---

## Body

## 1. 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 shear-thinning 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–8] or particle methods have been proposed in a Lagrangian framework. Among the various particle methods [9–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–25]. In this work, a periodic density reinitialization method based on 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 shear-thinning 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 Section2. 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.


## 2. 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:(1)
D
ρ
D
t
=
-
ρ
∂
v
β
∂
x
β
,
(2)
D
v
α
D
t
=
1
ρ
∂
σ
α
β
∂
x
β
+
g
α
,where ρ denotes the fluid density, v
β the βth component of the fluid velocity, σ
α
β the (
α
,
β
)th component of the total stress tensor, g
α is the αth component of the gravitational acceleration, and the D
/
D
t
=
∂
/
∂
t
+
v
β
·
∂
/
∂
x
β is the material derivative. The spatial coordinates x
β and time t are the independent variables.The total stress tensor in (2) is commonly made up of the isotropic pressure p and the components of extra stress tensor τ
α
β:(3)
σ
α
β
=
-
p
δ
α
β
+
τ
α
β
,where δ
α
β
=
1 if α
=
β and δ
α
β
=
0 if α
≠
β. In order to study a shear-thinning polymer material, the relating constitutive equation must be provided.


### 2.2. 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 (4)
τ
=
2
μ
γ
˙
d
,where μ
(
γ
˙
)
=
ρ
υ
(
γ
˙
) is the dynamic viscosity and the υ
(
γ
˙
) is the kinematic viscosity. The symmetric strain rate tensor d is defined as(5)
d
=
d
α
β
=
1
2
∂
v
α
∂
x
β
+
∂
v
β
∂
x
α
.The kinematic viscosity υ
(
γ
˙
) represents the shear-thinning nature of the fluid; it is defined as(6)
υ
γ
˙
=
υ
∞
+
υ
0
-
υ
∞
1
+
K
γ
˙
m
,where m, υ
0, υ
∞, and K are given positive constants. The γ
˙ is the local shear rate defined by (7)
γ
˙
=
2
t
r
d
2
1
/
2
,and the symbol “t
r” denotes the trace of matrix. Here, the positive constants K and m are all chosen equal to 1.0.


### 2.3. Equation of State

The incompressible flows were sometimes treated as slightly compressible flows by adopting a suitable equation of state in many previous works (see 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]:(8)
p
ρ
=
c
2
ρ
2
2
ρ
0
,where c 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 M (M
≡
V
/
c, where V 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.


## 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:(1)
D
ρ
D
t
=
-
ρ
∂
v
β
∂
x
β
,
(2)
D
v
α
D
t
=
1
ρ
∂
σ
α
β
∂
x
β
+
g
α
,where ρ denotes the fluid density, v
β the βth component of the fluid velocity, σ
α
β the (
α
,
β
)th component of the total stress tensor, g
α is the αth component of the gravitational acceleration, and the D
/
D
t
=
∂
/
∂
t
+
v
β
·
∂
/
∂
x
β is the material derivative. The spatial coordinates x
β and time t are the independent variables.The total stress tensor in (2) is commonly made up of the isotropic pressure p and the components of extra stress tensor τ
α
β:(3)
σ
α
β
=
-
p
δ
α
β
+
τ
α
β
,where δ
α
β
=
1 if α
=
β and δ
α
β
=
0 if α
≠
β. In order to study a shear-thinning polymer material, the relating constitutive equation must be provided.


## 2.2. 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 (4)
τ
=
2
μ
γ
˙
d
,where μ
(
γ
˙
)
=
ρ
υ
(
γ
˙
) is the dynamic viscosity and the υ
(
γ
˙
) is the kinematic viscosity. The symmetric strain rate tensor d is defined as(5)
d
=
d
α
β
=
1
2
∂
v
α
∂
x
β
+
∂
v
β
∂
x
α
.The kinematic viscosity υ
(
γ
˙
) represents the shear-thinning nature of the fluid; it is defined as(6)
υ
γ
˙
=
υ
∞
+
υ
0
-
υ
∞
1
+
K
γ
˙
m
,where m, υ
0, υ
∞, and K are given positive constants. The γ
˙ is the local shear rate defined by (7)
γ
˙
=
2
t
r
d
2
1
/
2
,and the symbol “t
r” denotes the trace of matrix. Here, the positive constants K and m are all chosen equal to 1.0.


## 2.3. Equation of State

The incompressible flows were sometimes treated as slightly compressible flows by adopting a suitable equation of state in many previous works (see 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]:(8)
p
ρ
=
c
2
ρ
2
2
ρ
0
,where c 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 M (M
≡
V
/
c, where V 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.


## 3. 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 m, velocity v, and other physical quantities depending on the problem. Any function f
(
r
) defined at the position r
=
x
,
y can be expressed by the following integral:(9)
f
(
r
)
=
∫
Ω
f
r
′
W
r
-
r
′
,
h
d
r
′
,where W represents the kernel function (or smoothing function) and h denotes the smoothing length defining the influence area of W. The kernel function W 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 f
(
r
), and its first derivative at the position r
=
x
,
y of the particle i can be written in the following condensed forms:(10)
f
i
=
∑
j
m
j
ρ
j
f
j
W
i
j
,
(11)
∂
f
∂
r
i
=
∑
j
m
j
ρ
j
(
f
j
-
f
i
)
∂
W
i
j
∂
r
i
,where m
j and ρ
j are the mass and density of the jth particle, and f
j
=
f
r
j. The m
j
/
ρ
j represents the occupied volume by the jth particle. The W
i
j
=
W
r
i
-
r
j
,
h, ∂
W
i
j
/
∂
r
i
=
-
∂
W
j
i
/
∂
r
j.In this paper, the cubic spline function is chosen as the smoothing function which is the function aboutr
i
j
=
r
i
-
r
j and q
=
r
i
j
/
h. Then it reads for 2-D as follows:(12)
W
i
j
=
W
(
r
,
h
)
=
15
7
π
h
2
2
3
-
q
2
+
1
2
q
3
,
0
≤
q
<
1
,
1
6
2
-
q
3
,
1
≤
q
<
2
,
0
,
q
≥
2
.In order to have an accurate interpolation, the smoothing length h should be chosen bigger than the mean interparticle distance. Here, the smoothing length h is given by 1.5
d
0 with d
0 as the initial distance between neighboring particles. The compact support domain size is 2
h.Considering the discrete gradient equation (11) and the following identity: (
1
/
ρ
)
(
∂
σ
α
β
/
∂
x
β
)
=
∂
σ
α
β
/
ρ
/
∂
x
β
+
(
σ
α
β
/
ρ
2
)
(
∂
ρ
/
∂
x
β
), the particle discretization schemes of the governing equations can be obtained at the particle i: (13)
D
ρ
D
t
i
=
ρ
i
∑
j
m
j
ρ
j
(
v
i
β
-
v
j
β
)
∂
W
i
j
∂
x
i
β
,
(14)
D
v
α
D
t
i
=
∑
j
m
j
σ
i
α
β
ρ
i
2
+
σ
j
α
β
ρ
j
2
∂
W
i
j
∂
x
i
β
+
g
α
.Introducing the velocity gradient(15)
κ
i
α
β
=
∂
v
α
∂
x
β
i
=
∑
j
m
j
ρ
j
(
v
j
α
-
v
i
α
)
∂
W
i
j
∂
x
i
β
,then the particle approximation schemes of constitutive equation for the Cross model can be obtained as (16)
σ
i
α
β
=
-
p
i
δ
α
β
+
μ
γ
˙
i
κ
i
α
β
+
κ
i
β
α
,
γ
˙
i
=
2
t
r
d
i
2
1
/
2
.


### 3.2. 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]:(17)
ρ
i
=
∑
j
m
j
W
i
j
,the particle approximations (17) do not have first order accuracy and C
0, C
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 second-order accurate particle approximation scheme based on Taylor series expansion (see [27, 31]) to periodically reinitialize the density field:(18a)
ρ
i
=
∑
j
m
j
W
i
j
T
a
y
,where the corrected kernel function W
i
j
T
a
y is given by(18b)
W
i
j
T
a
y
=
B
∑
j
A
r
i
V
j
-
1
W
,
(18c)
A
r
i
=
W
i
j
x
j
i
·
W
i
j
y
j
i
·
W
i
j
∂
W
i
j
∂
x
i
x
j
i
·
∂
W
i
j
∂
x
i
y
j
i
·
∂
W
i
j
∂
x
i
∂
W
i
j
∂
y
i
x
j
i
·
∂
W
i
j
∂
y
i
y
j
i
·
∂
W
i
j
∂
y
i
,
B
=
1,0
,
0
,
W
=
W
i
j
∂
W
i
j
∂
x
i
∂
W
i
j
∂
y
i
.Here x
j
i
=
x
j
-
x
i, y
j
i
=
y
j
-
y
i, and V
j is replaced by m
j
/
ρ
j. The particle approximations scheme (18a)–(18c) possesses C
0, C
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 of3
×
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 C
0 consistency on the whole domain, the matrix ∑
j
A
r
i
V
j may be replaced by ∑
j
W
i
j
V
j 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.


### 3.3. 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] (19a)
Π
i
j
=
-
α
Π
c
-
i
j
ϕ
i
j
+
β
Π
ϕ
i
j
2
ρ
-
i
j
,
v
i
j
·
r
i
j
<
0
,
0
,
v
i
j
·
r
i
j
≥
0
,where (19b)
ϕ
i
j
=
h
u
i
j
·
r
i
j
r
i
j
2
+
0.01
h
2
,
c
-
i
j
=
c
i
+
c
j
2
,
ρ
-
i
j
=
ρ
i
+
ρ
j
2
,
v
i
j
=
v
i
-
v
j
,
r
i
j
=
r
i
-
r
j
.The 0.01
h
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.


### 3.4. 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”: (20a)
f
i
j
n
S
i
α
β
+
S
j
α
β
,where n
=
W
0
,
h
/
W
d
0
,
h, f
i
j
=
W
i
j
r
i
-
r
j
,
h
/
W
i
j
d
0
,
h. The components of the artificial stress tensor S
i
α
β are given as (20b)
S
i
α
β
=
-
b
σ
i
α
β
ρ
2
,
i
f
∑
j
W
i
j
′
′
α
β
·
σ
i
α
β
>
0
,
0
,
o
t
h
e
r
,where b is a positive parameter (
0
<
b
<
1
), and the W
i
j
′
′
α
β
=
∂
2
W
i
j
/
∂
x
α
·
∂
x
β.Introducing theK
i
j
=
f
i
j
n
S
i
α
β
+
S
j
α
β, 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 (21)
D
v
α
D
t
i
=
∑
j
m
j
σ
i
α
β
ρ
i
2
+
σ
j
α
β
ρ
j
2
-
Π
i
j
δ
α
β
+
K
i
j
∂
W
i
j
∂
x
i
β
+
g
α
.Usually, the particle positions are updated by the following equation:(22)
D
x
i
α
D
t
=
v
i
α
.


### 3.5. 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 Figure1, 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:(23a)
p
i
=
∑
j
p
j
2
h
-
r
i
j
∑
j
2
h
-
r
i
j
,where i represents the index of a wall particles and j denotes the index of its neighboring fluid particles only.Figure 1
The sketch of wall particle and ghost particle at impact angle ofα.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 particleG, we assign a corresponding point W just on the rigid wall and point F 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 distancesd
G and d
F of the points G and F to the rigid wall, respectively.
(3)
The pressurep
G and velocity v
I
α for the ghost particles can be obtained through the following linear extrapolation: (23b)
L
G
=
L
F
+
1
+
d
G
d
F
L
W
-
L
F
,whereL represents the vector of variables p
,
v
α. To specify the values for L
F, the interpolation formulation (23a) is applied again. Here, we let the d
F
=
d
0.Moreover, the following total stress-free condition must be satisfied in the computational domain for surface particles:(24)
σ
·
n
=
0
,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 (25a)
X
~
i
n
+
1
=
X
i
n
+
Δ
t
2
Γ
i
n
,where X
i represents the vector of the unknown variables ρ
i
,
v
i
α and Γ
i denotes the vector of right-hand sides of (13), (21), and (22). In the corrected step, the updated value of X
i at the end of each time step is given by(25b)
X
i
n
+
1
=
X
i
n
+
Δ
t
2
Γ
i
n
+
Γ
~
i
n
.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:(25c)
Δ
t
≤
m
i
n
0.25
h
c
,
0.25
h
F
a
1
/
2
,
0.125
h
2
υ
0
,where F
a is the hydrodynamical force acting on the particle, and υ
0
=
η
/
ρ
0 is the kinematic viscosity.


## 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 m, velocity v, and other physical quantities depending on the problem. Any function f
(
r
) defined at the position r
=
x
,
y can be expressed by the following integral:(9)
f
(
r
)
=
∫
Ω
f
r
′
W
r
-
r
′
,
h
d
r
′
,where W represents the kernel function (or smoothing function) and h denotes the smoothing length defining the influence area of W. The kernel function W 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 f
(
r
), and its first derivative at the position r
=
x
,
y of the particle i can be written in the following condensed forms:(10)
f
i
=
∑
j
m
j
ρ
j
f
j
W
i
j
,
(11)
∂
f
∂
r
i
=
∑
j
m
j
ρ
j
(
f
j
-
f
i
)
∂
W
i
j
∂
r
i
,where m
j and ρ
j are the mass and density of the jth particle, and f
j
=
f
r
j. The m
j
/
ρ
j represents the occupied volume by the jth particle. The W
i
j
=
W
r
i
-
r
j
,
h, ∂
W
i
j
/
∂
r
i
=
-
∂
W
j
i
/
∂
r
j.In this paper, the cubic spline function is chosen as the smoothing function which is the function aboutr
i
j
=
r
i
-
r
j and q
=
r
i
j
/
h. Then it reads for 2-D as follows:(12)
W
i
j
=
W
(
r
,
h
)
=
15
7
π
h
2
2
3
-
q
2
+
1
2
q
3
,
0
≤
q
<
1
,
1
6
2
-
q
3
,
1
≤
q
<
2
,
0
,
q
≥
2
.In order to have an accurate interpolation, the smoothing length h should be chosen bigger than the mean interparticle distance. Here, the smoothing length h is given by 1.5
d
0 with d
0 as the initial distance between neighboring particles. The compact support domain size is 2
h.Considering the discrete gradient equation (11) and the following identity: (
1
/
ρ
)
(
∂
σ
α
β
/
∂
x
β
)
=
∂
σ
α
β
/
ρ
/
∂
x
β
+
(
σ
α
β
/
ρ
2
)
(
∂
ρ
/
∂
x
β
), the particle discretization schemes of the governing equations can be obtained at the particle i: (13)
D
ρ
D
t
i
=
ρ
i
∑
j
m
j
ρ
j
(
v
i
β
-
v
j
β
)
∂
W
i
j
∂
x
i
β
,
(14)
D
v
α
D
t
i
=
∑
j
m
j
σ
i
α
β
ρ
i
2
+
σ
j
α
β
ρ
j
2
∂
W
i
j
∂
x
i
β
+
g
α
.Introducing the velocity gradient(15)
κ
i
α
β
=
∂
v
α
∂
x
β
i
=
∑
j
m
j
ρ
j
(
v
j
α
-
v
i
α
)
∂
W
i
j
∂
x
i
β
,then the particle approximation schemes of constitutive equation for the Cross model can be obtained as (16)
σ
i
α
β
=
-
p
i
δ
α
β
+
μ
γ
˙
i
κ
i
α
β
+
κ
i
β
α
,
γ
˙
i
=
2
t
r
d
i
2
1
/
2
.


## 3.2. 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]:(17)
ρ
i
=
∑
j
m
j
W
i
j
,the particle approximations (17) do not have first order accuracy and C
0, C
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 second-order accurate particle approximation scheme based on Taylor series expansion (see [27, 31]) to periodically reinitialize the density field:(18a)
ρ
i
=
∑
j
m
j
W
i
j
T
a
y
,where the corrected kernel function W
i
j
T
a
y is given by(18b)
W
i
j
T
a
y
=
B
∑
j
A
r
i
V
j
-
1
W
,
(18c)
A
r
i
=
W
i
j
x
j
i
·
W
i
j
y
j
i
·
W
i
j
∂
W
i
j
∂
x
i
x
j
i
·
∂
W
i
j
∂
x
i
y
j
i
·
∂
W
i
j
∂
x
i
∂
W
i
j
∂
y
i
x
j
i
·
∂
W
i
j
∂
y
i
y
j
i
·
∂
W
i
j
∂
y
i
,
B
=
1,0
,
0
,
W
=
W
i
j
∂
W
i
j
∂
x
i
∂
W
i
j
∂
y
i
.Here x
j
i
=
x
j
-
x
i, y
j
i
=
y
j
-
y
i, and V
j is replaced by m
j
/
ρ
j. The particle approximations scheme (18a)–(18c) possesses C
0, C
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 of3
×
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 C
0 consistency on the whole domain, the matrix ∑
j
A
r
i
V
j may be replaced by ∑
j
W
i
j
V
j 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.


## 3.3. 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] (19a)
Π
i
j
=
-
α
Π
c
-
i
j
ϕ
i
j
+
β
Π
ϕ
i
j
2
ρ
-
i
j
,
v
i
j
·
r
i
j
<
0
,
0
,
v
i
j
·
r
i
j
≥
0
,where (19b)
ϕ
i
j
=
h
u
i
j
·
r
i
j
r
i
j
2
+
0.01
h
2
,
c
-
i
j
=
c
i
+
c
j
2
,
ρ
-
i
j
=
ρ
i
+
ρ
j
2
,
v
i
j
=
v
i
-
v
j
,
r
i
j
=
r
i
-
r
j
.The 0.01
h
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.


## 3.4. 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”: (20a)
f
i
j
n
S
i
α
β
+
S
j
α
β
,where n
=
W
0
,
h
/
W
d
0
,
h, f
i
j
=
W
i
j
r
i
-
r
j
,
h
/
W
i
j
d
0
,
h. The components of the artificial stress tensor S
i
α
β are given as (20b)
S
i
α
β
=
-
b
σ
i
α
β
ρ
2
,
i
f
∑
j
W
i
j
′
′
α
β
·
σ
i
α
β
>
0
,
0
,
o
t
h
e
r
,where b is a positive parameter (
0
<
b
<
1
), and the W
i
j
′
′
α
β
=
∂
2
W
i
j
/
∂
x
α
·
∂
x
β.Introducing theK
i
j
=
f
i
j
n
S
i
α
β
+
S
j
α
β, 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 (21)
D
v
α
D
t
i
=
∑
j
m
j
σ
i
α
β
ρ
i
2
+
σ
j
α
β
ρ
j
2
-
Π
i
j
δ
α
β
+
K
i
j
∂
W
i
j
∂
x
i
β
+
g
α
.Usually, the particle positions are updated by the following equation:(22)
D
x
i
α
D
t
=
v
i
α
.


## 3.5. 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 Figure1, 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:(23a)
p
i
=
∑
j
p
j
2
h
-
r
i
j
∑
j
2
h
-
r
i
j
,where i represents the index of a wall particles and j denotes the index of its neighboring fluid particles only.Figure 1
The sketch of wall particle and ghost particle at impact angle ofα.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 particleG, we assign a corresponding point W just on the rigid wall and point F 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 distancesd
G and d
F of the points G and F to the rigid wall, respectively.
(3)
The pressurep
G and velocity v
I
α for the ghost particles can be obtained through the following linear extrapolation: (23b)
L
G
=
L
F
+
1
+
d
G
d
F
L
W
-
L
F
,whereL represents the vector of variables p
,
v
α. To specify the values for L
F, the interpolation formulation (23a) is applied again. Here, we let the d
F
=
d
0.Moreover, the following total stress-free condition must be satisfied in the computational domain for surface particles:(24)
σ
·
n
=
0
,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 (25a)
X
~
i
n
+
1
=
X
i
n
+
Δ
t
2
Γ
i
n
,where X
i represents the vector of the unknown variables ρ
i
,
v
i
α and Γ
i denotes the vector of right-hand sides of (13), (21), and (22). In the corrected step, the updated value of X
i at the end of each time step is given by(25b)
X
i
n
+
1
=
X
i
n
+
Δ
t
2
Γ
i
n
+
Γ
~
i
n
.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:(25c)
Δ
t
≤
m
i
n
0.25
h
c
,
0.25
h
F
a
1
/
2
,
0.125
h
2
υ
0
,where F
a is the hydrodynamical force acting on the particle, and υ
0
=
η
/
ρ
0 is the kinematic viscosity.


## 4. 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 Section4.3.


### 4.1. 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
k
g
m
-
3, the viscosity η
=
0.001
k
g
m
-
1
s
-
1, and the speed of sound c
=
1400
m
s
-
1. The initial geometry of the water drop is a circle of radius R
=
1
m with its center located at the origin (x
=
0, y
=
0). There is no external forces but initial velocity field v
0
x
=
-
A
0
x, v
0
x
=
A
0
y with A
0
=
100
s
-
1 and the initial pressure field p
0
=
1
/
2
ρ
0
A
0
2
R
2
-
x
2
+
y
2. The number of fluid particles is 1961 and corresponding to the initial distance d
0
=
0.02
m, and the time step d
t
=
1
0
-
5
s. During the stretching process of water drop, the water drop remains elliptical shape and the value of a
~
·
b
~ (a
~ is the semiminor axis and b
~ is the semimajor axis) remains constant. We let “n
u
m” denote the interval time step of the periodic density reinitialization of PDRI-SPH, and the PDRI-SPH method becomes the SPH method if n
u
m
=
∞.Figure2 shows the positions of 1961 particles calculated by PDRI-SPH method with different interval time steps n
u
m
=
10,20,40,60
,
∞ at the time 0.01 s. The distributed particles of using PDRI-SPH (n
u
m
≠
∞) are all more uniform and the outer surfaces are all far smoother than the ones of SPH method (n
u
m
=
∞). From Figure 2, we can observe that the better results belong to n
u
m
=
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 (p
/
1
0
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.Figure 2
The position of particles calculated by PDRI-SPH method at the time 0.01 s.Figure 3
The pressure field distribution (p
/
1
0
7) obtained using the PDRI-SPH method with different “num” at time 0.003 s and 0.006 s, respectively.Figure 4
The comparisons of numerical results obtained using PDRI-SPH method about the semiminor axis varying with time.In a word, the effect of the density reinitialization method used in the standard SPH method is obvious. Through the results of Figures2–4, and considering the computational cost and the effect of PDRI-SPH with different “num,” we choose the interval time step n
u
m
=
20 in all the following numerical simulations.


### 4.2. 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 Figure5(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 X
b
a
c
k and the negative value of the elongation X
f
r
o
n
t (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.Figure 5
The initial state of a droplet impact (a) and side view of a droplet on an inclined rigid plate (b).
(a)
(b)All the physical parameter values for a Newtonian droplet are chosen as follows. Its initial diameter and velocity areD
=
0.02
m and U
0
=
-
1
m
s
-
1, respectively. The total viscosity is η
=
4
P
a
·
s, the reference density is ρ
0
=
1
0
3
k
g
m
-
3, the speed of sound is c
=
12
m
s
-
1, and the gravitational force acts downwards with g
y
=
-
9.81
m
s
-
2: in this subsection, 1961 fluid particles, 251 wall particles, and 753 ghost particles. The initial spacing d
0
=
0.0004
m and the time-step d
t
=
1
×
1
0
-
5. The height of dropping is H
=
0.04
m from the center of drop to the center (o) of inclined rigid wall (see Figure 5(a)). For the Cross model droplet, υ
0
=
4
×
1
0
-
3
m
2
s
-
1, υ
∞
=
4
×
1
0
-
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 (α
I
I
=
1, β
I
I
=
2) is adopted.Figure 6
The particles distribution predicted by standard SPH method for a Newtonian droplet impact on horizontal rigid plate at dimensionless timet
=
5.2. Without artificial viscosity (first row), artificial viscosity (second row).Figures7 and 8 show the effect of the artificial stress for simulating a Newtonian/Cross model droplet impact on horizontal/inclined (α
=
3
0
°) 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 (b
=
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 b
=
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 b 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 b is chosen appropriately for simulating a droplet impact on inclined rigid plate at low impact angles using PDRI-SPH method.Figure 7
The shape of a Newtonian droplet impact on horizontal (left) or inclined (α
=
3
0
°, right) rigid plate obtained by PDRI-SPH method with α
I
I
=
1, β
I
I
=
2 and different artificial stress parameter b, at dimensionless time t
=
6.2.Figure 8
The shape of a Cross droplet impact on horizontal (left) or inclined (α
=
3
0
°, right) rigid plate obtained by PDRI-SPH method with α
I
I
=
1, β
I
I
=
2 and different artificial stress parameter b, at dimensionless time t
=
4.3.


### 4.3. 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 (α
I
I
=
1, β
I
I
=
2), artificial stress (b
=
0.8), and boundary condition for simulating a Cross droplet impact on inclined rigid plate at different low impact angles α
=
1
5
°
,
3
0
°
,
4
5
°, 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 d
0
=
0.0002
m, 501 wall particles, and 1503 ghost particles. The time-step is 5
×
1
0
-
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 Figure9. 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 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.Figure10 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 at the short time (about dimensionless time t
≤
3) of droplet impact, due to the shear-thinning behavior of the Cross model droplet (see [29, 30]).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.
(a)
(b)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 lengthh and particles number N
x (along the x-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 h
=
1.5
d
0 are more accurate than those for h
=
1.2
d
0 by observing  Figures 10(a) and 11(a),  which implies that it is credible to adopt the h
=
1.5
d
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).Figure 11
Numerical convergence of the PDRI-SPH results with different smoothing lengthh and particles number N
x (along the x-axis direction): (a1)-(a2) Newtonian droplet; (b1)-(b2) Cross droplet (which corresponds to Figure 10).
(a)
(b)We can observe the shapes of a Newtonian drop spreading over inclined rigid surfaces (α
=
1
5
°
,
3
0
°
,
4
5
°) 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 α
=
1
5
°
,
3
0
°
,
4
5
° 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].Figure 12
The shape of a Newtonian droplet at impact angles of1
5
° (first column), 3
0
° (second column), and 4
5
° (third column) obtained using PDRI-SPH, varying with short dimensionless time.From Figure13, we can get that the shapes of the Cross drop spreading over inclined rigid surfaces (α
=
1
5
°
,
3
0
°
,
4
5
°) 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 (α
=
1
5
°) because of the shear-thinning of the Cross fluid.Figure 13
The shape of a Cross droplet at impact angles of1
5
° (first column), 3
0
° (second column), and 4
5
° (third column) obtained using PDRI-SPH, varying with short dimensionless time.Figure 14
Evolution of the elongation of a Newtonian (a) and a Cross (b) spreading droplet at different impact angles, obtained using PDRI-SPH.
(a)
(b)To further exemplify the reliability of the proposed method for simulating the droplet spreading over inclined rigid surface, Figure15 shows the numerical convergence of evolution of a Newtonian droplet and a Cross droplet spreading on an inclined plate with impact angle 3
0
° and smoothing length h
=
1.5
d
0. Observing Figure 15, the numerical results for N
x
=
81 are very close to those for N
x
=
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.Figure 15
Numerical convergence of evolution of a Newtonian droplet (a) and a Cross droplet (b) spreading on an inclined plate with impact angle3
0
° and smoothing length h
=
1.5
d
0 obtained using PDRI-SPH method.
(a)
(b)


## 4.1. 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
k
g
m
-
3, the viscosity η
=
0.001
k
g
m
-
1
s
-
1, and the speed of sound c
=
1400
m
s
-
1. The initial geometry of the water drop is a circle of radius R
=
1
m with its center located at the origin (x
=
0, y
=
0). There is no external forces but initial velocity field v
0
x
=
-
A
0
x, v
0
x
=
A
0
y with A
0
=
100
s
-
1 and the initial pressure field p
0
=
1
/
2
ρ
0
A
0
2
R
2
-
x
2
+
y
2. The number of fluid particles is 1961 and corresponding to the initial distance d
0
=
0.02
m, and the time step d
t
=
1
0
-
5
s. During the stretching process of water drop, the water drop remains elliptical shape and the value of a
~
·
b
~ (a
~ is the semiminor axis and b
~ is the semimajor axis) remains constant. We let “n
u
m” denote the interval time step of the periodic density reinitialization of PDRI-SPH, and the PDRI-SPH method becomes the SPH method if n
u
m
=
∞.Figure2 shows the positions of 1961 particles calculated by PDRI-SPH method with different interval time steps n
u
m
=
10,20,40,60
,
∞ at the time 0.01 s. The distributed particles of using PDRI-SPH (n
u
m
≠
∞) are all more uniform and the outer surfaces are all far smoother than the ones of SPH method (n
u
m
=
∞). From Figure 2, we can observe that the better results belong to n
u
m
=
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 (p
/
1
0
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.Figure 2
The position of particles calculated by PDRI-SPH method at the time 0.01 s.Figure 3
The pressure field distribution (p
/
1
0
7) obtained using the PDRI-SPH method with different “num” at time 0.003 s and 0.006 s, respectively.Figure 4
The comparisons of numerical results obtained using PDRI-SPH method about the semiminor axis varying with time.In a word, the effect of the density reinitialization method used in the standard SPH method is obvious. Through the results of Figures2–4, and considering the computational cost and the effect of PDRI-SPH with different “num,” we choose the interval time step n
u
m
=
20 in all the following numerical simulations.


## 4.2. 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 Figure5(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 X
b
a
c
k and the negative value of the elongation X
f
r
o
n
t (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.Figure 5
The initial state of a droplet impact (a) and side view of a droplet on an inclined rigid plate (b).
(a)
(b)All the physical parameter values for a Newtonian droplet are chosen as follows. Its initial diameter and velocity areD
=
0.02
m and U
0
=
-
1
m
s
-
1, respectively. The total viscosity is η
=
4
P
a
·
s, the reference density is ρ
0
=
1
0
3
k
g
m
-
3, the speed of sound is c
=
12
m
s
-
1, and the gravitational force acts downwards with g
y
=
-
9.81
m
s
-
2: in this subsection, 1961 fluid particles, 251 wall particles, and 753 ghost particles. The initial spacing d
0
=
0.0004
m and the time-step d
t
=
1
×
1
0
-
5. The height of dropping is H
=
0.04
m from the center of drop to the center (o) of inclined rigid wall (see Figure 5(a)). For the Cross model droplet, υ
0
=
4
×
1
0
-
3
m
2
s
-
1, υ
∞
=
4
×
1
0
-
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 (α
I
I
=
1, β
I
I
=
2) is adopted.Figure 6
The particles distribution predicted by standard SPH method for a Newtonian droplet impact on horizontal rigid plate at dimensionless timet
=
5.2. Without artificial viscosity (first row), artificial viscosity (second row).Figures7 and 8 show the effect of the artificial stress for simulating a Newtonian/Cross model droplet impact on horizontal/inclined (α
=
3
0
°) 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 (b
=
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 b
=
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 b 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 b is chosen appropriately for simulating a droplet impact on inclined rigid plate at low impact angles using PDRI-SPH method.Figure 7
The shape of a Newtonian droplet impact on horizontal (left) or inclined (α
=
3
0
°, right) rigid plate obtained by PDRI-SPH method with α
I
I
=
1, β
I
I
=
2 and different artificial stress parameter b, at dimensionless time t
=
6.2.Figure 8
The shape of a Cross droplet impact on horizontal (left) or inclined (α
=
3
0
°, right) rigid plate obtained by PDRI-SPH method with α
I
I
=
1, β
I
I
=
2 and different artificial stress parameter b, at dimensionless time t
=
4.3.


## 4.3. 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 (α
I
I
=
1, β
I
I
=
2), artificial stress (b
=
0.8), and boundary condition for simulating a Cross droplet impact on inclined rigid plate at different low impact angles α
=
1
5
°
,
3
0
°
,
4
5
°, 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 d
0
=
0.0002
m, 501 wall particles, and 1503 ghost particles. The time-step is 5
×
1
0
-
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 Figure9. 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 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.Figure10 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 at the short time (about dimensionless time t
≤
3) of droplet impact, due to the shear-thinning behavior of the Cross model droplet (see [29, 30]).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.
(a)
(b)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 lengthh and particles number N
x (along the x-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 h
=
1.5
d
0 are more accurate than those for h
=
1.2
d
0 by observing  Figures 10(a) and 11(a),  which implies that it is credible to adopt the h
=
1.5
d
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).Figure 11
Numerical convergence of the PDRI-SPH results with different smoothing lengthh and particles number N
x (along the x-axis direction): (a1)-(a2) Newtonian droplet; (b1)-(b2) Cross droplet (which corresponds to Figure 10).
(a)
(b)We can observe the shapes of a Newtonian drop spreading over inclined rigid surfaces (α
=
1
5
°
,
3
0
°
,
4
5
°) 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 α
=
1
5
°
,
3
0
°
,
4
5
° 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].Figure 12
The shape of a Newtonian droplet at impact angles of1
5
° (first column), 3
0
° (second column), and 4
5
° (third column) obtained using PDRI-SPH, varying with short dimensionless time.From Figure13, we can get that the shapes of the Cross drop spreading over inclined rigid surfaces (α
=
1
5
°
,
3
0
°
,
4
5
°) 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 (α
=
1
5
°) because of the shear-thinning of the Cross fluid.Figure 13
The shape of a Cross droplet at impact angles of1
5
° (first column), 3
0
° (second column), and 4
5
° (third column) obtained using PDRI-SPH, varying with short dimensionless time.Figure 14
Evolution of the elongation of a Newtonian (a) and a Cross (b) spreading droplet at different impact angles, obtained using PDRI-SPH.
(a)
(b)To further exemplify the reliability of the proposed method for simulating the droplet spreading over inclined rigid surface, Figure15 shows the numerical convergence of evolution of a Newtonian droplet and a Cross droplet spreading on an inclined plate with impact angle 3
0
° and smoothing length h
=
1.5
d
0. Observing Figure 15, the numerical results for N
x
=
81 are very close to those for N
x
=
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.Figure 15
Numerical convergence of evolution of a Newtonian droplet (a) and a Cross droplet (b) spreading on an inclined plate with impact angle3
0
° and smoothing length h
=
1.5
d
0 obtained using PDRI-SPH method.
(a)
(b)


## 5. 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.



---

*Source: 915973-2015-06-24.xml*