A Symmetry Particle Method towards Implicit Non-Newtonian Fluids

In this paper, a symmetry particle method, the smoothed particle hydrodynamics (SPHs) method, is extended to deal with non-Newtonian fluids. First, the viscous liquid is modeled by a non-Newtonian fluid flow and the variable viscosity under shear stress is determined by the Carreau-Yasuda model. Then a pressure correction method is proposed, by correcting density error with individual stiffness parameters for each particle, to ensure the incompressibility of fluid. Finally, an implicit method is used to improve efficiency and stability. It is found that the non-Newtonian behavior can be well displayed in all cases, and the proposed SPH algorithm is stable and efficient.


Introduction
Non-Newtonian fluid is widely presented in industry, mining engineering, food processing, and many other fields in cases such as pit filled paste, blood, honey, etc.However, the physical nature of non-Newtonian fluid is very complicated.It is pretty difficult to exactly simulate non-Newtonian fluid.A non-Newtonian fluid has characteristic properties of both solids and fluids [1].Take the shear thinning non-Newtonian fluid for example, when it flows with low velocity, its viscosity becomes higher, and the fluid starts to exhibit properties of solids; while as it flows with high velocity, its viscosity gets smaller, and the non-Newtonian fluid exhibits liquid properties.In recent years, smoothed particle hydrodynamics (SPHs) [2] has become an important Lagrangian method for industrial flow simulation [3] and computational fluid dynamics [4,5].Because of the symmetry of SPH method with respect to translations, rotations, and time, the exact conservation of momentum, angular momentum and energy are guaranteed [6][7][8][9].There are several difficulties to simulating non-Newtonian fluids.Firstly, how to model the non-Newtonian fluid.Due to the non-linear relationship between viscosity and shear force, the adequate physical equations and stable approximate simulation method should be re-established for non-Newtonian fluids.Secondly, how to ensure the incompressibility of non-Newtonian fluids.SPH method has some inherent problems which are the same as massless methods [10].Enforcing incompressibility for particle methods is a challenging problem [11].Thirdly, how to improve the efficiency of simulation.Due to the complicated physical equations of non-Newtonian fluids, a small time step is necessary, which proves to be rather time-consuming.
A novel non-Newtonian fluid simulation method for SPH is proposed in this article.The viscous liquid is modeled by a non-Newtonian fluid flow, and the variable viscosity under shear stress is achieved using a viscosity model known as the Carreau-Yasuda model.To improve numerical stability, a pressure correction method is used to correct density error by setting up individual stiffness Symmetry 2017, 9, 26; doi:10.3390/sym9020026www.mdpi.com/journal/symmetryparameters for each particle.Furthermore, an implicit viscosity solver is proposed for non-Newtonian SPH fluid.

Related Works
There are two major approaches to simulate non-Newtonian fluid: Eulerian method and Lagrangian method.The Lagrangian method treats fluid as a particle system.Each point in the fluid is labeled as a separate particle, with a position and a velocity.Instead of tracking each particle, the Eulerian method can look at fixed points in space and see how the fluid quantities (such as density, velocity, etc.) measured at those points change in time [12].
For the Eulerian method, the main problems are how to model non-Newtonian fluid and how to deal with free surface.Goatskin et al. [13] proposed a quasi-linear plastic model to control the change of viscosity in the gradual transformation from a solid to a non-Newtonian fluid.This model meets the von Mises's yield condition.They simulated the melting of a piece of wax with this model by using explicit Euler method.Losasso et al. [14] modeled unified fluids with different viscoelasticities and densities by extending the particle level set.There are many different level sets in the same grid region.The exchanges between different level sets can be achieved by adding the pressure correction process before calculating viscosity.Batty and Bridson [15] found out that the Laplace operator and Strain rate items of non-Newtonian fluid should be computed under the velocity field with no divergence; therefore, the pressure correction process was added after the calculation of viscosity.Besides, they proposed an implied unconditionally stable simulation method based on marker-and-cell.Even though they proposed high precision free surface boundary conditions, it still cannot simulate high viscosity Newtonian fluid.Bergou et al. [16] proposed a method based on discrete differential geometry to simulate one dimensional linear viscous fluid.This approach can simulate one dimensional phenomenon vividly, but distortion often appeared when the flow layer got thicker.Batty and Uribe [17] simulated thin viscous layer with dimensionality deduce technology.They combined nonlinear surface tension and smallest-discrete-surface-area-based formulas and reserved the mass of unit triangular mesh with the local grid reconstruction technology.This approach can simulate the physical conditions of viscous thin layer, such as bending and draping.
Particle-based methods are superior when free surface and complex interface are to be dealt with.So mesh free methods represented by SPH method have been more commonly applied in the field of non-Newtonian fluid simulation [18,19].Müller et al. [20] simulated fluids of high elasticity and high plasticity by structuring a point-based animation modeling method.He combined continuum mechanics with von Mises' yield condition and calculated the displacement and the velocity using moving least squares.Afonso et al. [21] also completed simulations with a General Newtonian Fluid model, and the viscosity of fluid was controlled by the Jump Number.The melting process of a heated non-Newtonian fluid [22], which later becomes a low viscosity fluid, was also modeled by the general Newtonian fluid model.Rafiee et al. [23] proposed a free surface SPH method.This method can be applied to Newton fluid and viscoplastic fluid, and simulate the bending phenomena of viscous fluids.The relationship between the pressure and the viscosity of non-Newtonian fluids was described with the real polymer state equation [24].Andrade et al. [25] can simulate the flow of polymer with high viscosity and high pressure.In addition, the unified molding of viscous Newtonian fluid and shear thinning non-Newtonian fluid were simulated, but this method can only simulate fluid with low viscosity; whereas unstable phenomenon would appear when the time step is larger.Zhang et al. [26] proposed an explicit method for shear thinning non-Newtonian fluid.
In addition to Euler and Lagrange methods, recently there is a novel alternative, the codimensional method [27], which is able to simulate a non-Newtonian fluid.This method enables the simulation of interaction between different dimensions of a non-Newtonian fluid, such as filamentous remaining marks left on a pool of paint after getting brushed over, which is difficult for Euler and Lagrange methods to achieve.
For Newtonian fluid simulation, the standard SPH has been gradually replaced by a predictive-corrective method to ensure incompressibility with larger time step.Predictive-corrective algorithm was first proposed by Solenthaler and Pajarola [28], which is to predict the state of a fluid in absence of pressure, then use pressure for its conduct correction.Ihmsen et al. [29] proposed an implicit SPH method that adjusts fluid's density by pressure according to the relationship between density and velocity in the next time step.Bender and Koschier [30] proposed a divergence-free velocity method, which is essentially also a predictive-corrective method that prevents volume compression and enforces a divergence-free velocity field formulation of the physical laws.However, none of them involved non-Newtonian fluid simulation.

Governing Equations
In the SPH scheme, the fluid is described as moving particles with physical attributes such as position x, velocity v, and mass m.Mass is used to determine the local geometry of the fluid and assumed to be constant.Non-Newtonian fluids are a continuous medium whose motion equations must follow the mass and momentum conservation principles.The motion equations are represented by simulating the advection of particles.Through the use of integral interplant, the dependent field variables are expressed by integrals which are approximated by summation interplant over neighboring particles.
Compared with Newtonian fluid [31], the strain tensor term should be added in non-Newtonian fluid.In the Lagrangian frame, the governing equations for the incompressible fluid with spatially constant viscosity can be written as the following where t denotes the time, ρ is the density, p is the pressure, σ is the viscous stress tensor and 1 m F ext denotes the acceleration due to external forces such as gravity, surface tension, and friction.
The Equation ( 1) is the continuity equation.Since the mass is carried by particles, the mass is conserved exactly.The Equation ( 2) is the momentum equation.The term − 1 ρ ∇p is related to particle acceleration due to pressure changes in the fluid.The term 1 ρ ∇ • σ describes the viscous acceleration.Viscosity, which is defined as the resistance of a fluid to flow, is essentially friction forces (Figure 1).For incompressible flow, the viscous stress tensor is σ = µγ, where µ is the viscosity and γ = ∇v + (∇v) T is the strain rate tensor.The shear rate γ is approximated by using Frobenius norm γ = .

Carreau-Yasuda Model
Both Newtonian and non-Newtonian fluid flows are considered in this paper, and the variable

Carreau-Yasuda Model
Both Newtonian and non-Newtonian fluid flows are considered in this paper, and the variable viscosity under shear stress is achieved using a rheological model known as Carreau-Yasuda model [32].The relationship between shear rate and viscosity is defined by: where K controls the shear rate and α governs smoothness of transitions, µ 0 and µ ∞ are given positive constants.When n = 1, the fluid is a Newtonian fluid with viscosity µ 0 .When K = 0, the kinematic viscosity is equivalent to µ 0 , which means under low shear rate, the viscosity of non-Newtonian fluid is governed by µ 0 .When n < 1, the second term approaches zero and the viscosity approaches to µ ∞ as the shear rate increases, the simulated fluid is shear thinning non-Newtonian fluid.When n > 1, the viscosity increases with the amount of scaled by µ 0 − µ ∞ as the shear rate increases, and the simulated fluid is shear thickening non-Newtonian fluid.

Smoothed Particle Hydrodynamics
In the SPH approach, the basic idea is to discretize the fluid domain into a finite number of particles with physical attributes (Figure 2).A continuous field A(x i ) can be obtained by summation over all the sampling points x j with mass m j and density p j in the support domain: where W ij = W x − x j , h is the kernel function with h the particle radius, and j iterating all the neighbors.The chosen kernel function is limited by several conditions, such as normalization, compact condition, and delta function property.In this paper, the popular quantic spline kernel is used for all simulations.
Symmetry 2017, 9, 26 5 of 14 In the SPH interpolation scheme, the derivatives of field quantities only affect the smoothing kernel.The gradient of A can be calculated as: By substituting density into Equation ( 5), the density at location i x is: Generally, there are two formulas to relate pressure and density, namely the Poisson equation and the gas state equation.Solving the Poisson equation guarantees the incompressibility of fluid, In the SPH interpolation scheme, the derivatives of field quantities only affect the smoothing kernel.The gradient of A can be calculated as: By substituting density into Equation ( 5), the density at location x i is: Generally, there are two formulas to relate pressure and density, namely the Poisson equation and the gas state equation.Solving the Poisson equation guarantees the incompressibility of fluid, but it is time-consuming.In contrast, the explicit calculation is adopted in an ideal gas equation, but results in the stiffness value having to be very large so that the density fluctuations will not exceed 1%.However, according to the Courant-Friedrichs-Levy (CFL) condition, a large stiffness value will restrict the time step and make the simulation inefficient.In this paper, the equation of the gas state is used for efficiency and the problem of the stiffness value will be solved by pressure projection: where k is a gas constant that depends on the temperature, and ρ 0 is the rest density.The value ρ 0 = 1000 kg/m 3 turns out to be suitable for all experiments.
To compute the pressure acceleration, substitute − 1 ρ ∇p into Equation (6): In order to compute the shear stress σ i at particle i, the same SPH approximation as [33] is adopted for the deformation tensor γ i = ∇v i + (∇v i ) T with: where v ji = v i − v j .After updating shear stress at all particles, the viscous acceleration is approximated by: 1

Pressure Projection Algorithm
In the non-Newtonian fluid model proposed above, due to a rather high compressibility caused by the gas state equation, when a non-Newtonian fluid flows, the volume of the fluid might change [31].In order to address this problem, a pressure projection algorithm is proposed to correct density by pressure.
Currently, there are many different pressure solvers that minimize the density error ρ − ρ 0 [28,29].The algorithm proposed in this paper is similar with [28].Firstly, only external force F ext is considered to obtain the first intermediate velocity v * .Secondly, to enforce fluid incompressibility, a pressure projection algorithm is proposed in this section by using a new solver that sets a separate stiffness parameter k i for each fluid particle i to satisfy local incompressibility, making the density of the entire fluid system remain unchanged and getting the second intermediate velocity v * * .Next, an implicit method is proposed to solve the viscous stress tensor (see Section 5).
Finally, particle velocity and position are integrated using the Euler scheme, and the velocity can be rewritten as: with unknown pressure forces F p i and viscosity forces F vis i and known external forces F ext i such as gravity and surface tension.Following the projection concept, intermediate (predicted) velocities are set to be v i * : The second intermediate velocity v * * is: The velocity in the next time step can be written as: The pressure gradient is computed using the SPH formulation [29]: The pressure force of particle i caused by pressure acceleration is determined by: F p j←i is the pressure forces that act from particle i on the neighboring particles j.According to Newton's third law, its value is equal to the total pressure of all neighboring particles j affects i, which satisfies Then we can get: The velocity changes ∆v i by pressure forces is: Similarly: According to [28], the intermediate density caused by intermediate velocity is: We now search for pressure forces to resolve the compression, i.e., the deviation from the rest density ∆ρ i : Solving for k i yields: Algorithm 1. Pressure Projection Algorithm for all particles i do update ρ i using Equation ( 7) update α i using Equation ( 21) end for for all particles i do update V i * using Equation ( 12) end for while (ρ avg − ρ 0 > η) or iter < 2 do for all particles i do update ρ i * using Equation (20) update k i using Equation ( 21) Algorithm 1 outlines the simulation using our novel method in detail.First, the particle neighborhoods N i using compact hashing [34] are determined.Then for each particle, the density ρ i and the factor α i are computed.Since this factor solely depends on the current positions, it is pre-computed before executing the solvers, which saves the solver from significant computational work.The third step in the simulation loop is to compute the strain tensor term of a non-Newtonian fluid, which is the main difference between a Newtonian fluid and a non-Newtonian fluid.
The forth and the fifth step of the algorithm are predictive-correction steps, whose purpose is to minimize the density error by the pressure solver, scilicet solving ρ − ρ 0 = f (p).In this step, the non-pressure forces are used to compute predicted velocities v i * .The constant density solver uses this prediction and the pre-computed factors α i to determine the pressure forces for each neighborhood in order to correct the density error.Finally, the resulting velocity changes are used to update the velocity of the particle without viscosity force v i * * .

Implicit Viscosity Solver
The viscous acceleration is commonly computed from the divergence of viscous stress tensor, which is mentioned in Section 3.2.In this paper, the proposed solver first computes an intermediate velocity v * without considering pressure and viscosity.Then, the second updated intermediate velocity v * * is computed by using the proposed pressure correction method, which is described in Section 4. Finally, the proposed viscosity solver computes the final velocity v(t + ∆t) and the position of particles.Algorithm 2 summarizes the proposed implicit non-Newtonian SPH fluid solver.
According to Algorithm 2, the viscosity solver starts with the strain rate tensor to compute the final velocity.Substituting Equation ( 15) into the strain rate tensor γ i = ∇v i + (∇v i ) T , and then getting the symmetric tensor for particle i as a six-dimensional vector yields: To solve the viscous acceleration 1 ρ ∇ • σ i , the gradient of viscous stress tensor should be computed: According to the Equation ( 4), µ i is the function of the strain rate tensor γ i , while K, α, µ 0 , and µ ∞ are given positive constants and will not influence the gradient computing.Thus, the gradient of viscous stress tensor can be written as: Due to the formulation, the required derivatives of the strain rate tensor γ i with respect to the velocity v j are determined by the following matrix: In the SPH scheme, the derivatives of the strain rate tensor γ i with respect to the velocity v i can be written as: Now that all the variable quantities in Equation ( 14) can be solved, the final velocity and position of each particle are gained.

Implementation and Results
All timings are given for an Intel 3.50 GHz CPU with four cores.The simulation software is parallelized with Open MP.The simulation and surface reconstruction are actualized with C++ language and multi-threading technology.The density fluctuation was set to 0.01.
Figure 3 shows the Newtonian fluid and non-Newtonian fluid according to our method.As described in Section 3.2, when n = 1, the non-Newtonian fluid model is simplified to a Newtonian fluid with constant kinematic viscosity (Figure 3a).When the fluid touches the container at the bottom of the scene (the container is set to be transparent for observers' convenience), the fluid particles that touch the container splash upward.In contrast, the non-Newtonian fluid particles (shown in Figure 3b) cannot splash when they collide with container due to their viscosity.The setting and statistics are illustrated in Table 1.
Symmetry 2017, 9, 26 10 of 14 In the SPH scheme, the derivatives of the strain rate tensor i γ with respect to the velocity i v can be written as: Now that all the variable quantities in Equation ( 14) can be solved, the final velocity and position of each particle are gained.

Implementation and Results
All timings are given for an Intel 3.50 GHz CPU with four cores.The simulation software is parallelized with Open MP.The simulation and surface reconstruction are actualized with C++ language and multi-threading technology.The density fluctuation was set to 0.01.
Figure 3 shows the Newtonian fluid and non-Newtonian fluid according to our method.As described in Section 3.2, when 1 n  , the non-Newtonian fluid model is simplified to a Newtonian fluid with constant kinematic viscosity (Figure 3a).When the fluid touches the container at the bottom of the scene (the container is set to be transparent for observers' convenience), the fluid particles that touch the container splash upward.In contrast, the non-Newtonian fluid particles (shown in Figure 3b) cannot splash when they collide with container due to their viscosity.The setting and statistics are illustrated in Table 1.   Figure 4 shows the free fall of a water column in particle state with 13,671 particles.In the same initial scene, both Newtonian fluid (Figure 4a) and non-Newtonian fluid (Figure 4b) could be achieved by governing the coefficients.We compared the non-Newtonian fluid simulation of our method with the method of Andrade et al. [25].Andrade et al. used explicit integration without pressure projection.During the experiment, when the scene was small or the fluid particles moved gently, running the algorithm with pressure correction was found to be time and memory-consuming.However, when the scene is large or the fluid particles collide intensely, the method of Andrade et al. will be unstable with the larger time step.As shown in Figure 4c, when the time step was set to 3.0 × 10 −4 s, fluid particles scattered in all directions with the method of Andrade et al. [25], while our algorithm still ran with the pressure projection (Figure 4d).1.14 min Figure 3 shows the free fall of a water column in particle state with 13,671 particles.In the same initial scene, both Newtonian fluid (Figure 4a) and non-Newtonian fluid (Figure 4b) could be achieved by governing the coefficients.We compared the non-Newtonian fluid simulation of our method with the method of Andrade et al. [25].Andrade et al. used explicit integration without pressure projection.During the experiment, when the scene was small or the fluid particles moved gently, running the algorithm with pressure correction was found to be time and memory-consuming.However, when the scene is large or the fluid particles collide intensely, the method of Andrade et al. will be unstable with the larger time step.As shown in Figure 3c, when the time step was set to 3.0 x 10 -4 s, fluid particles scattered in all directions with the method of Andrade et al. [25], while our algorithm still ran with the pressure projection (Figure 3d).      Figure 5 shows a water ball falling in a pool.The first line is Newtonian fluid (Figure 5a) and the second line is the non-Newtonian fluid (Figure 5b).As can be seen from the pictures, the wave of Newtonian fluid lashed against the water, sending up pearly spray.For the non-Newtonian fluid, there is no wave and the water sank into the pool directly.The setting and statistics are illustrated in Table 2.
Figure 5 shows a water ball falling in a pool.The first line is Newtonian fluid (Figure 5a) and the second line is the non-Newtonian fluid (Figure 5b).As can be seen from the pictures, the wave of Newtonian fluid lashed against the water, sending up pearly spray.For the non-Newtonian fluid, there is no wave and the water sank into the pool directly.The setting and statistics are illustrated in Table 2.At our test scene (Figure 4), we compared the simulation times and stability of our method with the method of Andrade et al. [25].The physical parameters and performance measurements are summarized in Table 3.For the explicit method without pressure projection (Andrade et al. [25]), the time step is set according to the CFL condition [35] and the viscosity of fluid [25].In the case of η = 0.01, the time step of our method could be running well with Δt = 3.0x10 -4 s, while the method of Andrade et al. [25] could not.As a result, for a two-second animation, our method needs 32 frames and 44.5 min for computing, while the method of Andrade et al. [25] needs 67 frames and 59.5 min for computing.The simulation parameters and performance of Figure 4 are also listed in Table 3.At our test scene (Figure 4), we compared the simulation times and stability of our method with the method of Andrade et al. [25].The physical parameters and performance measurements are summarized in Table 3.For the explicit method without pressure projection (Andrade et al. [25]), the time step is set according to the CFL condition [35] and the viscosity of fluid [25].In the case of η = 0.01, the time step of our method could be running well with ∆t = 3.0 × 10 −4 s, while the method of Andrade et al. [25] could not.As a result, for a two-second animation, our method needs 32 frames and 44.5 min for computing, while the method of Andrade et al. [25] needs 67 frames and 59.5 min for computing.The simulation parameters and performance of Figure 4 are also listed in Table 3.

Conclusions
In this paper, we proposed a novel symmetry particle-based method for non-Newtonian fluid, where the variable viscosity is decided by the Carreau-Yasuda Model.Our method includes a pressure correction step to reduce the density errors and ensure the incompressibility of fluid.An implicit method is used for the non-Newtonian fluid simulation to improve efficiency and stability.With our method, the fluid simulation runs well in large time steps compared with previous work.
One limitation of our method is that it does not consider temperature in our non-Newtonian fluid model.By taking this into fluid model, one could begin to simulate phenomena of solid-liquid conversion, like wax melt in high temperature.Additionally, work remains to be done on the topic of multi-phase coupling between different non-Newtonian fluids.This might be a potential direction for future work.

Figure 2 .
Figure 2. Illustration for smoothed particle hydrodynamics (SPH) method.It is a symmetry particle method and the physical attributions of each particle is weighted in summation by its neighbors.

Figure 2 .
Figure 2. Illustration for smoothed particle hydrodynamics (SPHs) method.It is a symmetry particle method and the physical attributions of each particle is weighted in summation by its neighbors.

Figure 3 .
Figure 3.A falling water column with our method.(a) is Newtonian fluid; (b) is the non-Newtonian fluid.Figure 3. A falling water column with our method.(a) is Newtonian fluid; (b) is the non-Newtonian fluid.

Figure 3 .
Figure 3.A falling water column with our method.(a) is Newtonian fluid; (b) is the non-Newtonian fluid.Figure 3. A falling water column with our method.(a) is Newtonian fluid; (b) is the non-Newtonian fluid.

Figure 4 .
Figure 4.The free fall of a water column.(a) is Newtonian fluid (our method).(b) is the non-Newtonian fluid (Andrade et al.[25]) with . (c) is the non-Newtonian fluid (Andrade et al.[25]) with , (d) is the non-Newtonian fluid (our method) with

Figure 5 .
Figure 5.A water block falling in a pool.(a)is Newtonian fluid.(b) is the non-Newtonian fluid.

Table 2 .
The simulation parameters for Figure 5. (Both Newtonian and non-Newtonian fluid) Simulation domain size 8 m × 8 m ×

Figure
Figure

Figure 5 .Table 2 .
Figure 5.A water block falling in a pool.(a) is Newtonian fluid; (b) is the non-Newtonian fluid.

Table 1 .
The simulation parameters for Figure3.(Both Newtonian and non-Newtonian fluid).

Table 1 .
The simulation parameters for Figure3.(Both Newtonian and non-Newtonian fluid).

Table 3 .
The simulation parameters and performance.