Open AccessArticle

A Symmetry Particle Method towards Implicit Non‐Newtonian Fluids

Yalan Zhang

Xiaojuan Ban

^*,

Xiaokun Wang

and

Xing Liu

University of Science and Technology Beijing, Xueyuan Road 30, Haidian District, Beijing, China, 100083

Author to whom correspondence should be addressed.

Symmetry 2017, 9(2), 26; https://doi.org/10.3390/sym9020026

Submission received: 17 January 2017 / Accepted: 10 February 2017 / Published: 17 February 2017

(This article belongs to the Special Issue Symmetry in Cooperative Applications II)

Download

Browse Figures

Figure 1
Viscous forces model. V: velocity. "> 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
A falling water column with our method. (a) is Newtonian fluid; (b) is the non-Newtonian fluid. "> Figure 4
The free fall of a water column. (a) is Newtonian fluid (our method); (b) is the non-Newtonian fluid (Andrade et al. [<a href="#B25-symmetry-09-00026" class="html-bibr">25</a>]) with <math display="inline"> <semantics> <mrow> <mo>Δ</mo> <mi>t</mi> <mo>=</mo> <mn>2.5</mn> <mo>×</mo> <msup> <mrow> <mn>10</mn> </mrow> <mrow> <mo>−</mo> <mn>6</mn> </mrow> </msup> <mo> </mo> <mi>s</mi> </mrow> </semantics> </math>; (c) is the non-Newtonian fluid (Andrade et al. [<a href="#B25-symmetry-09-00026" class="html-bibr">25</a>]) with <math display="inline"> <semantics> <mrow> <mo>Δ</mo> <mi>t</mi> <mo>=</mo> <mn>3.0</mn> <mo>×</mo> <msup> <mrow> <mn>10</mn> </mrow> <mrow> <mo>−</mo> <mn>4</mn> </mrow> </msup> <mo> </mo> <mi>s</mi> </mrow> </semantics> </math>; (d) is the non-Newtonian fluid (our method) with <math display="inline"> <semantics> <mrow> <mo>Δ</mo> <mi>t</mi> <mo>=</mo> <mn>3.0</mn> <mo>×</mo> <msup> <mrow> <mn>10</mn> </mrow> <mrow> <mo>−</mo> <mn>4</mn> </mrow> </msup> <mo> </mo> <mi>s</mi> </mrow> </semantics> </math>. "> Figure 5
A water block falling in a pool. (a) is Newtonian fluid; (b) is the non-Newtonian fluid. ">

Versions Notes

Abstract

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.

Keywords:

symmetry; cooperation visualization; smoothed particle hydrodynamics (SPHs); non-Newtonian fluid; implicit method

1. 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 parameters for each particle. Furthermore, an implicit viscosity solver is proposed for non-Newtonian SPH fluid.

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

3. Material Model

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

\frac{d ρ}{d t} = - ρ \nabla \cdot v

(1)

\frac{d v}{d t} = - \frac{1}{ρ} \nabla p + \frac{1}{ρ} \nabla \cdot σ + \frac{1}{m} F^{e x t}

(2)

where

t

denotes the time,

ρ

is the density,

p

is the pressure,

σ

is the viscous stress tensor and

\frac{1}{m} F^{e x t}

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

- \frac{1}{ρ} \nabla p

is related to particle acceleration due to pressure changes in the fluid. The term

\frac{1}{ρ} \nabla \cdot σ

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

γ = \nabla v + {(\nabla v)}^{T}

is the strain rate tensor. The shear rate

γ

is approximated by using Frobenius norm

γ = \sqrt{γ : γ}

to model non-Newtonian fluids.

3.2. 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:

μ (γ) = μ_{\infty} + (μ_{0} - μ_{\infty}) {(1 + {(K γ)}^{α})}^{\frac{n - 1}{α}}

(3)

where

K

controls the shear rate and

α

governs smoothness of transitions,

μ_{0}

and

μ_{\infty}

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

μ_{\infty}

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} - μ_{\infty}

as the shear rate increases, and the simulated fluid is shear thickening non-Newtonian fluid.

3.3. 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:

〈 A (x) 〉 = \sum_{j} \frac{m_{j}}{ρ_{j}} A_{j} W_{i j}

(4)

where

W_{i j} = 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.

In the SPH interpolation scheme, the derivatives of field quantities only affect the smoothing kernel. The gradient of

A

can be calculated as:

\nabla A (x) \approx 〈 \nabla A (x) 〉 = \sum_{j} \frac{m_{j}}{ρ_{j}} A_{j} \nabla W_{i j}

(5)

By substituting density into Equation (5), the density at location

x_{i}

is:

ρ_{i} = \sum_{j} m_{j} W_{i j}

(6)

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:

p_{i} = k (ρ_{i} - ρ_{0})

(7)

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

- \frac{1}{ρ} \nabla p

into Equation (6):

- \frac{1}{ρ_{i}} \nabla p_{i} = - \sum_{j} m_{j} \frac{p_{i} + p_{j}}{ρ_{j}} \nabla W_{i j}

(8)

In order to compute the shear stress

σ_{i}

at particle

i

, the same SPH approximation as [33] is adopted for the deformation tensor

γ_{i} = \nabla v_{i} + {(\nabla v_{i})}^{T}

with:

\nabla v_{i} = \sum_{j} \frac{m_{j}}{ρ_{j}} v_{j i} \otimes \nabla W_{i j}

(9)

where

v_{j i} = v_{i} - v_{j}

. After updating shear stress at all particles, the viscous acceleration is approximated by:

\frac{1}{ρ} \nabla \cdot σ_{i} = \sum_{j} m_{j} (\frac{σ_{i}}{{ρ_{i}}^{2}} + \frac{σ_{j}}{{ρ_{j}}^{2}}) \nabla W_{i j}

(10)

4. 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^{e x t}

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:

v_{i} (t + Δ t) = v_{i} (t) + Δ t \frac{F_{i}^{e x t} (t) + F_{i}^{p} (t) + F_{i}^{v i s} (t)}{m_{i}}

with unknown pressure forces

F_{i}^{p}

and viscosity forces

F_{i}^{v i s}

and known external forces

F_{i}^{e x t}

such as gravity and surface tension. Following the projection concept, intermediate (predicted) velocities are set to be

{v_{i}}^{*}

{v_{i}}^{*} = v_{i} (t) + Δ t \frac{F_{i}^{a d v} (t)}{m_{i}}

(11)

The second intermediate velocity

v^{* *}

is:

{v_{i}}^{* *} = {v_{i}}^{*} + Δ t \frac{F_{i}^{P} (t)}{m_{i}}

(12)

The velocity in the next time step can be written as:

v_{i} (t + Δ t) = {v_{i}}^{* *} + \frac{Δ t}{ρ} \nabla \cdot σ

(13)

The pressure gradient is computed using the SPH formulation [29]:

\nabla p_{i} = k_{i} \nabla ρ_{i} = k_{i} \nabla \sum_{j} m_{j} \nabla W_{i j}

(14)

The pressure force of particle

i

caused by pressure acceleration is determined by:

F_{i}^{p} = - m_{i} \frac{1}{ρ} \nabla p = - k_{i} \frac{m_{i}}{ρ_{i}} \sum_{j} m_{j} \nabla W_{i j}

(15)

F_{j \leftarrow i}^{p}

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

F_{i}^{p} + \sum_{j} F_{j \leftarrow i}^{p} = 0

. Then we can get:

F_{j \leftarrow i}^{p} = k_{i} \frac{m_{i}}{ρ_{i}} m_{j} \nabla W_{i j}

(16)

The velocity changes

Δ v_{i}

by pressure forces is:

Δ v_{i} = Δ t \frac{F_{i}^{P} (t)}{m_{i}} = - Δ t \frac{k_{i}}{ρ_{i}} \sum_{j} m_{j} \nabla W_{i j}

(17)

Similarly:

Δ v_{j} = Δ t \frac{F_{j \leftarrow i}^{p} (t)}{m_{i}} = Δ t \frac{k_{i}}{ρ_{i}} m_{j} \nabla W_{i j}

(18)

According to [28], the intermediate density caused by intermediate velocity is:

{ρ_{i}}^{*} = ρ_{i} (t) + Δ t \sum_{j} m_{j} ({v^{*}}_{i} - {v^{*}}_{j}) \nabla W_{i j}

(19)

We now search for pressure forces to resolve the compression, i.e., the deviation from the rest density

Δ ρ_{i}

\begin{array}{l} Δ ρ_{i} = Δ t \sum_{j} m_{j} (Δ v_{i} - Δ v_{j}) \nabla W_{i j} \\ = Δ t \sum_{j} m_{j} (- Δ t \frac{k_{i}}{ρ_{i}} \sum_{j} m_{j} \nabla W_{i j} - Δ t \frac{k_{i}}{ρ_{i}} m_{j} \nabla W_{i j}) \nabla W_{i j} \\ = - \frac{k_{i}}{ρ_{i}} Δ t^{2} ({(\sum_{j} m_{j} \nabla W_{i j})}^{2} + \sum_{j} {(m_{j} \nabla W_{i j})}^{2}) \\ = ρ_{0} - {ρ_{i}}^{*} \end{array}

(20)

Solving for

k_{i}

yields:

k_{i} = \frac{ρ_{0} - {ρ_{i}}^{*}}{Δ t^{2}} α_{i}, with α_{i} = \frac{ρ_{i}}{({(\sum_{j} m_{j} \nabla W_{i j})}^{2} + \sum_{j} {(m_{j} \nabla W_{i j})}^{2})}

(21)

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)

{ρ_{i}}^{*} = {ρ_{i}}^{*} = \frac{k_{i}}{α_{i}} Δ t^{2}

end for
end while
for all particles

i

{v_{i}}^{* *} = {v_{i}}^{*} - Δ t \sum_{j} m_{j} (\frac{k_{i}}{ρ_{i}} + \frac{k_{j}}{ρ_{i}}) \nabla W_{i j}

end for

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}}^{* *}

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

Algorithm 2. Implicit non-Newtonian SPH fluid solver

for all particles

i

do
find neighborhoods

N_{i}

compute time step size

Δ t