Table of Contents Author Guidelines Submit a Manuscript
Mathematical Problems in Engineering
Volume 2015, Article ID 915973, 17 pages
http://dx.doi.org/10.1155/2015/915973
Research Article

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

1Department of Mathematics, School of Hydraulic, Energy and Power Engineering, Yangzhou University, Yangzhou 225002, China
2Department of Applied Mathematics, Northwestern Polytechnical University, Xi’an 710129, China

Received 15 June 2014; Revised 22 October 2014; Accepted 23 October 2014

Academic Editor: Kim M. Liew

Copyright © 2015 Jinlian Ren et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

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.

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 [68] or particle methods have been proposed in a Lagrangian framework. Among the various particle methods [912], 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. It handles convection dominated flows and large deformation problems without any numerical diffusion. Complex free surfaces are modeled easily and naturally without the need of explicit surface tracking technique. 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 [2325]. 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 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.

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:where denotes the fluid density, 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 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 if and 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 where is the dynamic viscosity and the is the kinematic viscosity. The symmetric strain rate tensor is defined asThe kinematic viscosity represents the shear-thinning nature of the fluid; it is defined aswhere , , , and are given positive constants. The is the local shear rate defined by and the symbol “” denotes the trace of matrix. Here, the positive constants and are all chosen equal to .

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]:where is the speed of sound and 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.

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 , velocity , and other physical quantities depending on the problem. Any function defined at the position 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 , and its first derivative at the position of the particle can be written in the following condensed forms:where and are the mass and density of the th particle, and . The represents the occupied volume by the th particle. The , .

In this paper, the cubic spline function is chosen as the smoothing function which is the function about 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 with as the initial distance between neighboring particles. The compact support domain size is .

Considering the discrete gradient equation (11) and the following identity: , 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

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]:the particle approximations (17) do not have first order accuracy and , 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:where the corrected kernel function is given byHere , , and is replaced by . The particle approximations scheme (18a)–(18c) possesses , 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 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 consistency on the whole domain, the matrix 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.

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] where The 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”: where , . The components of the artificial stress tensor are given as where is a positive parameter , 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:

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, 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. 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.

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 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 for the ghost particles can be obtained through the following linear extrapolation:

where represents the vector of variables . To specify the values for , the interpolation formulation (23a) is applied again. Here, we let the .

Moreover, the following total stress-free condition must be satisfied in the computational domain for surface particles: where 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 represents the vector of the unknown variables and denotes the vector of right-hand sides of (13), (21), and (22). In the corrected step, the updated value of 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 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 Section 4.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 , the viscosity , and the speed of sound . The initial geometry of the water drop is a circle of radius with its center located at the origin (, ). There is no external forces but initial velocity field , with and the initial pressure field . The number of fluid particles is and corresponding to the initial distance , and the time step . During the stretching process of water drop, the water drop remains elliptical shape and the value of ( is the semiminor axis and is the semimajor axis) remains constant. We let “” denote the interval time step of the periodic density reinitialization of PDRI-SPH, and the PDRI-SPH method becomes the SPH method if .

Figure 2 shows the positions of 1961 particles calculated by PDRI-SPH method with different interval time steps at the time 0.01 s. The distributed particles of using PDRI-SPH () are all more uniform and the outer surfaces are all far smoother than the ones of SPH method (). From Figure 2, we can observe that the better results belong to and . 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 () 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 () 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 Figures 24, and considering the computational cost and the effect of PDRI-SPH with different “num,” we choose the interval time step 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 Figure 5(a). When a droplet impacts on the inclined rigid plate, the shape of the droplet distorts and spreads symmetrically () or asymmetrically () relative to the point of impact, which is shown in Figure 5(b). We define the positive value of the elongation and the negative value of the elongation (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).

All the physical parameter values for a Newtonian droplet are chosen as follows. Its initial diameter and velocity are and , respectively. The total viscosity is , the reference density is , the speed of sound is , and the gravitational force acts downwards with : in this subsection, fluid particles, wall particles, and ghost particles. The initial spacing and the time-step . The height of dropping is from the center of drop to the center () of inclined rigid wall (see Figure 5(a)). For the Cross model droplet, , , 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. The particles are more uniformly distributed than those without using it. The numerical accuracy and stability are improved. 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 (, ) is adopted.

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

Figures 7 and 8 show the effect of the artificial stress for simulating a Newtonian/Cross model droplet impact on horizontal/inclined () 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 (), 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 , 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. The problem of tensile instability occurs more evidently for the Cross droplet than the Newtonian case when the droplet impacts on horizontal rigid plate. A droplet impacts on rigid plate fracture more likely at low impact angles than . In fact, the tensile instability is also related to the ratio of the kinematic viscosity , for a Cross droplet impact on plate. Here, we can find that the fracture is avoided completely by increasing the value of up to 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.

Figure 7: The shape of a Newtonian droplet impact on horizontal (left) or inclined (, right) rigid plate obtained by PDRI-SPH method with , and different artificial stress parameter , at dimensionless time .
Figure 8: The shape of a Cross droplet impact on horizontal (left) or inclined (, right) rigid plate obtained by PDRI-SPH method with , and different artificial stress parameter , at dimensionless time .
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 (, ), artificial stress (), and boundary condition for simulating a Cross droplet impact on inclined rigid plate at different low impact angles , respectively. For the purpose of comparison, the Newtonian droplet case is also considered. The number of fluid particles is set to , corresponding to the initial spacing , wall particles, and ghost particles. The time-step is . 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 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 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 ) 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.

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 are more accurate than those for by observing Figures 10(a) and 11(a), which implies that it is credible to adopt the 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 length and particles number (along the -axis direction): (a1)-(a2) Newtonian droplet; (b1)-(b2) Cross droplet (which corresponds to Figure 10).

We can observe the shapes of a Newtonian drop spreading over inclined rigid surfaces () 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 is similar to that of impact angle ; 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 of (first column), (second column), and (third column) obtained using PDRI-SPH, varying with short dimensionless time.

From Figure 13, we can get that the shapes of the Cross drop spreading over inclined rigid surfaces () 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 1214. 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 () because of the shear-thinning of the Cross fluid.

Figure 13: The shape of a Cross droplet at impact angles of (first column), (second column), and (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.

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 and smoothing length . Observing Figure 15, the numerical results for are very close to those for , 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 angle and smoothing length obtained using PDRI-SPH method.

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.

Conflict of Interests

The authors declare that there is no conflict of interests regarding the publication of this paper.

Acknowledgments

Project is supported by the Postdoctoral Science Foundation of China (Grant no. 2014M550310), the Natural Science Foundation of Jiangsu Province, China (Grant no. BK20130436), and the Innovation Cultivation Funds of Yangzhou University, China (Grant no. 2013CXJ003).

References

  1. F. H. Harlow, “The particle-in-cell method for numerical solution of problems in fluid dynamics,” in Proceedings of the Symposium in Applied Mathematics, vol. 15, pp. 269–288, 1963.
  2. F. H. Harlow and J. E. Welch, “Numerical calculation of time-dependent viscous incompressible flow of fluid with free surface,” Physics of Fluids, vol. 8, no. 12, pp. 2182–2189, 1965. View at Publisher · View at Google Scholar · View at Scopus
  3. C. W. Hirt and B. D. Nichols, “Volume of fluid (VOF) method for the dynamics of free boundaries,” Journal of Computational Physics, vol. 39, no. 1, pp. 201–225, 1981. View at Publisher · View at Google Scholar · View at Scopus
  4. S. Osher and J. A. Sethian, “Fronts propagating with curvature-dependent speed: algorithms based on Hamilton-Jacobi formulations,” Journal of Computational Physics, vol. 79, no. 1, pp. 12–49, 1988. View at Publisher · View at Google Scholar · View at MathSciNet · View at Scopus
  5. M. Do-Quang and G. Amberg, “Numerical simulation of the coupling problems of a solid sphere impacting on a liquid free surface,” Mathematics and Computers in Simulation, vol. 80, no. 8, pp. 1664–1673, 2010. View at Publisher · View at Google Scholar · View at Zentralblatt MATH · View at MathSciNet · View at Scopus
  6. L. W. Zhang, Y. J. Deng, and K. M. Liew, “An improved element-free Galerkin method for numerical modeling of the biological population problems,” Engineering Analysis with Boundary Elements, vol. 40, pp. 181–188, 2014. View at Publisher · View at Google Scholar · View at Zentralblatt MATH · View at MathSciNet · View at Scopus
  7. L. W. Zhang and K. M. Liew, “An improved moving least-squares Ritz method for two-dimensional elasticity problems,” Applied Mathematics and Computation, vol. 246, pp. 268–282, 2014. View at Publisher · View at Google Scholar · View at MathSciNet
  8. L. W. Zhang, P. Zhu, and K. M. Liew, “Thermal buckling of functionally graded plates using a local Kriging meshless method,” Composite Structures, vol. 108, no. 1, pp. 472–492, 2013. View at Publisher · View at Google Scholar · View at Scopus
  9. R. A. Gingold and J. J. Monaghan, “Smoothed particle hydrodynamics theory and application to non-spherical stars,” Monthly Notices of the Royal Astronomical Society, vol. 181, pp. 375–389, 1977. View at Google Scholar
  10. S. Li, D. Qian, W. K. Liu, and T. Belytschko, “A meshfree contact-detection algorithm,” Computer Methods in Applied Mechanics and Engineering, vol. 190, no. 24-25, pp. 3271–3292, 2001. View at Publisher · View at Google Scholar · View at Zentralblatt MATH · View at MathSciNet · View at Scopus
  11. S. Li and W. K. Liu, “Mesh-free particle methods and their applications,” Applied Mechanics Review, vol. 54, pp. 1–34, 2002. View at Google Scholar
  12. V. P. Nguyen, T. Rabczuk, S. Bordas, and M. Duflot, “Meshless methods: a review and computer implementation aspects,” Mathematics and Computers in Simulation, vol. 79, no. 3, pp. 763–813, 2008. View at Publisher · View at Google Scholar · View at Zentralblatt MATH · View at MathSciNet · View at Scopus
  13. L. B. Lucy, “A numerical approach to the testing of the fission hypothesis,” The Astronomical Journal, vol. 83, pp. 1013–1024, 1977. View at Google Scholar
  14. J. J. Monaghan, “Simulating Free Surface Flows with SPH,” Journal of Computational Physics, vol. 110, no. 2, pp. 399–406, 1994. View at Publisher · View at Google Scholar · View at Zentralblatt MATH · View at Scopus
  15. J. P. Morris, P. J. Fox, and Y. Zhu, “Modeling low Reynolds number incompressible flows using SPH,” Journal of Computational Physics, vol. 136, no. 1, pp. 214–226, 1997. View at Publisher · View at Google Scholar · View at Zentralblatt MATH · View at Scopus
  16. S. J. Watkins, A. S. Bhattal, N. Francis, J. A. Turner, and A. P. Whitworth, “A new prescription for viscosity in smoothed particle Hydro dynamics,” Astronomy and Astrophysics Supplement Series, vol. 119, no. 1, pp. 177–187, 1996. View at Publisher · View at Google Scholar · View at Scopus
  17. S. J. Cummins and M. Rudman, “An SPH projection method,” Journal of Computational Physics, vol. 152, no. 2, pp. 584–607, 1999. View at Publisher · View at Google Scholar · View at MathSciNet · View at Scopus
  18. S. Shao and E. Y. M. Lo, “Incompressible SPH method for simulating Newtonian and non-Newtonian flows with a free surface,” Advances in Water Resources, vol. 26, no. 7, pp. 787–800, 2003. View at Publisher · View at Google Scholar · View at Scopus
  19. J. J. Monaghan and A. Kocharyan, “SPH simulation of multi-phase flow,” Computer Physics Communications, vol. 87, no. 1-2, pp. 225–235, 1995. View at Publisher · View at Google Scholar · View at Scopus
  20. J. P. Morris, “Simulating surface tension with smoothed particle hydrodynamics,” International Journal for Numerical Methods in Fluids, vol. 33, no. 3, pp. 333–353, 2000. View at Publisher · View at Google Scholar
  21. R. Ata and A. Soulaïmani, “A stabilized SPH method for inviscid shallow water flows,” International Journal for Numerical Methods in Fluids, vol. 47, no. 2, pp. 139–159, 2005. View at Publisher · View at Google Scholar · View at MathSciNet · View at Scopus
  22. L. Oger and S. B. Savage, “Smoothed particle hydrodynamics for cohesive grains,” Computer Methods in Applied Mechanics and Engineering, vol. 180, no. 1-2, pp. 169–183, 1999. View at Publisher · View at Google Scholar · View at Zentralblatt MATH · View at Scopus
  23. M. Ellero, M. Kröger, and S. Hess, “Viscoelastic flows studied by smoothed particle dynamics,” Journal of Non-Newtonian Fluid Mechanics, vol. 105, no. 1, pp. 35–51, 2002. View at Publisher · View at Google Scholar · View at Scopus
  24. M. Ellero and R. I. Tanner, “SPH simulations of transient viscoelastic flows at low Reynolds number,” Journal of Non-Newtonian Fluid Mechanics, vol. 132, no. 1–3, pp. 61–72, 2005. View at Publisher · View at Google Scholar · View at Scopus
  25. J. Fang, R. G. Owens, L. Tacher, and A. Parriaux, “A numerical study of the SPH method for simulating transient viscoelastic free surface flows,” Journal of Non-Newtonian Fluid Mechanics, vol. 139, no. 1-2, pp. 68–84, 2006. View at Publisher · View at Google Scholar · View at Scopus
  26. W. Benz, “Smooth particle hydrodynamics: a review,” in The Numerical Modelling of Nonlinear Stellar Pulsations: Problems and Prospects, J. R. Buchler, Ed., vol. 302 of NATO ASI Series, pp. 269–288, Kluwer Academic Publishers, Boston, Mass, USA, 1990. View at Google Scholar
  27. J. K. Chen and J. E. Beraun, “A generalized smoothed particle hydrodynamics method for nonlinear dynamic problems,” Computer Methods in Applied Mechanics and Engineering, vol. 190, no. 1-2, pp. 225–239, 2000. View at Publisher · View at Google Scholar · View at Zentralblatt MATH · View at Scopus
  28. J. W. Swegle, D. L. Hicks, and S. W. Attaway, “Smoothed particle hydrodynamics stability analysis,” Journal of Computational Physics, vol. 116, no. 1, pp. 123–134, 1995. View at Publisher · View at Google Scholar · View at MathSciNet · View at Scopus
  29. M. F. Tomé, B. Duffy, and S. McKee, “A numerical technique for solving unsteady non-Newtonian free surface flows,” Journal of Non-Newtonian Fluid Mechanics, vol. 62, no. 1, pp. 9–34, 1996. View at Publisher · View at Google Scholar · View at Scopus
  30. M. F. Tome, L. Grossia, A. Casteloa et al., “A numerical method for solving three-dimensional generalized Newtonian free surface flows,” Journal of Non-Newtonian Fluid Mechanics, vol. 123, pp. 85–103, 2004. View at Google Scholar
  31. M. B. Liu and G. R. Liu, “Restoring particle consistency in smoothed particle hydrodynamics,” Applied Numerical Mathematics, vol. 56, no. 1, pp. 19–36, 2006. View at Publisher · View at Google Scholar · View at MathSciNet · View at Scopus
  32. J. J. Monaghan, “Smoothed particle hydrodynamics,” Reports on Progress in Physics, vol. 68, no. 8, pp. 1703–1759, 2005. View at Publisher · View at Google Scholar · View at MathSciNet · View at Scopus
  33. J. Fang, A. Parriaux, M. Rentschler, and C. Ancey, “Improved SPH methods for simulating free surface flows of viscous fluids,” Applied Numerical Mathematics, vol. 59, no. 2, pp. 251–271, 2009. View at Publisher · View at Google Scholar · View at MathSciNet · View at Scopus
  34. J. J. Monaghan, “Smoothed particle hydrodynamics,” Annual Review of Astronomy and Astrophysics, vol. 30, no. 1, pp. 543–574, 1992. View at Publisher · View at Google Scholar · View at Scopus
  35. J. P. Gray, J. J. Monaghan, and R. P. Swift, “SPH elastic dynamics,” Computer Methods in Applied Mechanics and Engineering, vol. 190, no. 49-50, pp. 6641–6662, 2001. View at Publisher · View at Google Scholar · View at Zentralblatt MATH · View at Scopus
  36. J. J. Monaghan, “SPH without a tensile instability,” Journal of Computational Physics, vol. 159, no. 2, pp. 290–311, 2000. View at Publisher · View at Google Scholar · View at Scopus
  37. J. J. Monaghan and A. Kos, “Solitary waves on a cretan beach,” Journal of Waterway, Port, Coastal and Ocean Engineering, vol. 125, no. 3, pp. 145–154, 1999. View at Publisher · View at Google Scholar · View at Scopus
  38. M. Yildiz, R. A. Rook, and A. Suleman, “SPH with the multiple boundary tangent method,” International Journal for Numerical Methods in Engineering, vol. 77, no. 10, pp. 1416–1438, 2009. View at Publisher · View at Google Scholar · View at MathSciNet · View at Scopus
  39. M. F. Tomé, N. Mangiavacchi, J. A. Cuminato, A. Castelo, and S. McKee, “A finite difference technique for simulating unsteady viscoelastic free surface flows,” Journal of Non-Newtonian Fluid Mechanics, vol. 106, no. 2-3, pp. 61–106, 2002. View at Publisher · View at Google Scholar · View at Zentralblatt MATH · View at Scopus
  40. S. F. Lunkad, V. V. Buwa, and K. D. P. Nigam, “Numerical simulations of drop impact and spreading on horizontal and inclined surfaces,” Chemical Engineering Science, vol. 62, no. 24, pp. 7214–7224, 2007. View at Publisher · View at Google Scholar · View at Scopus
  41. Š. Šikalo, C. Tropea, and E. N. Ganić, “Impact of droplets onto inclined surfaces,” Journal of Colloid and Interface Science, vol. 286, no. 2, pp. 661–669, 2005. View at Publisher · View at Google Scholar · View at Scopus