Mathematical Problems in Engineering

Volume 2015 (2015), Article ID 915973, 17 pages

http://dx.doi.org/10.1155/2015/915973

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

^{1}Department of Mathematics, School of Hydraulic, Energy and Power Engineering, Yangzhou University, Yangzhou 225002, China^{2}Department 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 [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. 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 [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 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.