Dependence of Sediment Suspension Viscosity on Solid Concentration : A Simple General Equation

In this study, a simple and new parametric equation is proposed that describes the relative viscosity of a suspension as a function of suspended solid concentration, covering a range from very dilute to highly concentrated states, based on Costa (2005). The proposed viscosity law depends mainly on the solid volumetric concentration φ and maximum packing fraction φm at which a regime transition occurs. Furthermore, the viscosity of dilute as well as a highly-concentrated kaolinite suspension (in excess of 40% and lower than 50%) was measured in a coaxial cylinder rheometer. The proposed formula shows a capability of fitting the measured bulk viscosity into the hyper-concentrated regime in the experiment. Finally, the proposed equation could also be found to show a good fitting relationship with published experimental results on partially crystallised Li2Si2O5 melt and partially melted granite with solid contents from 75% to 95%.


Introduction
Cohesive sediment commonly exists in harbour, estuary and coastal waters all over the world [1].Understanding cohesive sediment transport is important for engineering projects related to dredging operations of navigation channels, morphological changes of shorelines and water quality control in estuarine and coastal areas [2][3][4][5][6].Transport of cohesive sediment occurs both in the water column and in the near-bed water layer.For the latter, the sediment concentration is often very high and the term "fluid mud" has been used to describe it [7][8][9].Motion of fluid mud layer, which is caused by water current and wave, may lead to a high sediment transport rate [6].In order to predict the onset and motion characteristics of the fluid mud layer, it may be important to understand their rheological properties (including bulk viscosity and yield stress) [6,10,11].Furthermore, the viscosity of the fluid mud layer has been shown to be important for evaluation of the entrainment rate of sediment between the fluid mud layer and the overlying water layer [12,13].Meanwhile, the hindering settling velocities of cohesive sediment flocs, which is an important variable for calculation of suspended sediment flux, are also related to the viscosity of water-sediment suspension when sediment concentration is high enough [14,15].Thus, investigating rheological properties of the cohesive sediment mixture especially in the case of a highly concentrated suspension, where no experimental data are seemingly available [2,3,14], could be helpful to better understand cohesive sediment transport in estuarine and coastal waters.This study focuses on the suspension viscosity.The suspension viscosity, not limited to the sediment mixture, is affected by particle concentration, shear rate, particle shape, particle size distribution and particle deformability [16][17][18].Among these factors, the dependence of suspension viscosity on the particle concentration is of great interest in the sediment field [2,3,19,20].Experiments have been performed to study this relationship, and some models have been used to describe it [21][22][23][24][25][26].However, a simple relationship may not exist because the particle in a suspending medium can Water 2017, 9, 474 2 of 14 be subject to hydrodynamic, Brownian, colloidal force and other effects, which are not easily modelled in a simple form [16,[27][28][29].
Some proposed models characterizing the suspension viscosity-solid particle concentration relationship are simply summarized as follows.In the case of a notably dilute suspension of spherical non-interacting mono-disperse particles, Einstein [30] derived an analytical solution to describe the relative viscosity of the suspension (η r ) (defined as the ratio of the suspension viscosity and viscosity of the suspended medium) as a function of solid fraction, where φ is the volumetric concentration (or fraction) of particles, defined as the volume occupied by particles and total volume of the suspension, and B is the Einstein coefficient with a theoretical value of 2.5 for hard spheres.Batchelor [31] theoretically extended this relationship to the second order, where k = 6.2 for Brownian suspensions in any flow, and k = 7.6 for non-Brownian suspensions in a straining flow.When particle concentrations are high, particle crowding causes hydrodynamic interactions among particles, which results in large positive deviations from the viscosity value of highly concentrated suspensions from Equations ( 1) and ( 2) [16].
For the case of a highly concentrated suspension, Roscoe [32] proposed the following expression starting from the Einstein equation, where the new parameter φ m is the maximum packing concentration of particles, at which there is no longer sufficient fluid to lubricate relative particle motion, and the incompressible solid would prevent any flow [21].Its value varies depending on various factors, such as the variation of particles from mono-size, and the geometric arrangement of particles [16][17][18].Equations ( 1)-(3) apply only to the suspension in which non-cohesive particles are suspended.In such cases, Boyer et al. [25] proposed the following φ-dependent constitutive laws of dense suspensions, by modelling the friction law of dense suspension as the sum of two contributions coming from contact between particles and hydrodynamic stresses between particle and surrounding fluid, an µ 1 = 0.32, µ 2 = 0.7 and I 0 = 0.005 were suggested in Boyer et al. [25].Later, by conducting numerical simulation of an athermal dense suspension under shear, Trulsson et al. [26] showed that the Newtonian and Bagnoldian regimes of dense suspension rheology could be unified as a function of a single dimensionless number.
On the other hand, some rheologists have proposed many semi-empirical or empirical models to fit the observed experimental results, considering Equation (1) as a starting point [21][22][23][24].A frequently adopted model that fits experimental data well is the semi-empirical expression, introduced by Krieger and Dougherty [21] for mono-dispersed suspensions, For very dilute suspensions ( φ → 0 ), this expression can be reduced to the Einstein equation: η r = 1 + Bφ.Some experiments have shown that the exponent of Equation (5) (i.e., Bφ m ) is always approximately two for many suspensions [33,34].The exact particle concentration (φ m ) at which there Water 2017, 9, 474 3 of 14 is no fluid to lubricate particle motion for a system of particles remains a matter of debate [16,34], and its determination appears to be an open problem.In many experiments, φ m was commonly used as an adjustable parameter, so that Equation (5) can provide an excellent fitting relationship with the experimental data.Other fitting expressions proposed by some authors also used φ m as an adjustable parameter, like for example, the equation proposed by Chong et al. [22], and the proposed equation by Dabak and Yucel [23], After summarizing the common form of some semi-empirical fitting expressions, Liu [24] proposed a more general formula to describe the viscosity behaviour of suspensions, with φ m = (1 − q)/p, where p and q are the slope and intercept values of the straight line 1 − η −1/2 r vs. φ relation, which must be experimentally determined, and n is a flow-dependent suspension-specific parameter (the author used n = 2).The term (φ m − φ) is clearly interpreted as the effective space available for particles to move in the matrix fluid.
The majority of the proposed models for a highly concentrated suspension assume that there is a critical solid concentration φ m at which the incompressible solid would prevent any flow.It may be reasonable to consider that when the solid concentration increases across φ m , the suspension experiences a transition of rheological regime: the suspension rheology changes from being primarily determined by the liquid phase to exhibiting remarkably non-Newtonian behaviour with a yield stress [19,20,28,29,35].In the latter regime, the rheology is much more complex than the former regime.Costa noticed this transition of rheological regime of the suspension, and further proposed a general parameterization to characterize the suspension viscosity-solid concentration relationship as follows, where α, β and γ are three adjustable parameters.According to Equation ( 9), the suspension viscosity initially increases (here we term it an initial linear increasing regime), and then experiences an intermediate increasing process (here called an intermediate increasing regime) with the increasing solid particle concentration.However, with the continual increase of solid concentration, the suspension viscosity predicted by Equation ( 9) exhibits an asymptotic behaviour and finally reaches a steady state for a very high solid concentration.The assumption that the highly-concentrated suspension has a fixed viscosity value seems questionable.In fact, our experimental works (as discussed in Section 4) on the viscosity of kaolinite suspension across a certain range in particle concentration shows that the suspension viscosity increases with increasing solid concentration.Moreover, measured data for viscosity-crystal fraction relationship in partially crystallised Li 2 Si 2 O 5 and observation data of partially melted granite with solid content from 75% to 95% also clearly indicated this increasing behaviour when crystal fraction exceeds the maximum packing concentration (as also discussed in Section 4) [36,37].Costa [35] has used these data to verify Equation ( 9), but neglecting the increasing behaviour of the suspension viscosity at a high solid concentration.Therefore, a more general η r − φ parameterization may be necessary to describe the complex behaviour of the suspension, covering a range of solid concentration from very dilute to highly concentrated states.This study attempts to identify a more general relationship between η r and φ, based on Costa [35].The proposed equation is intended to not only describe the aforementioned rheological transitions of the suspension (especially the increasing trend at high solid concentration), but also to provide adjustable parameters that can be easily interpreted from a physical perspective and that reduce to the Einstein equation for a very dilute suspension ( φ → 0 ), similar to previously proposed models.This paper is arranged as follows.Section 2 proposes a simple new η r − φ parameterization.An experiment that measures viscosity-solid concentration values of a water-kaolinite suspension, performed to verify the validity of the proposed equation, is detailed in Section 3. Section 4 presents the comparisons of the proposed equation with experimental results in this study and previously published experimental results on partially crystallised Li 2 Si 2 O 5 and partially melted granite.Finally, some concluding remarks are presented in Section 5.

Simple Newly-Proposed Equation
In order to reproduce the continually increasing behaviour of η r for φ higher than φ m , we simply add a power function of solid concentration φ (up to second order) on the right-handed side of Costa equation, so that the final form of the suspension viscosity η r as a function of solid concentration φ becomes where δ is an exponent (δ > 1), and φ/φ m is the non-dimensional form of the solid concentration.The error function is the integral of the Gaussian distribution: x 0 exp(−m 2 )dm at m ≥ 0, and it shows a rapid and an asymptotic increasing behaviour with continually increasing x.Its Taylor expansion near x = 0 is: For very dilute suspensions ( φ → 0 ), the argument of the error function approaches Meanwhile, the square bracket part of the right-handed side of Equation (10) approaches unity when φ → 0 .Therefore, the η r −φ parameterisation for very dilute suspensions becomes: η r (φ) ≈ 1 + (1 + β)Bφ, which is similar to well-established Einstein equation valid for dilute suspensions, where the Einstein coefficient becomes (1 + β)B, and B could be regarded as a correction coefficient due to the non-linear term in the argument of the error function.For highly concentrated solid-liquid systems (large φ), the error function approaches unity with increasing φ, and the curly bracket part of the right-handed side of Equation ( 10) becomes (1 − α) −B/α .Moreover, the square bracket part exhibits an increasing behaviour with φ.As a result, η r for highly concentrated solid-liquid systems shows an expected increasing trend with increasing φ.
A graphical analysis is employed to explore the physical explanations of all parameters incorporated in Equation ( 10), as shown in Figure 1.As discussed previously, and shown in Figure 1a, the parameter α controls the viscosity value from which the suspension exhibits remarkably non-Newtonian behaviour with a yield stress, and the continually increasing regime of η r with φ for large φ begins.Parameter β (Figure 1b), controls the slope of the initial linear regime of viscosity (as shown in reduced form of Equation ( 10) as φ → 0 ), and partially determines the slope curve between the initial linear and continually increasing regime.Parameter γ (Figure 1c), similar to parameter β, determines the width of the intermediate increasing regime and can be regarded as a measure of the rapidity to reach the continually increasing regime from the initial linear increasing regime.A smaller parameter γ leads to a larger width of the intermediate increasing regime.The newly-added parameter δ in this study primarily controls the rate at which the viscosity increases with particle concentration in the continually increasing regime, as shown in Figure 1d.
to parameter  , determines the width of the intermediate increasing regime and can be regarded as a measure of the rapidity to reach the continually increasing regime from the initial linear increasing regime.A smaller parameter  leads to a larger width of the intermediate increasing regime.The newly-added parameter  in this study primarily controls the rate at which the viscosity increases with particle concentration in the continually increasing regime, as shown in Figure 1d.

Experiment Introduction
In order to avoid discussing the effect of particle mineral composition on the suspension [38], this study simply adopted kaolinite (China clay) suspended in deionized water to form the solid-liquid suspension.The grain size distribution of the kaolinite sample was measured using a laser diffraction particle size analyser (Helos & Rodos; produced by Sympatec Corporation, Clausthal-Zellerfeld, Germany), and more information about measurement method of this instrument (including dry dispersion for sample powder) can be found on the website [39].As shown in Figure 2, the median size of this sample was 3.71 µm with a range of 0.5-18.5 µm.The kaolinite concentration in the suspension was chosen to be in the range of 1.88 × 10 −3 to 0.48 in this study, covering a very dilute to highly concentrated state [2][3][4].The upper limit value of kaolinite concentration (i.e., 0.48) is determined based on the experimental finding that the suspension shear stress value, when particle concentrations are higher than 0.48, exceeds the upper limit value of measuring scope of the rheometer in this experiment.
size of this sample was 3.71 μm with a range of 0.5-18.5 μm.The kaolinite concentration in the suspension was chosen to be in the range of 1.88 × 10 −3 to 0.48 in this study, covering a very dilute to highly concentrated state [2][3][4].The upper limit value of kaolinite concentration (i.e., 0.48) is determined based on the experimental finding that the suspension shear stress value, when particle concentrations are higher than 0.48, exceeds the upper limit value of measuring scope of the rheometer in this experiment.Rheological properties of the water-kaolinite suspension were measured using a Brookfield RST-CCT40 rheometer(Brookfield Engineering Laboratories Incorporation, Middleboro, MA, USA) This rheometer mainly consists of a touch screen, an electronic unit and a measuring drive, which is coupled to coaxial cylinder systems and the Rheo3000 software (Rheo3000 1.2.1395.1,Brookfield Engineering Laboratories Incorporation, Middleboro, MA, USA.When a coaxial cylinder spindle, which consists of a shaft and a measuring bob, is immersed in the sample cup (filled with the waterkaolinite sample, volume 68.5 mL) and driven at different angular velocities, the torques experienced by the spindle are detected.The angular velocity of the spindle is converted to shear rate in the suspension using a shear rate factor, and the torque corresponds to shear stress using the shear stress factor.A comprehensive introduction to the rheometer configuration and measuring system can be found on the website [40].The Rheo 3000 software coupled with a personal computer was used to remotely manipulate the rheometer.The shear rate was set to be 2.16-644.41 1 s  , covering the low to high shear conditions in this experiment [3].
After the water-kaolinite suspension was prepared for a kaolinite concentration of  , it was gently injected into the sample cup, and the rheological measurement was initiated.The angular velocity of the spindle was adjusted to gradually increase over 100 s, resulting in a shear rate range of 2.16-644.41 1 s  The torque (corresponding to the shear stress) that the spindle experienced was measured.The measuring system plotted the shear stress-shear rate relationship in real time, and presented the suspension viscosity value.For every kaolinite concentration, we repeated the procedure, from preparing the sample to measuring rheological characteristics, three times.The viscosity value of the water-kaolinite suspension was finally determined as the mean of the three measured viscosity values.For some water-kaolinite suspensions, we carefully monitored the solid Rheological properties of the water-kaolinite suspension were measured using a Brookfield RST-CCT40 rheometer(Brookfield Engineering Laboratories Incorporation, Middleboro, MA, USA) This rheometer mainly consists of a touch screen, an electronic unit and a measuring drive, which is coupled to coaxial cylinder systems and the Rheo3000 software (Rheo3000 1.2.1395.1,Brookfield Engineering Laboratories Incorporation, Middleboro, MA, USA.When a coaxial cylinder spindle, which consists of a shaft and a measuring bob, is immersed in the sample cup (filled with the water-kaolinite sample, volume 68.5 mL) and driven at different angular velocities, the torques experienced by the spindle are detected.The angular velocity of the spindle is converted to shear rate in the suspension using a shear rate factor, and the torque corresponds to shear stress using the shear stress factor.A comprehensive introduction to the rheometer configuration and measuring system can be found on the website [40].The Rheo 3000 software coupled with a personal computer was used to remotely manipulate the rheometer.The shear rate was set to be 2.16-644.41s −1 , covering the low to high shear conditions in this experiment [3].
After the water-kaolinite suspension was prepared for a kaolinite concentration of φ, it was gently injected into the sample cup, and the rheological measurement was initiated.The angular velocity of the spindle was adjusted to gradually increase over 100 s, resulting in a shear rate range of 2.16-644.41s −1 The torque (corresponding to the shear stress) that the spindle experienced was measured.The measuring system plotted the shear stress-shear rate relationship in real time, and presented the suspension viscosity value.For every kaolinite concentration, we repeated the procedure, from preparing the sample to measuring rheological characteristics, three times.The viscosity value of the water-kaolinite suspension was finally determined as the mean of the three measured viscosity values.For some water-kaolinite suspensions, we carefully monitored the solid particle concentration before the experiment was initiated and again after the experiment was completed, and found that there was no obvious variation for particle concentration due to the possible appearance of the sedimentation effect of the kaolinite particle [19].Before all experiments, the viscosity of deionized water was also measured.The experimental temperature was 21 ± 0.5 • C.

Comparison with Experimental Results
Figure 3 shows the measured relative viscosity values of the water-kaolinite suspension for particle concentrations from 1.88 × 10 −3 to 0.48.With increasing φ, η r first increased linearly (0 < φ < 0.15) (as shown by blue background colour in Figure 3), then showed intermediate increasing behaviour (0.15 < φ < 0.42) (as shown by the green background colour), and finally a continuous increasing trend (0.42 < φ < 0.48) (as shown by the red background colour).The presence of the continually increasing regime during this experiment may indicate complicated rheological properties of the highly concentrated suspension, which exhibited obvious non-Newtonian characteristics and behaved as a solid with applied stress.We attempt to fit these measured data using the proposed equation and the formula already proposed by Costa [35], as well as published formulae proposed by other authors, as shown in Figure 4.


, where m and o are the modelled and observed points, respectively, and N is the number of observed points.Goodness of fit increases as 2 R increases and NRMSE decreases.As presented in Table 1, in general, the proposed equation provides the best fit with the measured data, because all four adjustable parameters have been incorporated, and Costa formula also exhibits a better fit result.They are followed by the Krieger and Dougherty formula [21], which gives a high 2 R of 0.73; the Boyer formula [25] provides a small NRMSE value equal to 0.14.As observed in Figure 4, with increasing  , all formulae predict the viscosity value well for 0 <  < 0.05; when  reaches 0.1, the Roscoe formula [32] and Liu formula [24] clearly deviate from the measured viscosity values.As  approaches 0.15, the Chong formula [22] begins to overestimate the viscosity value, and then underestimates it from  = 0.33.For We attempt to fit these measured data using the proposed equation and the formula already proposed by Costa [35], as well as published formulae proposed by other authors, as shown in Figure 4. Goodness of fit is measured with the correlation coefficient, R 2 , defined as R 2 = [cov(m, o)/σ m σ o ] 2 , and normalized root mean square error, NRMSE, defined as , where m and o are the modelled and observed points, respectively, and N is the number of observed points.Goodness of fit increases as R 2 increases and NRMSE decreases.As presented in Table 1, in general, the proposed equation provides the best fit with the measured data, because all four adjustable parameters have been incorporated, and Costa formula also exhibits a better fit result.They are followed by the Krieger and Dougherty formula [21], which gives a high R 2 of 0.73; the Boyer formula [25] provides a small NRMSE value equal to 0.14.As observed in Figure 4, with increasing φ, all formulae predict the viscosity value well for 0 < φ < 0.05; when φ reaches 0.1, the Roscoe formula [32] and Liu formula [24] clearly deviate from the measured viscosity values.As φ approaches 0.15, the Chong formula [22] begins to overestimate the viscosity value, and then underestimates it from φ = 0.33.For 0.2 < φ < 0.25, the Dabak [23] and Krieger and Dougherty [21] formulae slightly overestimate the suspension viscosity.The Dabak formula [23] begins to underestimate the viscosity from φ = 0.3 on, and from φ = 0.25 on, the Krieger and Dougherty formula [21] starts to underestimate the viscosity.For 0 < φ < 0.25, the Batchelor formula [31] still provides a better fit effect to measured data, whereas the Einstein equation [30] still underestimates the suspension viscosity value from φ = 0.1 on.Across the entire range of solid concentration, the Boyer formula [25] is still very close to the Chong formula [22] except for a narrow range of 0.45 < φ < 0.48, during which the viscosity values predicted by the Chong formula [22] are slightly larger than those by the Boyer formula [25].
The Costa formula [35] that contains three adjustable parameters can fit the initial increasing regime and intermediate increasing regime of suspension viscosity value as a function of solid concentration well.However, as solid concentration exceeds 0.42 (the highly-concentrated suspension starts to exhibit obvious non-Newtonian characteristics during the experiment), it simply tackles the suspension viscosity-solid concentration relationship as a steady state.In fact, when 0.42 < φ < 0.48, note that the vertical axis is logarithmic coordinate in Figure 4, the suspension viscosity η r rises rapidly with increasing solid concentration φ.The proposed formula is capable of modelling this continual increasing regime of suspension viscosity with solid concentration with the cost of adding another new adjustable parameter.
Chong formula [22] are slightly larger than those by the Boyer formula [25].
The Costa formula [35] that contains three adjustable parameters can fit the initial increasing regime and intermediate increasing regime of suspension viscosity value as a function of solid concentration well.However, as solid concentration exceeds 0.42 (the highly-concentrated suspension starts to exhibit obvious non-Newtonian characteristics during the experiment), it simply tackles the suspension viscosity-solid concentration relationship as a steady state.In fact, when 0.42 <  < 0.48, note that the vertical axis is logarithmic coordinate in Figure 4, the suspension viscosity r  rises rapidly with increasing solid concentration  .The proposed formula is capable of modelling this continual increasing regime of suspension viscosity with solid concentration with the cost of adding another new adjustable parameter.[30], the Batchelor formula [31], the Roscoe formula [32], the Krieger and Dougherty formula [21], the Chong formula [22], the Dabak formula [23], the Liu formula [24] and the Boyer formula [25].m  is taken as 0.48 in this study.Fitting values are:  = 0.8,  = 0.0096,  = 10.3, and  = 3.Table 1.Goodness of fit for the proposed equation, Costa formula and other published formulae (compared with measured data in Figure 3)., the Batchelor formula [31], the Roscoe formula [32], the Krieger and Dougherty formula [21], the Chong formula [22], the Dabak formula [23], the Liu formula [24] and the Boyer formula [25].φ m is taken as 0.48 in this study.Fitting values are: α = 0.8, β = 0.0096, γ = 10.3, and δ = 3.
Table 1.Goodness of fit for the proposed equation, Costa formula and other published formulae (compared with measured data in Figure 3).

Comparison with Published Observational Data
To the best of our knowledge, few measured viscosity values for water-sediment suspensions at high solid concentrations are available in the sediment research field, and it is possible that the rheometer does not work properly at high volumetric concentrations because of the critical yield stress, particle migration and non-homogeneous suspension [16,41].To further confirm the validity of the proposed equation, we attempt to use two groups of previously-published measurement data for high solid concentrations that were introduced by Costa [35].
The first group of published data was conducted by Lejeune and Richet [37].They measured the viscosity of partially crystallized Li 2 Si 2 O 5 under uniaxial compression as a function of crystal volumetric concentration.The Li 2 Si 2 O 5 melt contains small ellipsoidal crystals.Their experimental results suggested that the rheology of the melts is primarily determined by the crystal fraction, even though additional measurements would be helpful to determine the possible influence of other factors, such as the size distribution or shape of inclusions.
Figure 5 shows the measured viscosity-crystal fraction relationship for a partially crystallized Li 2 Si 2 O 5 melt.A comparison between measured data and predicted values by the proposed equation and that between measured data and predicted values by the Costa formula are also presented in this figure.The agreement with the proposed formula is fairly good, with an R 2 value of 0.99 and an NRMSE value of 0.03, whereas the Costa formula also gives a good fit effect with an R 2 value of 0.92 and an NRMSE value of 0.16.However, the Costa formula neglects the increasing behaviour of the melt viscosity as the crystal volumetric concentration increases from 0.47 to 0.67, and simply regards it as a plateau of viscosity with solid concentration.In contrast, the proposed formula captures the increasing characteristic of the melt viscosity-crystal concentration relationship.Finally, it should be pointed out that there are very limited experimental data points in this figure to calibrate the proposed equation; in particular, there are no available measured results for φ larger than 0.67.More experimental results at high solid concentrations to calibrate the proposed equation are necessary in future studies.
The first group of published data was conducted by Lejeune and Richet [37].They measured the viscosity of partially crystallized Li2Si2O5 under uniaxial compression as a function of crystal volumetric concentration.The Li2Si2O5 melt contains small ellipsoidal crystals.Their experimental results suggested that the rheology of the melts is primarily determined by the crystal fraction, even though additional measurements would be helpful to determine the possible influence of other factors, such as the size distribution or shape of inclusions.
Figure 5 shows the measured viscosity-crystal fraction relationship for a partially crystallized Li2Si2O5 melt.A comparison between measured data and predicted values by the proposed equation and that between measured data and predicted values by the Costa formula are also presented in this figure.The agreement with the proposed formula is fairly good, with an 2 R value of 0.99 and an NRMSE value of 0.03, whereas the Costa formula also gives a good fit effect with an 2 R value of 0.92 and an NRMSE value of 0.16.However, the Costa formula neglects the increasing behaviour of the melt viscosity as the crystal volumetric concentration increases from 0.47 to 0.67, and simply regards it as a plateau of viscosity with solid concentration.In contrast, the proposed formula captures the increasing characteristic of the melt viscosity-crystal concentration relationship.Finally, it should be pointed out that there are very limited experimental data points in this figure to calibrate the proposed equation; in particular, there are no available measured results for  larger than 0.67.
More experimental results at high solid concentrations to calibrate the proposed equation are necessary in future studies.We also compare these measured data with some published formulae proposed by other authors, as shown in Figure 6, and Table 2 shows goodness values of fitting for these published formulae.Since the Roscoe formula, the Krieger and Dougherty formula, the Chong formula, the We also compare these measured data with some published formulae proposed by other authors, as shown in Figure 6, and Table 2 shows goodness values of fitting for these published formulae.Since the Roscoe formula, the Krieger and Dougherty formula, the Chong formula, the Dabak formula, the Liu formula and the Boyer formula deviate greatly from the measured viscosity value at φ = 0.67, in general the Batchelor formula provides the best fit to measured data with an R 2 value of 0.94 and an NRMSE value of 0.45 among such formulae.For 0 < φ < 0.31, the Liu formula fits the measured data very well, followed by the Roscoe formula, whereas the rest formulae still underestimate the measured viscosity values.
Water 2017, 9, 474 10 of 14 Dabak formula, the Liu formula and the Boyer formula deviate greatly from the measured viscosity value at  = 0.67, in general the Batchelor formula provides the best fit to measured data with an 2 R value of 0.94 and an NRMSE value of 0.45 among such formulae.For 0 <  < 0.31, the Liu formula fits the measured data very well, followed by the Roscoe formula, whereas the rest formulae still underestimate the measured viscosity values.[37] with some published formulae, including the Einstein equation [30], the Batchelor formula [31], the Roscoe formula [32], the Krieger and Dougherty formula [21], the Chong formula [22], the Dabak formula [23], the Liu formula [24] and the Boyer formula [25].Here m  is taken as 0.67 in this study.
Table 2. Goodness of fit for some published formulae (compared with published measured data in Figure 5).The second group of published data was mainly conducted by van der Molen and Paterson [36].They investigated the deformation of partially melted granite with solid concentrations from 75% to 95% at 800 and 300 MPa confining pressure and measured viscosity-solid fraction values for partially melted granite.Their measurement results suggested that the critical crystal concentration differentiating granular-framework-controlled flow behaviour from the suspension-like behaviour could be around 0.72.Measured viscosity-crystal fraction values at low solid concentrations for partially crystallized Mg3Al2Si3O12 melt by Lejeune and Richet [37] (for the Mg3Al2Si3O12 melt, crystals are well rounded), and other viscosity-solid fraction values for uniform spherical particles by Thomas [42] were also incorporated into this group of published data in order to make a simple comparison between the available experimental results and the proposed equation [35], although these experimental data are limited and their sources are not coherent.
Figure 7 shows the comparisons of these measurement data with the proposed formula in this study and also with the Costa formula.Fitting by the proposed formula has an  [37] with some published formulae, including the Einstein equation [30], the Batchelor formula [31], the Roscoe formula [32], the Krieger and Dougherty formula [21], the Chong formula [22], the Dabak formula [23], the Liu formula [24] and the Boyer formula [25].Here φ m is taken as 0.67 in this study.
Table 2. Goodness of fit for some published formulae (compared with published measured data in Figure 5).The second group of published data was mainly conducted by van der Molen and Paterson [36].They investigated the deformation of partially melted granite with solid concentrations from 75% to 95% at 800 and 300 MPa confining pressure and measured viscosity-solid fraction values for partially melted granite.Their measurement results suggested that the critical crystal concentration differentiating granular-framework-controlled flow behaviour from the suspension-like behaviour could be around 0.72.Measured viscosity-crystal fraction values at low solid concentrations for partially crystallized Mg 3 Al 2 Si 3 O 12 melt by Lejeune and Richet [37] (for the Mg 3 Al 2 Si 3 O 12 melt, crystals are well rounded), and other viscosity-solid fraction values for uniform spherical particles by Thomas [42] were also incorporated into this group of published data in order to make a simple comparison between the available experimental results and the proposed equation [35], although these experimental data are limited and their sources are not coherent.
Figure 7 shows the comparisons of these measurement data with the proposed formula in this study and also with the Costa formula.Fitting by the proposed formula has an R 2 value of 0.98 and an NRMSE value of 0.03, whereas the Costa formula only has an R 2 value of 0.82 and an NRMSE value of 0.10.As solid concentration φ increases from 0 to 0.3, the Costa formula agrees with the measured viscosity value slightly better than the proposed formula.For 0.3 < φ < 0.69, the Costa formula overestimates the viscosity values, whereas it underestimates the viscosity value for 0.69 < φ < 0.83.When solid concentration continually increases from 0.83 to 0.95, the Costa formula simply adopts a steady state value of the highly-concentrated suspension viscosity, which do not characterize the continually increasing trend of the viscosity with solid concentration.In this respect, the proposed formula is capable of characterizing this behaviour.
The comparisons of measured data with some published formulae are shown in Figure 8. Table 3 shows goodness values of fitting for these published formulae, among which, the Krieger and Dougherty formula [21] provides the best fit to published experimental data, with an R 2 value of 0.83 and a lowest NRMSE value of 0.12.
Water 2017, 9, 474 11 of 14  < 0.83.When solid concentration continually increases from 0.83 to 0.95, the Costa formula simply adopts a steady state value of the highly-concentrated suspension viscosity, which do not characterize the continually increasing trend of the viscosity with solid concentration.In this respect, the proposed formula is capable of characterizing this behaviour.The comparisons of measured data with some published formulae are shown in Figure 8. Table 3 shows goodness values of fitting for these published formulae, among which, the Krieger and Dougherty formula [21] provides the best fit to published experimental data, with an 2 R value of 0.83 and a lowest NRMSE value of 0.12.7 and other published formulae, including the Einstein equation [30], the Batchelor formula [31], the Roscoe formula [32], the Krieger and Dougherty formula [21], the Chong formula [22], the Dabak formula [23], the Liu formula [24] and the Boyer formula [25].Here m  is taken as 0.95.Water 2017, 9, 474 11 of 14  < 0.83.When solid concentration continually increases from 0.83 to 0.95, the Costa formula simply adopts a steady state value of the highly-concentrated suspension viscosity, which do not characterize the continually increasing trend of the viscosity with solid concentration.In this respect, the proposed formula is capable of characterizing this behaviour.The comparisons of measured data with some published formulae are shown in Figure 8. Table 3 shows goodness values of fitting for these published formulae, among which, the Krieger and Dougherty formula [21] provides the best fit to published experimental data, with an 2 R value of 0.83 and a lowest NRMSE value of 0.12.7 and other published formulae, including the Einstein equation [30], the Batchelor formula [31], the Roscoe formula [32], the Krieger and Dougherty formula [21], the Chong formula [22], the Dabak formula [23], the Liu formula [24] and the Boyer formula [25].Here m  is taken as 0.95.7 and other published formulae, including the Einstein equation [30], the Batchelor formula [31], the Roscoe formula [32], the Krieger and Dougherty formula [21], the Chong formula [22], the Dabak formula [23], the Liu formula [24] and the Boyer formula [25].Here φ m is taken as 0.95.Finally, although the proposed equation in this study has a good fitting relationship with the measured data and some published experimental data, there is an intrinsic limitation: there are four adjustable parameters incorporated within it, and it has no simple form, whereas the previously proposed published formulae do.Future work should be aimed at simplifying the proposed equation and further improving its applicability to fit experimental or observational results.

Concluding Remarks
In this study, a general equation was proposed that describes the relative viscosity of a suspension as a function of suspended solid concentration, covering a range from very dilute to highly concentrated states, based on Costa [35].The viscosity of dilute as well as a highly concentrated kaolinite suspension (in excess of 40% and lower than 50%) was measured in a coaxial cylinder rheometer.A comparison with measured viscosity-solid concentration values showed the validity of the proposed equation.The proposed equation also showed a good fitting relationship with published rheological behaviours of partially crystallised Li 2 Si 2 O 5 and Mg 3 Al 2 Si 3 O 12 .Simplifying the proposed equation and further improving its fit with experimental or observational results should be the goal of future research.

Figure 2 .
Figure 2. Cumulative particle-size distribution for the kaolinite sample.

Figure 2 .
Figure 2. Cumulative particle-size distribution for the kaolinite sample.

Figure 3
Figure 3 shows the measured relative viscosity values of the water-kaolinite suspension for particle concentrations from 1.88 × 10 −3 to 0.48.With increasing  , r  first increased linearly (0 <  < 0.15) (as shown by blue background colour in Figure 3), then showed intermediate increasing behaviour (0.15 <  < 0.42) (as shown by the green background colour), and finally a continuous increasing trend (0.42 <  < 0.48) (as shown by the red background colour).The presence of the continually increasing regime during this experiment may indicate complicated rheological properties of the highly concentrated suspension, which exhibited obvious non-Newtonian characteristics and behaved as a solid with applied stress.

Figure 3 .
Figure 3. Measured relative viscosity values as a function of kaolinite particle concentration.

2 R
root mean square error, NRMSE, defined as

Figure 3 .
Figure 3. Measured relative viscosity values as a function of kaolinite particle concentration.

Figure 5 .
Figure 5.Comparison of measured values (indicated by black circles) for the viscosity-crystal fraction relationship for partially crystallized Li2Si2O5 melt from Lejeune and Richet [37], predicted values from the proposed equation (indicated by red dashed line) and those from the Costa formula (indicated by blue dashed line).Here m  is taken as 0.45 from the experiment results of Lejeune and Richet [37].The experimental temperature is 770.8K. Fitting values of the proposed formula are:  = 0.38,  = 0.00001,  = 20, and  = 1.3.

Figure 5 .
Figure 5.Comparison of measured values (indicated by black circles) for the viscosity-crystal fraction relationship for partially crystallized Li 2 Si 2 O 5 melt from Lejeune and Richet [37], predicted values from the proposed equation (indicated by red dashed line) and those from the Costa formula (indicated by blue dashed line).Here φ m is taken as 0.45 from the experiment results of Lejeune and Richet [37].The experimental temperature is 770.8K. Fitting values of the proposed formula are: α = 0.38, β = 0.00001, γ = 20, and δ = 1.3.

Figure 7 .
Figure 7.Comparison of measured viscosity-solid fraction values for partially melted granite (indicated by blue circle) from van der Molen and Paterson [36], partially crystallized Mg3Al2Si3O12 melt (green circle) from Lejeune and Richet [37], uniform spherical particles (red circle) from Thomas [42] at low solid concentration, and predicted values from the proposed equation (black dashed line) and those from the Costa formula (pink dashed line).Here m  is taken as 0.72 from the experiment results of van der Molen and Paterson [36], Lejeune and Richet [37] and Thomas [42].The experimental temperature in Lejeune and Richet [37] is 1175 K. Fitting values of the proposed equation are:  = 0.9975,  = 0.014,  = 3.5 and  = 6.

Figure 7 .
Figure 7.Comparison of measured viscosity-solid fraction values for partially melted granite (indicated by blue circle) from van der Molen and Paterson [36], partially crystallized Mg 3 Al 2 Si 3 O 12 melt (green circle) from Lejeune and Richet[37], uniform spherical particles (red circle) from Thomas[42] at low solid concentration, and predicted values from the proposed equation (black dashed line) and those from the Costa formula (pink dashed line).Here φ m is taken as 0.72 from the experiment results of van der Molen and Paterson[36], Lejeune and Richet[37] and Thomas[42].The experimental temperature in Lejeune and Richet[37] is 1175 K. Fitting values of the proposed equation are: α = 0.9975, β = 0.014, γ = 3.5 and δ = 6.

Figure 7 .
Figure 7.Comparison of measured viscosity-solid fraction values for partially melted granite (indicated by blue circle) from van der Molen and Paterson [36], partially crystallized Mg3Al2Si3O12 melt (green circle) from Lejeune and Richet [37], uniform spherical particles (red circle) from Thomas [42] at low solid concentration, and predicted values from the proposed equation (black dashed line) and those from the Costa formula (pink dashed line).Here m  is taken as 0.72 from the experiment results of van der Molen and Paterson [36], Lejeune and Richet [37] and Thomas [42].The experimental temperature in Lejeune and Richet [37] is 1175 K. Fitting values of the proposed equation are:  = 0.9975,  = 0.014,  = 3.5 and  = 6.

Table 3 .
Goodness of fit for published formulae based on the comparison with measured viscosity-solid fraction values presented in Figure7.