Parameter determination for the Cross rheology equation and its application to modeling non-Newtonian flows using the WC-MPS method

ABSTRACT

A weakly compressible moving particle semi-implicit (WC-MPS) method is utilized to simulate non-Newtonian free surface flows due to the advantages of particle methods with respect to handling large deformation and fragmentation. The Cross rheology equation was selected in order to capture the viscous features of the mixture flows. To numerically implement the Cross equation, an experiment-based method was proposed to determine the four rheology parameters in the equation. The method of using a WC-MPS model to study non-Newtonian dam break flow problems was then adopted. The capabilities of the proposed method were tested by simulating different materials with the proposed method in modeling non-Newtonian free surface flows. Significant viscous features were reproduced by the proposed model.

KEYWORDS:

• Particle method
• MPS
• non-Newtonian flow
• rheology
• Cross equation

1. Introduction

In engineering applications, the computation of free surface deformation is often challenging, especially for non-Newtonian flows. As a significant example, the failure of a non-Newtonian material column, or dam break, is not only an important issue in understanding debris flows but also a popular topic in fluid mechanics. The deposition of debris flows and runout distance becomes a critical issue (Whipple, 1997), and is crucial for predicting debris behavior in terms of hazard assessment and mitigation (Roche, Montserrat, Nino, & Tamburrino, 2008).

Dam break problems with different kinds of non-Newtonian materials have been widely studied, both experimentally and numerically. For example, Major and Pierson (2010) analyzed debris flow rheology in their experiments for fine-grained slurries; Komatina and Jovanovic (1997) investigated a dam-break-type flow of a set of water-clay mixtures with different solid volume concentrations; Wang, Morgenstern, and Chan (2010) developed a Lagrangian finite difference model to numerically study flow slides and debris flows; Lube, Herbert, Sparks, and Freundt (2011) performed an experimental study for sand column collapses onto a rough inclined channel; Lai and Khan (2012) used one-dimensional shallow water equations to numerically study dam break flows on wet and dry beds; Spinewine and Capart (2013) studied a sheet flow over an erodible granular bed due to a thick and rapidly sheared layer of water and grains that is driven by water; and Pu, Shao, Huang, and Hussain (2013) adopted an improved shallow water equation with a surface gradient upwind method in their incompressible smoothed particle hydrodynamics model to investigate dam break flows along both wet and dry channels.

As an alternative to mesh-based simulation methods, mesh-free methods have been flourishing since their emergence a few decades ago and are expected to be very practical in engineering applications. The core idea and development of mesh-free methods is to establish an alternative for solving numerical equations with accurate and stable solutions. Particles methods offer an easier way to acquire the solutions to problems with free surface, deformable boundaries, and large deformation, to which the conventional grid-based methods are difficult to apply.

The moving particle semi-implicit (MPS) method is one of the mesh-free methods and was first developed by Koshizuka, Tamako, and Oka (1995). It utilizes a set of arbitrarily distributed particles (or nodes) to represent the problem domain without introducing any mesh or grid which requires connectivity between nodal points. The MPS method was originally designed to model fluid flow problems, and it provides approximations of the partial differential flow equations on the basis of integral interpolations. In short, the MPS method offers a set of simple models for the derivatives in the flow equations according to a local weighted averaging algorithm. In the past decade, the MPS method has been applied to a wide range of practical problems related to fluid flows (Gotoh, Ikari, Memita, & Sakai, 2005; Gotoh & Sakai, 2006; Shakibaeina & Jin, 2010; Shibata & Koshizuka, 2007; Sueyoshi, Kashiwagi, & Naito, 2008). The MPS method has also been applied to other engineering problems. For example, Koshizuka and Oka (2001) implemented the MPS method in simulating thermal fragmentation processes, which are common in nuclear engineering problems; Heo, Koshizuka, and Oka (2002) studied shear transfers numerically in mechanical engineering problems; Z. Sun, Xi, and Chen (2009) studied the coalescence and separation mechanisms for the binary collision of liquids, which are related to both mechanical engineering and chemical engineering; and Tsubota et al. (2006) investigated the bioengineering problems associated with the flow of red blood cells and platelets in blood flows.

On the other hand, the MPS method has been further investigated and improved by researchers since its introduction, such as in the investigation of the behaviors of different kernel functions (Ataie-Ashtiani & Farhadi, 2006), momentum conservation studies (Khayyer & Gotoh, 2008), artificial pressure fluctuation analysis (Gotoh & Sakai, 2006), the sub-particle scale (SPS) turbulence model (Gotoh, Shibahara, & Sakai, 2001), and the hybrid particle-mesh method multiphase fluid flows (J. Liu, Koshizuka, & Oka, 2005). More recently, Li, Sun, Chen, Xi, and Liu (2015) developed an improved surface recognition method and collision model for avoiding particle clustering, Q. R. Wu, Tan, and Xing (2014) developed a new fluid-wall interface to replace traditional solid boundary definitions in MPS; and Chen, Xi, and Sun (2014) introduced new surface detection techniques for surface particles.

Recently, Fu and Jin's (2014) study on water dam break flows indicated the advantage of a weakly compressible MPS (WC-MPS) method in tracking the flow profiles and transient flow features. To the current knowledge of the authors, little attention in the MPS model was given to the non-Newtonian flows. To improve the accuracy in the discretization of the divergence of the shear stress tensor, Xiang and Chen (2015) used an SPH (Smoother Particle Hydrodynamics) kernel to simulate non-Newtonian flows in MPS. There are several examples of using SPH to study non-Newtonian flows: Fan, Tanner, and Zhang (2010) studied an implicit SPH model for non-Newtonian flows; Tao, Ren, Lei, and Lu (2014) used a corrected kernel derivative scheme for non-isothermal, non-Newtonian viscous fluid; Pan, Tartakovsky, and Monaghan (2013) studied coupled ice-sheet and ice-shelf dynamics with a non-Newtonian SPH model; Martys, George, Chun, and Lootens (2010) modeled a non-Newtonian fluid with a spatially varying viscosity; Riviere, Khelladi, Farzaneh, Bakir, and Tcharkhtchi (2013) simulated polymer behaviors as non-Newtonian fluid on the mold surface; H. Sun and Han (2010) improved the model for elastic and plastic parameter values; and Basu, Das, Jenetzke, and Green (2011) used Bingham rheology model for mud material to simulate deformable landslide on a mild slope. Since the non-Newtonian materials show very different flow patterns compared to those of Newtonian materials, rheology equations (or constitutive equations) are usually introduced to study their flows. The current study aims at facilitating a rheology model, i.e., the Cross rheology equation, into the WC-MPS model to capture the viscosity changes of the non-Newtonian materials due to the shear stress variation. Moreover, an experiment-based method is proposed to determine the four parameters required in the Cross rheology equation. The method is validated by the measured rheological relations of three different non-Newtonian materials, and the determined parameters together with the Cross equation are then incorporated into the WC-MPS model and validated by experimental tests for different non-Newtonian flows. The simulation results using the parameters determined by the proposed method are compared with measurements as well. The advantage of using WC-MPS with parameters determination technique for Cross model are shown through the comparisons between the numerical and experimental results.

2. Governing equations and MPS theory

The mass and momentum conservation laws are considered in the current WC-MPS method. The corresponding equations are the continuity equation and the Navier–Stokes momentum equation: (1) (2) where ρ is the material density, u is the velocity vector, t is the time, p is the pressure, τ is the shear stress, and f is the external forces. It should be noted that the convection term is directly included in the governing equations since the Lagrangian form of the momentum equation is used. The numerical error induced by the convection term is minimized since the fluid movement is tracked in the Lagrangian particle method.

2.1. Particle interaction

Particle methods pose another computational regime in the flow of materials without the use of a mesh or grid system. The problem domain is discretized and represented by numerous unique and independent particles with pre-determined sizes – thus, the flow process can then be studied by tracking the movement of these particles. These particles carry certain properties of interest, of which velocity and pressure are those which are most commonly examined. In the MPS method, the properties of each particle are influenced by others, and the assessment of the influence can be monitored by a kernel function. In this study, a third-order polynomial spiky function is used as follows: (3) where R_ij=|r_j-r_i| and is the distance between particles i and j, and r_e is the interaction radius surrounding particle i. Obviously, the interaction area surrounding particle i is dominated by the interaction radius r_e, within which all the other particles will be considered as affecting particle i. In other words, all the particles falling in the interaction area of i contribute to the property computation of particle i through the kernel function while considering their distances in relation to particle i. The kernel function is a measurement of the influence which the surrounding particles have on the target particle, and the kernel function is actually a weighted average function to approximate and smooth the physical quantity of a particle through its surroundings. Koshizuka et al. (1995) studied the influence of various r_e values and found that the simulation results are in good agreement with the analytical solution as long as r_e is greater than three times the particle size Δl. Therefore, r_e = 3Δl is employed in the current study.

2.2. MPS discretization

In order to simulate the whole flow process, the governing equations need to be solved over time. Therefore, the differential terms in the governing equations have to be approximated and solved. In the current study, using the MPS method, the approximation of these differential terms is referred to as MPS discretization or MPS differential modeling. The MPS model uses a smoothed quantity to approximate the gradient and divergence term in the governing equation; and the second-order derivative (or the Laplacian term) was approximated through the study of a diffusion model (Koshizuka et al., 1995). The MPS method was originally designed for Computational Fluid Dynamics (CFD), and the discretization model includes an approximation of the first- and second-order derivatives in the governing equations.

The first-order derivative (Koshizuka et al., 1995) is computed as: (4) where φ is a physical property possessed by a particle, d is the dimension, and n₀ is the initial particle number density.

The second-order derivative, or Laplacian term, is defined as: (5) where λ is the correction and is defined as: (6)

2.3. Weak compressibility

In the MPS method, the real fluid density is expressed by the particle number density and defined as: (7) The particle number density is used to identify the number of neighboring particles (the particles in the interaction area of the particle of interest) and it is related to the real fluid density as follows: (8) where m is the mass of the fluid particle. In this study, the weak compressibility model is used to handle particle number density (or density). The use of such a model dramatically decreases the computation time of the whole simulation, since it introduces an equation of state for the pressure calculation in an explicit form (Shakibaeinia & Jin, 2010). The equation of state is defined as: (9) where c₀ is the sound speed in the reference medium, ψ is the polytrophic constant 1 < ψ < 7 (a typical value of 7 is used in the current study; see Crespo, 2008; Narayanaswamy, 2008), <N*>_i is the temporary particle number density, and n₀ is the initial particle number before the flow starts.

In practical simulations, a smaller artificial sound speed is used to avoid small time steps (Shakibaeinia & Jin, 2010). Since the fluid is assumed to be weakly compressible, a temporary change of the particle number density variation is allowed in the simulation. However, the particle number density variation should not be greater than 1%, and c₀ was assigned to be 10 times greater than the maximum fluid velocity (Dalrymple & Roger, 2006).

2.4. MPS solution algorithm

The WC-MPS method used in the current study employs a prediction-correction method for one complete step of time marching.

In order to implement the prediction-correction method, the velocity at the next time step u^t+Δt is split into two parts, i.e., the predicted velocity u* and the velocity correction u'. Then from the finite difference discretization and mathematical transformation Du/Dt in the momentum, Equation (1) can be rewritten as (10) where u^t is the velocity at the current time step. By substituting the above into Equation (1), the following is obtained: (11) Consequently, the discretization of the momentum equation is then divided into two steps: the prediction step and the correction step, respectively: (12) (13) The Courant–Friedrichs–Lewy (CFL) stability condition should be satisfied due to an explicit time-splitting algorithm used in the present model. The CFL condition is a relation between the average particle distance Δl and the time step Δt : (14) where C_CFL is the Courant number with 0 < C_CFL ≤ 1, and c₀ is the sound speed as described in the equation of state. The viscous diffusion (Shao & Lo, 2003) should also be considered in the time step determination since the viscosity is larger in the non-Newtonian flows. However, the sound speed is employed in the calculation instead of the flow velocity normally used in the mesh-based method, the diffusion number factor is not included in the calculation in this study, and time steps are decided through Equation (14) and numerical trials.

2.5. Boundary conditions

Near the solid boundary or the free surface boundary, the particle number density of a given particle will decrease, as there are no particles outside the boundaries. Hence, the specific boundary conditions are required when computation approaches the solid boundary and the free surface boundary. Obviously, they cannot be treated in the same manner since the flow patterns are different between the solid boundary and free surface boundary.

2.5.1. Solid boundary

Close to the solid boundaries, such as the side walls and channel beds, the boundary conditions will influence the flows. As discussed previously, the particle number density decreases near the solid boundary, therefore dummy particles are generated outside the solid boundaries to compensate for the particle deficiency. These dummy particles ensure the correct neighbor search near the solid boundary and contribute to the computation of the properties for the fluid particles near the boundary. However, the dummy particles do not move and remain in their preset locations throughout the simulation. Several layers of dummy particles are defined according to the neighbor search radius, with the number of layers is usually set as the same as the neighbor search radius. The neighbor search radius r_e = 3Δl is used in the present study and therefore three layers of dummy particles are required outside the solid boundary. It should be noted that the pressure of the first layer of dummy particles is computed and then copied to the remaining two layers of dummy particles.

Practical boundary effects can be reflected by employing the dummy particles and assigning the corresponding properties to the dummy particles, such as the commonly used free-slip and no-slip boundary conditions.

2.5.2. Free surface boundary

Similar to the solid boundary, there is no particle outside the free surface boundary. However, the free surface boundary cannot be treated the same as the solid boundary. In MPS, the free surface particles are recognized by a critical particle number density defined as a certain percentage of the initial particle number density. The model is controlled to calculate the particle number density at every time step of the simulation for every particle movement, and the critical particle number density is defined as (15) where β is the threshold coefficient. Obviously, β ranges between 0 and 1. In general, the value of β is between 0.8 and 0.99 (Koshizuka, Nobe, & Oka, 1998), and in the current study, β = 0.96 is used. Figure 1 shows the typical setup of the fluid particles, free surface particles, and solid boundary particles. As previously stated, wall particles are the first layer of dummy particles in Figure 1.

Figure 1. Setup of different particles.

3. Rheology model

Most non-Newtonian fluids show abrupt transformations in flow regimes over a relatively narrow range of the stress they undergo. In other words, these fluids will not show liquid-like behaviors when the stresses are not large enough (i.e., to generate continued deformations). In short, the viscosity of non-Newtonian fluids is shear rate dependent.

Those non-Newtonian fluids that exhibit a decreased viscosity with increasing shear rate are defined as shear-thinning fluids, which account for the majority of non-Newtonian fluid flows. The behaviors of shear-thinning materials have drawn attention over the last few decades due to their common presence in everyday activities and industries (Barnes, 1999). The flows of shear-thinning fluids are usually referred to as pseudoplastic flows.

To study the relationship between deformation and stress, rheology equations are commonly used. In the current WC-MPS study, the Cross rheology equation is used to investigate the non-Newtonian behaviors. In theory, an ideal rheology equation to relate viscosity and shear rate should provide an accurate fit for most experimental measurements over a wide range of the shear rate change (Cross, 1965). Moreover, the equation should include a minimum number of independent constants as well, since these constants will have real physical meaning and must be firmly evaluated accordingly. At least four parameters are required to predict the general flow behaviors of non-Newtonian fluids.

3.1. Cross equation

The Cross rheology equation (Cross, 1965) for non-Newtonian fluids was developed based on the assumption that the pseudoplastic flow is related to the formation and rupture of structural linkages of the materials: (16) where γ is the shear rate, μ₀ is the viscosity when the shear rate is close to zero, μ_∞ is the viscosity when the shear rate is infinity, n is the flow behavior index, and α is the consistency index with the dimension of sⁿ. For shear-thinning materials, the value of n is between 0 and 1.

Under the Bingham model (Shao & Lo, 2003), the fluid at rest can resist the shear stress less than yield stress without deformation. If the shear stress exceeds the yield stress, the fluid starts to flow and behaves as a Newtonian fluid. With the condition of αγⁿ >> 1, n = 1, and using Bingham's theory to acquire the relationship of μ_∞ = μ_B, τ_B =  μ₀/K, Shao and Lo (2003) obtained a simplified Cross equation as: (17) where μ₀ = 1000μ_∞, K can be determined by their method and μ_∞ can be experimentally measured according to Bingham's theory. By certain approximations of the Cross equation, several other popular and rheology models can be obtained, such as the power law model, the Sisko model, and the Bingham's model (Barnes, Hutton, & Walters, 1989). However, to implement the Cross equation, the above-mentioned parameters need to be pre-determined before the implementation of the equation. In the current study, an experimental method is proposed to determine these parameters.

3.2. Parameter determination for the Cross equation

As indicated by Cross (1965), the flow behavior index n needs to be determined first. The following section introduces a hybrid method with four steps to determine the four required parameters in the Cross equation. To demonstrate the proposed method for determining the parameters for the Cross equation, water-kaolinite mixtures with a solid volume concentration of 20% are used.

To determine the four parameters, a measured viscosity and shear rate relationship (Figure 2) is required, which can be usually measured using rheometers. A shear stress and shear rate relationship can also be measured by the rheometers and used to check the accuracy of the parameter estimates.

Figure 2. Viscosity and shear rate relationship of water-kaolinite mixture at C_v = 20%.

The first step is to determine the flow behavior index n, which is determined at first by a variation of the graphical method proposed by Cross (1965). At a high shear rate αγⁿ >> 1, therefore the Cross equation can be reduced to: (18)

Converting Equation (18) into logarithmic form, the following relationship can be obtained: (19)

Since μ_∞ is the viscosity at a very high shear rate, it is extremely small compared to the effective viscosity. The material property changes dramatically; hence, the relationship of μ >> μ_∞ can be used, and μ_∞ on the left-hand side of Equation (19) can be ignored: (20)

Then, the relationship between log μ and log γ can be drawn from the measurement. It should be noted that log μ is linear with log γ (Figure 3), and n can be obtained as the slope of the linear relationship of log (u) and log (γ).

Figure 3. Linear relationship between log μ and log γ to determine the flow behavior index n.

The second step is to determine the viscosity at the very high shear rate μ_∞ by the graphical method. After the flow behavior index n is determined, μ_∞ can then be determined. Converting Equation (16) into the form (21) it can be seen that the relationship between μ and γ⁻ⁿ is linear with the interception of μ_∞, as shown in Figure 3. Therefore, μ_∞ can be determined (Figure 4).

Figure 4. Linear relationship between μ and γ-n to determine the viscosity at the very high shear rate μ∞.

The third step is to determine the viscosity at the very low shear rate μ₀ by the relationship of μ₀ = 1000μ_∞. Owing to the findings of Hammad and Vradis (1994) that μ₀ has nearly no effect on the numerical accuracy when μ₀ is 1000 times greater than μ_∞, a fixed value of μ₀ = 1000μ_∞ is selected and used.

The fourth step is to determine the consistency index α by regressing fitting according to the Cross equation using the values of μ₀, μ_∞, and n from steps 1 to 3. The unique value of α can be obtained as the best-fitting curve with the measurement.

The final obtained values for the four parameters of the water-kaolinite mixtures with solid volume concentration C_v of 20% are shown in Figure 5(a), along with the values for C_v = 22.5%, C_v = 17.5%, and C_v = 15%. From Figure 5, it can be seen that the viscosity and shear rate relationship using the estimated parameters is in reasonably agreement with the experimental results. In Figure 5(c), the Cross fittings have slight discrepancies between the shear rate range of 0 1/s and 15 1/s. Though the method to determine the parameters is based on the viscosity and shear rate relationships at high shear rate regions, the estimated parameters are nonetheless sufficiently accurate in depicting the non-Newtonian dam break flows, as indicated in the next section.

Figure 5. The relationship between the viscosity and the shear rate by Cross fitting using the proposed method for the water-kaolinite mixtures at different values of C_v.

4. Applications

In order to validate the proposed method for determining the parameters for the Cross equation and its performance in the WC-MPS model, three different groups of dam break tests with different materials are used.

Figure 6 shows a sketch of the dam break test simulated in the current study. The initial channel length L, initial channel height H, initial mixture column length l₀, and initial mixture column height h₀ vary for different groups of tests; the channel width varies as well, but it is not considered here since the simulations are 2D. Water-kaolinite mixtures with a solid volume concentration of 25% were used for the first group of simulations (Xie, Tai, & Jin, 2013), water-clay mixtures (Komatina & Jovanovic, 1997) with different solid volume concentrations were used for the second group, and Carbopol 940 aqueous solutions with different Carbopol volume concentrations (Minussi & Maciel, 2012) were used for the last group.

Figure 6. Sketch of the dam break flow problem.

4.1. Water-kaolinite mixtures of 25% solid volume concentration

The first group of experiments was carried out in a horizontal channel made of transparent Plexiglas of 1392 mm length, 51 mm width, and 200 mm depth. A high-speed CMOS camera was used to capture the flow movement. The camera was located 1.5 m away from the channel. The initial water-kaolinite column was set as 356 mm long and 100 mm high.

For the numerical simulations, 8900 fluid particles were used with a total number of 11,014 particles including the dummy particles. The particle size (or the initial particle distance) was set as 0.002 m. Before the dam break flow was initiated, all the particles were distributed evenly and the hydrostatic pressure conditions were assigned to the particles.

According to the experimental measurements, the time span for the removal of the confining gate was 0.1 s. The movement of the confining gate was considered during the simulation. The speed of removing the gate was assumed to be constant as 1 m/s, and consequently the gate was fully removed 0.1 s after the simulation was initiated.

Figure 7(a) shows the estimated Cross parameters and the corresponding shear stress and shear rate relationships for the 25% water-kaolinite mixture. From Figure 7(a) it can be seen that the viscosity and shear rate relationship using the estimated Cross parameters agree well with the measurements, with a slight underestimation of the viscosity at low shear rates. Moreover, Figure 7(b) shows the comparison of the shear stress and shear rate relationship between the estimated value and the measurements. The underestimation of the shear stress at lower shear rates is due to the underestimated viscosities, which further confirms the current method for determining the four parameters for the Cross equation. These estimated parameters are then used in the following simulation of the dam break flow of the 25% water-kaolinite mixture.

Figure 7. The relationship between (a) the viscosity and the shear rate and (b) the shear stress and the shear rate for water-kaolinite mixture at C_v = 25%.

The relative root mean square error (RRMSE) between the calculated viscosity using the proposed method and the measured ones is shown in Figures 5 and 7. The small RRMSE values indicate good agreement between the calculated viscosities and the experimental measurements. However, as shown in Figure 7(a), when the shear rate is less than 5 1/s, the measured values of viscosities in this shear rate range are not stable. This may be because the more viscous nature of the mixture at higher solid volume concentrations leads to difficulties in accurately measuring the viscosities at lower shear rates (0 1/s to 5 1/s).

Figure 8 indicates the flow regimes of both the laboratory measurements and the WC-MPS simulations for the water-kaolinite mixtures at a concentration of 25%. The simulations show reasonable agreement with the experiments. It can been seen from Figure 8 that, due to the more viscous nature of the water-kaolinite mixture which is different from the water dam break, the dam break of the water-kaolinite mixture has its flow depth and front propagating velocity decaying more rapidly over time. This important flow feature of a typical non-Newtonian dam break flow was observed in the simulations as well.

Figure 8. Surface profile comparison of the water-kaolinite mixture at C_v = 25%.

More specifically, since the confining gate was considered to be removed at a constant speed of 1 m/s and was opened completely at 0.1 s following the simulation initiation, the flow pattern was considerably affected before 0.1 s. There is a relatively slow fracturing avalanche of the initial reservoir heaps due to the partially opened confining gate before 0.1 s (Figure 8(a)), which transits to more violent cascading collapses of the taller part of the material column (Figure 8(b)). At 0.17 s after the complete gate removal, a rapid drop of the flow depth due to the influence of gravity was observed near the leading edge of the mixture column. The flow depth reduced to approximately 50 mm (half of the initial column height) at x = 0.356 m, after which the flow depths were almost fixed at around 50 mm with very little fluctuation. The numerical simulations of WC-MPS have shown this feature clearly during the flow (Figure 8).

On the other hand, at the beginning of the flow, right after the initial hydrostatic state (Figure 8(a)), the numerical simulation shows discrepancies with the experimental measurements. This is largely attributed to the high viscosity of the 25% water-kaolinite mixture. More specifically, owing to the highly viscous nature of the mixture and the high viscosity, the mixture will exhibit a more solid-like behavior due to being attached to the confining gate and moving upwards with the gate for a short distance. Since the gate removal in the experiments is not instantaneous, the dam break is modeled by removing gate particles with a constant speed of 1 m/s in the simulation process. The gate opening was also simulated by assigning no-slip boundary conditions to the gate modeling process. There are some improvements with no-slip boundary condition, as shown in the simulation results. Fu and Jin (2014) studied the effect of using a smaller particle size in their dam break cases. They indicated that by using a smaller particle size, a more accurate simulation was obtained but with the downside of requiring an extremely time-consuming calculation. As shown in Figure 8, the simulated water surface profiles are similar for the particle sizes d_s = 0.002 m and d_s = 0.001 m, even though the total number of fluid particles are four times less and the simulation time is at least eight times less.

A small overestimation of the flow depths was observed in the simulations. As mentioned above, due to the highly viscous nature of the mixture and the corresponding high viscosity, the mixture is more difficult to propel as it sticks to the gate, so the amount of material mass passing through the gate decreases during the removal of the confining gate, leading to the slightly overestimated flow depths.

4.2. Water-clay mixtures

Komatina and Jovanovic (1997) conducted experiments with water-clay mixtures of dam break flow that included 69 tests with different clay volume concentrations ranging between 0.0% and 36.1% for the initial reservoir height from 100 mm to 300 mm and channel slopes of 0.0% and 0.1%. Three cases, with an initial reservoir height of 100 mm, a channel slope of 0.1%, and clay volume concentrations of 13.9%, 20.1%, and 27.4%, respectively, were selected from their work for the numerical investigations in the current study. Table 1 illustrates the variation of the properties of the water-clay mixture with the clay volume concentration.

Table 1. Physical properties of the water-clay mixtures.

Fluid	Concentration Cv (%)	Density (km/m^3)	Yield stress τc (Pa)	Plastic viscosity µ(Pa·s)
1	13.9	1230	0.2	0.0069
2	20.1	1332	1.96	0.0201
3	27.4	1452	25.0	0.07

In the numerical simulations, a total of 59,632 particles 0.002 m in size were constructed, including 50,000 fluid particles. All particles were uniformly distributed initially and the hydrostatic pressure condition was assigned. As with the same treatment in the previous section, a moveable confining gate was set during the simulation and the gate was assumed to move upwards completely within 0.1 s when simulation was initiated.

Table 2 shows the parameters estimated for the Cross equation by the method proposed in this study for the water-clay mixtures, and Figure 9 shows the calculated shear stress; the shear rate relationships with the viscosities computed from the Cross equation in comparison to the Bingham theory are shown as well. The Bingham model curve in Figure 9 refers to Komatina and Jovanovic's (1997) findings for water-clay mixtures, which state that Bingham's model can be used to study the behaviors of the mixtures by neglecting the non-linearity of the viscosity at low shear rates. However, the non-linearity of the viscosity at low shear rates will affect the flow features of the mixtures at the beginning of the flows, as the viscosity changes dramatically during the early stage of the flows. This is the reason the Cross equation is used to perform the numerical tests for the dam break flow of the water-clay mixtures in the current study, since the dam break flows usually occur over a very short timescale and in the narrow shear rate region, during which the viscosity of the material decreases sharply. It can be seen from Figure 9 that the calculated relationships using Cross equations are in good agreement with the experimental results. The zoomed plots at the right-hand side of each plot show the comparisons at lower shear rates from 0 1/s to 10 1/s. It is obvious that at the lower shear rate, the Cross theory provides a non-linear relationship between the shear stress and the shear rate rather than the fixed value of the Bingham theory.

Figure 9. Shear stress and shear rate relationships for the water-clay mixtures for different values of C_v.

Note: The graphs on the right-hand side show the enlarged plots at a low shear rate from 0 1/s to 10 1/s.

Table 2. Estimated Cross parameters for the water-clay mixtures.

Fluid	Cv (%)	µ∞ (Pa·s)	µ0 (Pa·s)	α (Pa·s^n)	n
1	13.9	0.0069	6.9	32.0	0.988
2	20.1	0.0196	19.6	10.5	0.98
3	27.4	0.066	66	2.75	0.986

Figure 9 also shows the RRMSE values between the calculated shear stress from the proposed method and the calculations from Bingham's equation with the estimated Bingham parameters from experiments (Komatina & Jovanovic, 1997). The small RRMSE values also show good agreement between the calculated shear stress and the shear stress from Bingham's equation.

However, there are discrepancies between the two curves in Figure 9 at a very low shear rate (<5 1/s), and these discrepancies contribute to the increased RRMSE values even when the RRMSE is small. These discrepancies are due to the different descriptions of flows at a very low shear rate. Bingham's model can be considered as a simplified form of the Cross equation (Barnes et al., 1989), while Bingham's rheology law assumes no movement below the yield stress even when the material is undergoing small shears. Therefore, when a material is experiencing smaller shears than the Bingham's yield stresses, the actual deformation or movement cannot be reflected by using Bingham's model. As such, Bingham's law may not accurately describe the flow features at a very low shear rate. It is worth noting that the discrepancies become smaller as the shear rate increases (Figure 9(a) to (c)), indicating that both the Cross equation and Bingham's equation are comparable at higher shear rates.

The front depths at different cross-sections, namely, the stage hydrographs, are shown in Figure 10. It is clear that the MPS simulations take a longer time to reach the same flow depth of the measurements in all three cross-sections. The major reason for this latency is the particles reaching the cross-sections.

Figure 10. Stage hydrographs at different locations for the water-clay mixtures.

In Figure 10(a), x/H = 0.83 refers to the cross-section immediately after the initial dam site. The MPS simulation of the flow depth at this depth reaches the same value as the measurement at around 0.4 s for all concentrations, while further away from the initial dam site (Figure 10(b)), x/H = 2.50, the MPS simulations take 0.6 s to reach the measured depths for all concentrations. However, with a lower concentration, the simulation results of the flow depth take less time to reach the measured depths due to the larger viscosities of the higher concentrations and the existence of the confining gate; fluid of a higher viscosity is more prone to sticking to the gate in the real flows, causing a delay at the initiation of the flow – but this delay not included in the current MPS model.

At the furthest distance measured from the initial dam site (Figure 10(c)), x/H = 4.17, no flow reaches the location for the concentration 27.4% due to the slower movement caused by the high viscosity of the mixture. The MPS simulation shows better agreement here than with the previous two cross-sections.

Figure 11 shows the surface comparisons between the simulations using the method proposed and the experiments at different times, and the simulation results of the simplified Cross equation with WC-MPS (Xie et al., 2013) are shown as well. It is clear that the surface profiles of the current simulations are in good agreement with the experimental results. This is the advantage of the particle method in tracking free surfaces. In Figure 11(b), slight discrepancies are observed between x = 1.6 m and x = 2 m in the MPS simulations compared with the experiments. It is suspected that this is due to the local particles at these locations, as the weakly compressible nature of the current MPS model is for weakly compressible fluid that may lead to particles clustering at these locations. In the current model, particle clustering is avoided by using repulsive force in the pressure calculation (Koshizuka et al., 1998). Snapshots of the WC-MPS results at t = 0.27 s and t = 0.32 s for 25% simulation are shown in Figure 12. It can been seen that there are a few particles which slightly overlap with each other, but a more reasonable overall distribution of the fluid particle is obtained due to the existence of the repulsive forces.

Figure 11. Surface profiles for water-clay mixtures of C_v = 27.4% at two different times.

Figure 12. Snapshots of typical particles clustering in the simulations.

In order to further study the capabilities of the MPS method for simulating non-Newtonian flows with different concentrations, the wave front propagation was compared for different concentrations. Good agreement was found between the simulations and measurements (Figure 13). The speed of the wave front propagation decreases as the mixture concentration increases, and this viscous flow feature can be seen in Figure 13 as well. This also confirms the validity of the method proposed in the current study in determining the parameters for the Cross equation and the ability of the current WC-MPS model with the Cross equation in simulating non-Newtonian fluid flows.

Figure 13. Wave front propagation for the water-clay mixtures for different values of C_v.

4.3. Carbopol 940 aqueous solutions

The last group of experiments selected in this study is the work of Minussi and Maciel (2012). In their experimental study, rheological tests were performed for Carbopol 940 aqueous solutions of five different concentrations. Among the experimental studies, three groups of dam break tests were performed with different initial reservoir heights h₀. Two initial reservoir heights, h₀ = 100 mm and h₀ = 130 mm, were selected in the current study to validate the numerical model. Table 3 gives the details of the selected experiments. The density of the Carbopol 940 solutions was selected as 1007 kg/m³ as suggested by Minussi and Maciel (2012) and kept the same for all five fluids, since the concentration variations of the five fluids are small.

Table 3. Setup of the Carbopol 940 dam break tests.

Group 1			Group 2
Test	h0 ± 5 (mm)	Fluid	Test	h0 ± 5 (mm)	Fluid
1	100	1	2	130	1
4	100	2	5	130	2
7	100	3	8	130	3
9	100	4	10	130	4
12	100	5	13	130	5

Table 4 lists the properties of the five Carbopol 940 aqueous solutions for different concentrations.

Table 4. Properties of the Carbopol 940 solutions and the corresponding rheological parameters for the Herchel–Bulkley model.

Fluid	Cv (%) ± 10^−4	t ± 0.1(°C)	τc (Pa)	k (Pa·s^n)	θ
1	0.116	28	30.002	4.297	0.479
2	0.099	28	18.242	1.904	0.531
3	0.139	30	49.179	7.837	0.442
4	0.112	30	31.366	4.550	0.470
5	0.139	18	39.170	7.263	0.446

Note: The last three columns show the rheological properties of the fluids for the Herchel–Bulkley model, where τ_c is the yield stress, k is the consistency, and θ is the flow index.

During the numerical simulations, the particle size for both groups was set as 0.002 m. There were 12,500 fluid particles and a total of 16,954 particles (including dummy particles) for the tests in Group 1 and 16,250 fluid particles and a total of 20,704 particles for the tests in Group 2. All particles were uniformly distributed initially and the hydrostatic pressure condition was assigned to the particles.

Figure 14 shows the comparisons between the Cross equation using the fitted parameters and the calculation by the Herchel–Bulkley model with the estimated parameters (Minussi & Maciel, 2012). The RRMSE values are shown as well and indicate reasonable agreement between the calculated viscosities and the viscosities from Minussi and Maciel (2012). There are discrepancies between the calculated viscosities and the viscosities from Minussi and Maciel (2012) at the shear rate range of 10 1/s and 50 1/s. Since the Cross equation and the Herchel–Bulkley model provide different descriptions at a low shear rate, the discrepancies for these two rheology laws become smaller as the shear rate increases.

Figure 14. The relationship between the viscosity and the shear rate for the Carbopol 940 aqueous solutions with different constituents.

Therefore, as shown in the figures, the relationships estimated from the proposed method show reasonable agreement with the apparent viscosity obtained by Minussi and Maciel (2012), with slight underestimations of the viscosity between the shear rate of 25 1/s and 50 1/s. On the other hand, the estimated relationship shows good agreement with the Herchel–Bulkley model, confirming the validity of the method in describing the rheological properties for non-Newtonian materials.

Table 5 shows the estimated Cross parameters for the Carbopol 940 aqueous solutions using the proposed method, and Figures 15 and 16 illustrate wave front comparisons between the WC-MPS results and the experiments of Minussi and Maciel (2012) for Group 1 and Group 2, respectively. The numerical results of Minussi and Maciel (2012) using ANSYS CFX software are shown as well.

Figure 15. Wave front locations for the Group 1 tests.

Figure 16. Wave front locations for the Group 2 tests.

Table 5. Estimated Cross parameters for the Carbopol 940 solutions.

Fluid	Cv (%) ± 10^−4	µ∞ (Pa·s)	µ0 (Pa·s)	α (Pa·s^n)	n
1	0.116	0.1037	103.7	2.80	0.9
2	0.099	0.0781	78.1	3.60	0.92
3	0.139	0.1786	178.6	3.02	0.91
4	0.112	0.1295	129.5	3.68	0.9
5	0.139	0.1674	167.4	3.45	0.9

For both of the test groups (100 mm and 130 mm initial mixture height for Group 1 and Group 2, respectively), the WC-MPS simulations underestimated the wave front locations compared with the measurements. The major reason for this is the errors that are introduced in the parameter fitting process for the Cross equation and the errors induced by the use of the Cross equation, since the Cross equation itself is a constitutive equation derived based on the assumption of the pseudoplastic flow of non-Newtonian materials in connection with the formation and rupture of structural linkages.

From Figures 15 and 16, it can be seen that the current simulations show overestimations of the wave front locations at the early stage of the flows for several tests, i.e., test 1 (Figure 15(a)), test 2 (Figure 16(a)), test 4 (Figure 15(b)), test 7 (Figure 15(c)), and test 8 (Figure 16(c). By analyzing the relationships of the materials in Figure 14, the µ₀ (viscosity at a very low shear rate) from the Cross fitting is larger than the apparent viscosity, and the estimated viscosity decreases sharply with the increase of the shear rate ranging through 0 1/s to 50 1/s. Then, the estimated viscosity agrees well with the apparent viscosity. Since the WC-MPS method adopts the estimations of the Cross equation and the higher µ₀, the material is slower to initiate flow at lower shear rates (i.e., at the beginning of the collapse). Thus, the sharp decrease in viscosity with the shear rate indicates that the material has a lower viscosity under the same shear rate compared with the experimental results. Subsequently, the lower viscosity induces a larger front propagating velocity and, added to the absence of the confining gate in the simulations, this finally leads to the overestimation of the wave front locations. In addition, the discrepancy between the numerical and experimental results may due to the application of the WC-MPS method, and modifications such as using an enhanced MPS technique (Khayyer & Gotoh, 2013) may improve the fit between the numerical and experimental results.

Nonetheless, the model still underestimated the wave front locations. From the relationship between the viscosity and the shear rate fitted from the Cross equation, the viscosity drops sharply at first and then maintains a smooth and steady decrease over the shear rate change between 0 1/s and 50 1/s, which indicates that the model overestimates the propagating velocity at first, and then the velocity decreases at a very slow pace while the apparent viscosity drops dramatically within the same shear rate range, and this sharp decrease in viscosity, at last, leads to a faster increase in propagating velocity than the model predicts. Consequently, the model has an underestimated velocity and thus underestimates the wave front locations.

5. Conclusion

In this study, the WC-MPS model was used to study the free surface flow of non-Newtonian fluid. The flow development process and regimes are directly represented by particle interactions and movements. The non-Newtonian nature and flow behaviors are accounted for by the Cross rheology equation coupled with an experiment-based method to estimate the parameters required by the Cross equation. The viscosity and the shear rate relationships obtained from the Cross equation and the proposed parameter-determination method were compared with the experimental measurements for three different sets of tests with different non-Newtonian materials. Good agreement was found between the proposed method and the measurements. The estimated parameters together with the Cross equation were then applied within the WC-MPS model to numerically study the dam break flows of the materials. The surface profiles and wave front locations of the simulations are in good agreement with the experimental results, indicating that the proposed method to estimate the Cross parameters and the use of the Cross equation can capture the flow regimes of non-Newtonian materials, as well as showing the capability of the weakly compressible particle method to simulate non-Newtonian free surface flows.

However, the current WC-MPS model has several limitations. Comparisons of the stage hydrographs of the water-clay mixtures show that in order to accurately track the front depth of the flows, more particles with decreased particle sizes are required, which will increase the computation cost. A more advanced computing method and technology to accompany it will greatly improve the efficiency when large amounts of particles are employed. On the other hand, a more sophisticated interaction model is absent from the current model to further improve the accuracy, including the physical mechanisms between the fluid and the confining gate, as well as the interaction between the non-Newtonian material and the solid boundaries. The current model is limited to 2D simulations, thus a 3D model considering the effects of the lateral directions along the channel width will improve the applicability of the WC-MPS method to more practical engineering problems with non-Newtonian materials.

Nonetheless, due to the advantage of the WC-MPS method in handling flows with large deformation and fragmentation, and its capability of treating different types of fluid flows, the WC-MPS method can be expanded to study more practical problems, such as sophisticated 3D flows (C. L. Wu & Chau, 2006), salinity and temperature (Chau & Jiang, 2001, 2004), and multiphase flow problems (T. Liu & Yang, 2014). Moreover, using the MPS method to solve the shallow water equation (SWE; Lai & Khan, 2012; Pu et al., 2013) also hold great promise in improving the devising of solutions for practical problems as well.

Acknowledgements

The authors would like to thank Dr Y. C. Tai and Mr Bo-Jun Chen for providing experimental data.

Disclosure statement

No potential conflict of interest was reported by the authors.

Additional information

Funding

This work was supported by the Natural Sciences and Engineering Research Council of Canada [grant number 429672-2012].