Macroscopic rheology of non-Brownian suspensions at high shear rates: the influence of solid volume fraction and non-Newtonian behaviour of the liquid phase

Modelling the macroscopic rheology of non-Brownian suspensions is complicated by the non-linear behaviour that originates from the interaction between solid particles and the liquid phase. In this contribution, a model is presented that describes suspension rheology as a function of solid volume fraction and shear rate dependency of both the liquid phase, as well as the suspension as a whole. It is experimentally validated using rotational rheometry (φ\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varphi$$\end{document} ≤ 0.40) and capillary rheometry (0.55 ≤ φ\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varphi$$\end{document}  ≤ 0.60) at shear rates > 50 s−1. A modified Krieger-Dougherty relation was used to describe the influence of solid volume fraction on the consistency coefficient, K\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$K$$\end{document}, and was fitted to suspensions with a shear thinning liquid phase, i.e. having a flow index, n\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$n$$\end{document}, of 0.50. With the calculated fit parameters, it was possible to predict the consistency coefficients of suspensions with a large variation in the shear rate dependency of the liquid phase (n\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$n$$\end{document} = 0.20–1.00). With increasing solid volume fraction, the flow indices of the suspensions were found to decrease for Newtonian and mildly shear thinning liquid phases (n\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$n$$\end{document} ≥0.50), whereas they were found to increase for strongly shear thinning liquid phases (n\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$n$$\end{document} ≤0.27). It is hypothesized that this is related to interparticle friction and the relative contribution of friction forces to the viscosity of the suspension. The proposed model is a step towards the prediction of the flow curves of concentrated suspensions with non-Newtonian liquid phases at high shear rates.


Introduction
Concentrated suspensions of solid particles in Newtonian and non-Newtonian liquid phases are common in nature and industry. These multiphase systems are often subjected to transport or shaping processes where shear rates typically vary between 10 0 and 10 3 s −1 ( Barnes et al. 1989). Examples include extrusion of confectionaries (Altan and Maskan 2012), 3D-printing of food (Wilms et al. 2021a) and non-food pastes (Buswell et al. 2018), pumping of cement slurries that are used in the drilling of oil wells (Tao et al. 2020) or the pumping of blood through arteries (Stroev et al. 2007). The prediction and control of these processes are complicated by the variation in composition of the suspensions. Minor differences in, for instance, the solid volume fraction can lead to a substantial increase in pressure requirements, leading to complications in both humans (Jóźwik-Plebanek et al. 2020) and machines (El Cheikh et al. 2017). The modelling and optimization of these processes not only requires a description of the flow behaviour of the suspensions at process relevant conditions, but also as a function of their composition.
In dilute suspensions, the presence of particles not only distorts the flow field of the liquid but also introduces friction at the particle's surface, both contributing to an increase in the viscosity (Mewis and Wagner 2012). This increase is linear below a volume fraction of ≈0.01 (Einstein 1911;Rutgers 1962), but rises exponentially when the particle concentration increases further, as a result of the hydrodynamic, multibody interaction between neighbouring particles. At a certain solid volume fraction, the viscosity diverges to infinity, a point that is referred to as the maximum packing fraction of the material (Krieger and Dougherty 1959).
Numerous equations have been proposed to describe the viscosity as a function of solid volume fraction, often relative to the viscosity of the pure liquid phase (Mooney 1951;Maron and Pierce 1956;Krieger and Dougherty 1959;Mendoza and Santamaría-Holek 2009;Pal 2020). The majority of these models has been applied to suspensions with a Newtonian liquid phase (Konijn et al. 2014;Pal 2020) and if non-Newtonian behaviour was present, it was only qualitatively described by incorporating it in fit parameters (Mueller et al. 2010). By using qualitative descriptions, the effect of nonlinear rheology of the liquid phase on the rheology of the suspension tends to be forgotten. Suspensions with a shear thinning phase, for example, show a smaller increase in the relative viscosity upon addition of solid particles. Since solid particles are resistant to deformation, an increase in their volume fraction will limit the shear within the system to increasingly small volumes (Barnes 2003;Vázquez-Quesada et al. 2016). The resulting local shear rates are higher, which leads to a more pronounced shear thinning of the liquid phase and hence a smaller increase in the relative viscosity of the suspension (Pal 2015). As the shear rate dependency of the liquid phase varies with formulation (e.g. thickener concentration (Wyatt and Liberatore 2009)), the influence on the viscosity can only be predicted a priori if a generalized function is known. Currently, only a small number of predictive models include non-linear rheology of a liquid phase in their description of the viscosity as a function of solid volume fraction. These models are often theoretical, with limited support of experimental data at solid volume fractions that approach the maximum packing (Jarzebski 1981;Tanner et al. 2010;Pal 2015).
The nature of the shear rate dependency of concentrated multiphase systems complicates the description into a single function. Not only the liquid phase, but also the solid particles contribute to non-linearities of the rheological response (Mueller et al. 2010). This non-linear behaviour is generally linked to changes in the material's microstructure, i.e. the position and orientation of suspended particles relative to each other (Pabst 2004;Stickel and Powell 2005). As the microstructure varies depending on, among others, the flow conditions, suspensions often show a different shear rate dependency than their pure liquid phases (Mueller et al. 2010;Denn et al. 2018;Papadopoulou et al. 2020;Wilms et al. 2021b). Considering the range of shear rates observed in practical flow situations (Barnes et al. 1989), the accuracy of predictive models will increase by including a shear rate dependency of not only the liquid phase, but also of the suspension as a whole.
In this contribution, a model describing the viscosity as a function of both the solid volume fraction as well as the rheological behaviour of the liquid phase is verified and a step towards the incorporation of a shear rate dependency of the resulting suspension is made. After quantifying the yield stress, the description of the macroscopic rheology of the suspension is reduced to an Ostwald-de Waele function. The influence of solid volume fraction and liquid phase rheology on the consistency coefficient and shear rate dependency, i.e. flow index, are then discussed in turn. The former is modelled using a modified Krieger-Dougherty model, based on a proposition by Laven and Stein (1991) and the latter using an empirical equation with a single fit parameter. The models are fitted to experimental data that is obtained using a rotational rheometer, for low solid volume fractions, and a high-pressure capillary rheometer (HPCR), for high solid volume fractions. The HPCR is equipped with a 3D-printed die with a grooved interior surface, following the method proposed by Wilms et al. (2021b). The grooved surface prevents wall slip and allows for the direct determination of the shear viscosity. The accuracy of the proposed model is verified for liquid phases having flow indices between 0.20 and 1.00. Brownian motion can be neglected over the entire range of measurement conditions (Peclet number > 10 3 ) and no shear thickening has been observed.

Theory
The basis of the rheological characterisation is the description of the flow curve, i.e. the shear stress as a function of shear rate ̇ . To include a shear rate dependency, the Ostwald-de Waele relationship (power law) is used; where K is the consistency coefficient and n the flow index. Equation (1) can be used to describe the complete flow behaviour, yet it is often converted to the shear viscosity, = ∕̇ , which found more practical use; To relate the viscosity to the solid volume fraction, , i.e. ( ) , the Krieger-Dougherty equation is used (Krieger and Dougherty 1959); where r is the relative viscosity, and l the viscosities of the suspension and liquid phase, respectively, m the maximum packing fraction and [ ] the intrinsic viscosity. The intrinsic viscosity depends on, among others, the particle shape and is generally assumed to be 2.5 in case of rigid, non-interacting spheres (Wierenga and Philipse 1998).
(1) = K̇n Equation (3) can be extended to include the non-linear behaviour of the liquid phase, by incorporating its flow index n l ; Equation (4) is one of a set of equations proposed by Pal (2015), albeit slightly modified, now including the intrinsic viscosity as a fit parameter, rather than a fixed value of 2.5. The term between brackets has been originally derived by Laven and Stein (1991) and is related to the reference viscosity of the liquid phase that is used to calculate the relative viscosity. In practice, the relative viscosity is defined as the ratio of the suspension viscosity at a given shear rate to the viscosity of the liquid phase at the same shear rate (constantshear-rate relative viscosity) (Pal 2015). In reality, the local shear rates of the liquid phase are significantly higher inside the suspension and the true relative viscosity is different. The term between brackets is required to correct the intrinsic viscosity of the constant-shear-rate relative viscosity for this difference.
As the viscosity depends on the shear rate (Eq. (2)), any difference in the flow index between the suspension and the liquid phase will not be accounted for in Eq. (3) or Eq. (4). Knowing that the flow index can decrease with increasing solid volume fraction (Mueller et al. 2010;Denn et al. 2018;Papadopoulou et al. 2020;Wilms et al. 2021b), Eq. (4) is not a unique function. The influence of solid volume fraction is, therefore, best compared when the viscosity is independent of shear rate, i.e. at a shear rate of 1 s −1 , in which case the viscosity is equal to the consistency coefficient, K , Eq. (2) (Mueller et al. 2010;Konijn et al. 2014). Within this document, the relative consistency coefficient K r is therefore used instead of r in the (modified) Krieger-Dougherty equation, where K r = K∕K l , with K consistency coefficient of the suspension and K l the consistency coefficient of the liquid phase.
In suspensions, a yield stress may arise when particles form a percolating network that is able to resist gravitational forces (Mahaut et al. 2013). The yield stress y can be included into Eq. (1) to form the Herschel-Bulkley equation; Similar to the viscosity, the yield stress can be described as a function of solid volume fraction, as proposed by Heymann et al. (2002).
where * is a fitting parameter that describes the yield stress at = m Krieger-Dougherty equation, or, to be more precise, a modified Maron-Pierce equation, where the exponent −[ ] m is a fixed value of − 2. The yield stress is included in the analysis, not only because it could influence the flow curve, but also because it may affect the velocity profile and consequently the choice of equations to determine the wall shear rate from experimental data (Kelessidis et al. 2006).

Suspension preparation
All aqueous liquid phases were made by dispersing the thickeners in regular tap water and gently stirring for 4 h, except for HPMC. In the case of HPMC, the powder was first dispersed into half of the desired amount of water (warm, ~ 75 °C), before adding the remaining water (cold, ~ 20 °C). For all aqueous liquid phases, clear samples were obtained after 48 h at 8 °C. The wax dispersion was made by melting and recrystallizing the oil-wax mixture inside a ploughshare mixer, as described by Wilms et al. (2020). The concentration and density of the liquid phases are shown in Table 1. Suspensions (150 mL) with solid volume fractions > 0.54 were made using a KitchenAid (Artisan 5KSM150, Whirlpool Company, USA), operated at speed 1. After 5 min of mixing, the homogeneity of the suspensions was increased by subjecting them to 5 s of high shear inside a hand blender (Multiquick 5, Braun GmbH, Germany) equipped with a chopper. All but the HPMC suspensions with a solid volume fraction < 0.55 were first mixed manually, before performing an identical high-shear step inside the hand blender. In the case of HPMC, to prevent incorporation of air, the suspensions were made by mixing the solids and thickener before adding hot water. The suspensions were stirred until the material did not phase separate at rest. For PEG, sodium alginate and HPMC suspensions, the solid volume fraction was varied between 0.10 and 0.60, whereas for xanthan, gum arabic and wax suspensions only 0.60 was made. To prevent the evaporation of moisture, the suspensions are immediately stored in an airtight container and contact with surrounding air is minimized during measurements. Two batches were produced of each suspension, to account for production variation. The air content of the suspensions with solid volume fractions ≥ 0.58 was determined using the Archimedes principle (Wilms et al. 2020) (Table 1). For suspensions with solid volume fractions < 0.58, the samples were too liquid to perform the measurements.

Rotational rheometry
The viscosity of the liquid phase was measured using a rotational rheometer (Kinexus ultra, Malvern Instruments GmbH, Germany) equipped with a Couette geometry consisting of a serrated cup (diameter, 27.5 mm) and serrated bob (diameter, 25.0 mm) with a conical tip (15 0 ) at a gap height of 9.15 mm. The samples were tempered at 22 °C for 15 min, before being pre-sheared at 80% of the maximum shear rate that was measured for the particular sample, to avoid an extensive period of rotation close to/at the maximum torque for the concentrated suspensions. For most samples, this corresponded to a pre-shear at 400 s −1 for 30 s. The flow curve was determined using a shear rate ramp with a stepwise reduction of shear rates from 500 to 0.01 s −1 , following a logarithmic scale. Only for the HPMC suspension with a solid volume fraction of 0.40, a pre-shearing at 80 s −1 and maximum shear rate of 100 s −1 were used instead, due to torque limitations. Each shear rate was maintained until a steady state was reached, with a minimum time of 1 / ( Mezger 2006). No thixotropic effects were observed. To correct for non-Newtonian behaviour of the material, the true shear rate at the wall of the bob, ̇w , was calculated using Eq. (7) (Steffe 1996); where Ω is the angular velocity of the bob in radians per second, n the slope of plots describing the logarithm of the shear stress versus the logarithm of the angular velocity (straight line) and is the ratio between the radius of the cup and the radius of the bob. The measurement points with a true shear rate > 50 s −1 were fitted to Eq. (1), as sodium alginate suspensions were still in their first Newtonian plateau around 10 s −1 , at low solid volume fractions. The presence of flow instabilities, and in particular Taylor vortices, can be estimated using the Taylor number, Ta , calculated as , with the density and R i and R o the radius of the inner and outer cylinder respectively (Ewoldt et al. 2015). For the liquid phase with the lowest viscosity at the highest shear rate, i.e. gum arabic, the Taylor number never exceeded 600. Assuming a critical Taylor number of 1708 (Dou et al. 2008), flow instabilities are considered to be absent. Note that for HPMC samples, rod climbing was observed at higher shear rates (> 800 s −1 ), due to the Weissenberg effect. Following Stokes law, the settling velocity of the limestone particles was in the order of 10 −8 m/s and sedimentation was considered to be negligible during measurements. A possible way to strengthen this assumption involves the determination of the so-called Shields number, which is originally used to determine the motion of sediment in currents. The Shields number has already been applied to flows between parallel-plates (Tanner and Dai 2016; Martone et al. 2020), but not for concentric cylinders. Lacking boundary values, this is considered to be out of scope for this contribution.

Oscillation rheometry
The yield stress was measured using oscillation rheometry (Kinexus ultra, Malvern Instruments GmbH, Germany), using parallel plates (diameter of 40.0 mm) with a serrated surface. A gap height of 1 mm was used and the sample exterior was coated with silicon oil (10,000 cSt) to prevent water evaporation. After a 5-min waiting step, stress sweeps were performed between 0.1 and 1000 Pa at a frequency of 10 rad/s (~ 1.6 Hz). The yield stress was defined as the point at which the storage modulus has decreased by 30% from that observed in the linear viscoelastic regime. This point is chosen to prevent the ambiguity of choosing the range that is used for the power law fit to the decreasing storage modulus, which does not apply to the linear viscoelastic regime. Measurements were performed in duplicate on suspensions with a solid volume fraction between 0.40 and 0.60.

Capillary rheometry
The shear viscosity was measured using a twin-bore highpressure capillary rheometer (HPCR) (RH 2000; Malvern Instruments GmbH, Germany) according to Wilms et al. (2021b). The barrels were equipped with a 10 MPa (left) and 1.5 MPa (right) pressure transducer. To perform the measurements a 'rough' die was used with a radius of 0.50 mm and length to radius ratio of 32. The die was made using Direct Metal Laser Sintering (DMLS) of stainless steel (GP1) (Speedpart GmbH, Germany). The internal surface of the die consisted of an equilateral saw tooth pattern along the die length. The teeth had a depth of 0.4 mm and were 0.4 mm apart. It is assumed that the roughness completely prevents wall slip. The parallel alignment of roughness elements is thought to prevent any additional flow in between the teeth, i.e. the flow cannot extend over more than the distance between successive teeth (Wilms et al. 2021b). This is in contrast to other roughened surfaces that present a porous structure, for which an additional slip velocity has been proposed in the case of Newtonian (Carotenuto et al. 2015) and non-Newtonian materials (Paduano et al. 2019). For a more elaborate discussion on the use of a roughened extrusion die and the no-slip condition please refer to Wilms et al. (2021b). A drilling tip with a thickness of 1.00 mm and a metal thread with a thickness of 1.005 mm (Malvern Instruments GmbH, Germany) were used to inspect the radius of the internal flow channel prior to every set of measurements.
For all experiments, the drilling tip could penetrate the die, whereas the metal thread did not, indicating little to no surface wear. The entry pressure loss was determined using an orifice die, with a length to radius ratio of ~ 0 and all Bagley plots are thus assumed to be linear (Bagley 1957). Material is filled into 20 mL disposable syringes by punching, thereby minimizing the inclusion of air, before transferring them into the temperature-controlled barrels (22 °C). After pre-compression at a rate of 10 mm s −1 until material discharge is observed at the die, the material is allowed to rest for 5 min. For the measurements, a table of apparent shear rates was used (20, 40, 80, 160, 320, 640, 960, 1280 s −1 ), which was controlled by the volumetric flow rate ( Q ) through the die, following Eq. (8).
The term 'apparent shear rate' indicates that the corresponding value is the expected shear rate, when assuming Poiseuille flow. For the PEG and gum arabic suspensions, the shear rate was increased until the maximum extrusion pressure (10 MPa) was reached, e.g. for PEG (0.60) the maximum shear rate was 320 s −1 . As all suspensions were prepared in duplicate, they are separately measured once. Prior to every measurement, the rough die was cleaned using an interdental brush and a small layer of oil was applied to the orifice die's surface. The oil layer reduced the sticking of material to the orifice die at low shear rates. The final 10 values (corresponding to 5 s) at each shear rate stage were averaged and used for further calculations. To obtain the true pressure drop over the rough die, the results from the orifice die were subtracted from the rough data. The wall shear stress w was then calculated from this pressure drop ( ΔP ) using Eq. (9); With R is the radius and L is the length of the die. The wall shear rate is calculated from the apparent shear rate (Eq. (8)) using the Weissenberg-Rabinowitsch correction (Eq. (10)) (Rabinowitsch 1929;Steffe 1996). This correction is required to account for the non-parabolic velocity profile as a result of non-Newtonian behaviour of the material, similar to Eq. (7).
The shear viscosity then follows from = w ∕̇w and flow curve was fitted to Eq. (1).
Fits to exponential functions are performed by using the logarithm of the dependent variable, to avoid overfitting to datapoints at higher solid volume fractions. All data was analysed in Origin (OriginPro 2020, OriginLab, USA). The relative standard error of the fitted parameters to the data obtained using the capillary rheometer ( ≤0.55) was on average 4.9% for K and 2.2% for n , with wax and gum arabic showing significantly higher values of 18% and 16% for K and 7.6% and 4.0% for n respectively. For the data obtained using the rotational rheometer ( ≤0.40), the average relative standard error was lower, i.e. 2.8% for K and 1.5% for n . In addition, the relative standard deviation of the pressure reading between two batches at a specific shear rate never exceeded 13% and was on average 3.1%.

Liquid phase rheology -K l and n l
The used thickeners vary in their rheological properties (Table 2). In the case of PEG and gum arabic, the viscosity is independent of shear rate, indicating Newtonian behaviour with a flow index n l of 1.0. For all other thickeners, varying degrees of shear thinning are observed, being the most pronounced in the case of HPMC ( n l = 0.20). Note that the coefficient of determination (R 2 ) was > 0.99 for all thickeners, apart from xanthan (R 2 = 0.98) and that the values are different from a previous study with similar liquid phases (Wilms et al. 2021b), due to the correction for the non-Newtonian behaviour of the data obtained with the rotational rheometer.
In this contribution, it is initially assumed that the parameters described in Table 2 are valid within the entire range of experimental conditions. This includes the high local shear rates of the liquid phase in between solid particles. For suspensions, the local shear rate can be roughly estimated by dividing the macroscopic, apparent, shear rate by the difference between total volume and the packing of the solid particles ( ̇l ocal =̇a pparent ∕(1 − ∕ m ) (Windhab 1986). At a solid volume fraction of 0.60 and assuming a m of 0.64, the microscopic, local, shear rate of the liquid phase would be approximately 15 times the apparent shear rate of the suspension. Dakhil et al. (2019) have found an infinite shear viscosity at shear rates ~ 6 × 10 4 s −1 for 1 wt.% xanthan dispersions, increasing with increasing thickener concentration. Although no information is available for the other thickeners, extrapolating the data of Table 2 to apparent shear rates of the suspension in the order of ~ 10 3 and liquid phases in the order of ~ 10 4 appears reasonable. This assumption is reexamined later in this document.

Yield stressy
The importance of a yield stress is evident at low shear rates, when it describes a large part of the total shear stress (Eq. (5)). Although the importance decreases with increasing shear rates, the yield stress can still influence the flow profile, especially in pressure driven flows (Kelessidis et al. 2006). As the shear stress linearly decreases from the channel wall towards the interior, a shear discontinuity can develop in the form of an unsheared, plug flow, region (Kelessidis et al. 2006;Ovarlez et al. 2009). Quantification of the yield stress is required to choose the correct equation to calculate the shear rates at the wall of the specific flow gap (Steffe 1996;Kelessidis et al. 2006;Kelessidis and Maglione 2008).
The yield stress as a function of solid volume fraction for PEG and HPMC suspensions is well-described by Eq. (6) (Fig. 1). The fitting parameters * and m become 0.361 and 0.660, in the case of HPMC, and 0.310 and 0.659, in the case of PEG, both having an R 2 of 0.97. The values are similar to values reported by Heymann et al. (2002) and Mueller et al. (2010), which was already discussed in a previous study for the same suspensions (Wilms et al. 2021b). The yield stress increases exponentially with solid volume fraction and absolute values are in the same order of magnitude, with the apparent exception of sodium alginate and gum arabic. For sodium alginate and gum arabic, the relaxation time between gap adjustment and measurement start was found to have a significant influence on the measured yield stress. Only after increasing the waiting step from 5 to 30 min, the yield stress of sodium alginate resembled the values Table 2 Parameters of the Ostwald-de Waele function (Eq. (1)), fitted to rotational rheometry data between 50 and 500 s −1 for liquid phases with different thickeners. The standard error of the parameters is included and the coefficient of determination (R 2 ) was > 0.99 for all thickeners, apart from xanthan (R 2 = 0.98) found for the other suspensions (star, Fig. 1), whereas gum arabic only increased to 9.45 ± 0.54 Pa after 150 min (not shown). The choice of relaxation time depends on the practical conditions to which the yield stress values are thought to apply. It is likely that the materials have experienced stresses greater than the yield stress in the case of capillary extrusion and many other applications that involve a converging flow into the flow channel (Ardakani et al. 2011). Longer relaxation times are not further explored as they have little practical relevance for the flow during target processes. The yield stress of the suspensions did not exceed 80 Pa at the highest solid volume fraction (Fig. 1) and was well below the shear stresses that are observed during HPCR measurements (e.g. Fig. 2). The description of the flow curve of the suspensions can thus be simplified to an Ostwald-de Waele function Eq. (1), at the used measurement conditions. As the relative contribution of unsheared regions to the velocity profile will be limited, Eq. (10) applies to the suspensions in this study. Although it is possible to determine the rheological parameters from the pressure drop using constitutive equations that include the yield stress (Kelessidis et al. 2006;Ardakani et al. 2011;Ahuja et al. 2018), this involves a numerical approximation for materials that can be described by Eq. (5). The analytical solution of Eq. (10) is therefore preferred.

Function of solid volume fraction
The measurement data of the rotational and capillary rheometers is used to determine the flow curves of the suspensions, after correcting the apparent shear rates for the non-Newtonian behaviour (Eq. (7) and Eq. (10)). The flow curves are fitted to Eq. (1) at shear rates > 50 s −1 , focussing on the power law region whilst ensuring the use of data points over at least one decade (Fig. 3). An overview of the fitted consistency coefficients and flow indices is shown in Table 3. The coefficient of determination for all fits was > 0.98 with the exception of wax ( m = 0.60) having an R 2 of 0.94 and HPMC with m = 0.10, 0.20 and 0.40 having an R 2 of 0.94, 0.96 and 0.84 respectively.
The relative consistency coefficients are plotted as a function of solid volume fraction and fitted to the original Krieger-Dougherty equation (Eq. (3)) (Fig. 2). The individual datasets are well-described (R 2 > 0.99), as expected when using two fit parameters. The height of the curve shifts downward with decreasing flow index of the liquid phase, as the high local shear rates lead to a more pronounced shear thinning and a smaller increase in consistency coefficient at larger solid volume fractions (Pal 2015). In the case of PEG and sodium alginate, the maximum packing fractions are similar, i.e. 0.642 and 0.638 respectively, suggesting that  the effect of a shear thinning liquid phase on the consistency coefficient can be described by the intrinsic viscosity, which is the basis for Eq. (4). The intrinsic viscosity for the PEG suspensions (4.64, Fig. 2) is significantly larger than 2.5 that is generally assumed for suspensions with Newtonian liquid phases. The particle shape has often been linked to changes in the intrinsic viscosity, but the observed difference would require a large shape anisotropy (Wierenga and Philipse 1998). The aspect ratio of the limestone particles was only 1.55 ± 0.38, as determined with the image analysis of micrographs that are made using a scanning electron microscope (JSM-IT100, Jeol Ltd. Japan). Approximately 4000 particles are analysed using a Matlab script (Matlab R2020b, MathWorks, USA) with an advanced watershed function that is described by van Dalen and Koster (2011). A few other studies have described similarly high intrinsic viscosities for near-spherical particles, both experimental (Pabst et al. 2006), as well as theoretical (Happel 1957). Pabst et al. (2006) obtained intrinsic viscosities of 4.72 and 4.77 for wheat-and corn starches respectively and argued that attractive interactions can lead to aggregation at higher solid volume fractions. For limestone particles (~ 5 µm), the formation of aggregates cannot be excluded, possibly enhanced by the depletion attraction resulting from the high concentration of macromolecules (Bessaies-Bey et al. 2018). As the intrinsic viscosity in the dilute limit will naturally describe the interaction of a single particle in a fluid, aggregation implies a changing value with increasing particle concentration and even shear rate (Brenner 1974). The intrinsic viscosity, as used in Eq. (4), can thus be regarded as a best fit average.
Whilst the intrinsic viscosity is a characteristic of the particle system on a small scale (Pabst et al. 2006), the maximum packing fraction is a measure for the suspensions microstructure as a whole (Stickel and Powell 2005). The microstructure naturally changes with flow conditions, and the maximum packing is thus no fixed value (Wildemuth and Williams 1984). Considering the maximum packing fraction as the threshold above which the particles form a percolating network, i.e. analogue to the yield stress, its value varies between the random loose packing and the random close packing (Song et al. 2008). Assuming similar flow conditions and shear stresses that are significantly above the yield  stress, the maximum packing fraction of suspensions should be independent of the liquid phase. This is confirmed for PEG and sodium alginate suspensions, whereas HPMC suspensions show a significantly lower m , i.e. 0.614 (Fig. 2). In the case of HPMC, a possible source of error is the large volume fraction of air (Table 1) and its changing volume under the influence of the extrusion pressure, following the ideal gas law. Although the incorporation of air is often associated with an increase in the viscosity, e.g. for foams, it was found to decrease the viscosity for concentrated suspensions (Kalyon et al. 1991;Kalyon and Aktas 2014). The exact origin of this decrease is unknown, but one could argue that it is related to the decreasing effective volume fraction of the solids, due to the presence of air. As the interparticle distance increases with decreasing solid volume fraction, the viscosity of the suspension decreases for solid fractions where interparticle friction is important. The effect of air is not limited to the consistency coefficient, but affects the complete rheology, including the shear rate dependency. This is briefly discussed later in this document. By mixing the pastes under vacuum, the influence of air can be confirmed in future studies. An additional source of error is presented by the small cavities that are inherent to the design of the capillary rheometer. Small cavities exist between the barrel and (i) the capillary die, (ii) the piston and (iii) the pressure transducer. These cavities can result in drainage of the liquid phase, which can be directly observed by inspecting the die and piston in between measurements. Any drainage will lead to drying of the matrix and growth of the static zone in front of the die entry, thereby increasing the extrusion pressure. Following Darcy's law, this drainage will depend on the absolute pressure and the viscosity of the draining liquid. The PEG suspensions are most susceptible to this drainage, given the low viscosity of the liquid phase and quick increase of extrusion pressure with increasing shear rate, as it has the highest flow index. Although it could not be prevented, nor quantified, the possibility of an error in the consistency coefficients of the PEG suspensions in Fig. 2 should be considered. Drainage of the liquid phase will lead to a larger solid volume fraction, i.e. the true measurement points could be shifted to the right, and to an increase in entry pressure that is not accounted for by the orifice die, i.e. the true measurement point could be shifted downward.
An important implication of determining the shear rate at the wall using Eq. (10) is that it assumes a homogeneous bulk. It is known that pressure-driven flows are subject to shear-induced migration of particles from the wall towards the interior of the flow channel, because of the inhomogeneous shear field (Leighton and Acrivos 1987;Hampton et al. 1997;Ramachandran 2007). This particle migration results in a concentration gradient that could possibly lead to an apparent slip effect and lower apparent viscosity (Wilms et al. 2021b). The rheological parameters in Table 3 thus describe the macroscopic rheology of the material under the considered flow conditions. As the migration is inherent to flow and the dimensions of the die are commonly used for HPCR measurements, further elaboration on this topic is considered to be out of scope for this contribution. For additional information on shear-induced migration, please refer to Semwogerere et al. (2007) and Guazzelli and Pouliquen (2018).

Function of solid volume fraction and liquid phase rheology
To include the non-linear behaviour of the liquid phase into the function of the relative consistency, the validity of Eq. (4) is tested. The initial fit of the maximum packing fraction and the intrinsic viscosity is performed for sodium alginate, being the least affected by air and drainage (Fig. 4). The relative consistency of the PEG and HPMC suspensions is then predicted from the flow index of their liquid phase ( n l, Table 2). Considering the effects of drainage and air, the predicted curves satisfactorily describe the experimental data, for both PEG and HPMC. This is true for the entire range of solid volume fraction, including solid volume fractions → m . Although it has been theoretically predicted, this is the first time that the accuracy at such high solid volume fractions has been verified experimentally (Laven and Stein 1991;Pal 2015). To further improve the fit, the experimental difficulties of drainage and air inclusion into the suspension need to be addressed, which will be part of future work.  Table 2 The universal applicability of Eq. (4) is tested for three additional suspensions with a solid volume fraction of 0.60. The used thickeners, i.e. xanthan, gum arabic and wax (in oil), add additional variation to polymer concentration, flow index and mode of thickening. Compared to the polymers, the thickening of the wax is based on the crowding of solid fat particles, possibly affecting the packing of other solid particles inside the suspension. The relative consistency coefficients of all suspensions are well-described (Fig. 5). The results suggest that Eq. (4) can be used to calculate the consistency coefficient of the suspension when the liquid phase changes, provided the rheology of the liquid phase is known and the solid phase properties remain the same.

Flow index -n
The consistency coefficient is only part of the description of the flow curve (Eq. (1)) and does not describe the shear rate dependency of the suspension. This shear rate dependency is described by the flow index, which is found to decrease with increasing solid volume fraction for PEG, gum arabic, wax and sodium alginate, i.e. the suspensions become more shear thinning (Fig. 6). Part of decrease in the flow index can be explained by the particle shape anisotropy, which allows for alignment of particles in the direction of flow (Barnes 2000), and the breakdown of agglomerates under the influence of shear (Papadopoulou et al. 2020). It is, however, unlikely that both will lead to such pronounced shear thinning behaviour, especially in the case of the Newtonian PEG (n~ 0.75 at φ = 0.60), as the aspect ratio was only 1.55 and the particles are non-Brownian, i.e. Peclet numbers are large (> 10 3 ). Similar shear thinning behaviour has been observed for non-Brownian spheres in Newtonian media and its origin is subject to debate (Mueller et al. 2010;Vázquez-Quesada et al. 2017). Different hypotheses with a hydrodynamic origin have been put forward, suggesting that the high local shear rates could; (i) induce a slip length at the surface (Kroupa et al. 2017), (ii) result in "hidden" non-Newtonian behaviour of the liquid phase (Vázquez-Quesada et al. 2016), (iii) lead to shear-induced organization of particles (Lin et al. 2015b) or (iv) lead to a reduction of the liquid phase viscosity due to viscous heating (Mueller et al. 2010). More recently, a nonhydrodynamic hypothesis has been proposed that is related to the non-linear behaviour of friction forces between particles with respect to pressure and, more specifically, to the elastic deformation of surface asperities as a function of normal load (Chatté et al. 2018;Lobry et al. 2019;Arshad et al. 2021). Papadopoulou et al. (2020) have tested this hypothesis by varying the polarity of the liquid phase, thereby changing the interaction between the solid particles and the liquid. Suspensions that have a strong interaction between solid particles and the liquid phase did not show shear thinning (e.g. smooth-surfaced glass beads in glycerol), whereas non-interacting systems (e.g. smooth-surfaced glass beads in mineral oil) did. It is hypothesized that, in interacting systems, the particles are covered by a thick solvation layer that prevents direct solid-solid contact and thereby prevents shear thinning behaviour, provided the layer covers the surface asperities (Papadopoulou et al. 2020). The absence of contacting surface asperities could also explain the similarity between the flow index of the suspension and the liquid phase, as observed by, for instance, Ovarlez et al. (2015). Their liquid phase consisted of an emulsion, of which the droplets were large compared to the particle roughness, thereby preventing frictional interactions. From Fig. 6, two observations can be made; (i) with increasing solid volume fraction, the flow index decreases for suspensions with liquid phases that have a large flow index (PEG, gum arabic, wax and sodium alginate), but shows little change and possibly a slight increase for suspensions with liquid phases that have a small flow index (xanthan and HPMC) and (ii) the absolute change of the flow index with increasing solid volume fraction is larger for suspensions with a large flow index, i.e. PEG > sodium alginate > HPMC. Let us assume that interparticle friction indeed strongly contributes to the viscosity of the suspensions over a wide range of volume fractions (Tanner 2018) and that the friction coefficient decreases with increasing normal load (Chatté et al. 2018;Lobry et al. 2019). The viscosity of the suspension then becomes a combination of two contributions; the contribution of the liquid phase (hydrodynamic) and the contribution of interparticle contact and resulting friction (frictional). Similarly, the total stress experienced by the suspension becomes a combination of hydrodynamic stress and contact stress (Lin et al. 2015a). As the hydrodynamic-and contact stress both have a unique relation that describes their shear rate dependency, the shear rate dependency of the total stress varies depending on the relative importance of both contributions. More specifically, depending on the relative importance of both contributions, the flow index of the suspension takes on a value that is intermediate between the flow index of the hydrodynamic dissipation and the relation describing the friction coefficient of solid-solid contact.
Starting from zero-friction, in the case of pure liquid, the nature of frictional forces implies an increasing relevance with increasing solid volume fraction and decreasing interparticle distance (Mari et al. 2014). With increasing importance of interparticle friction, the flow index of the suspension can then either increase or decrease, depending on the relative magnitude of the flow indices of the hydrodynamic ( n h ) and frictional ( n f ) contributions. The flow index of the suspension would thus decrease with increasing solid volume fraction when n h > n f and vise versa. This is clarified by plotting the relative flow index, n r , i.e. flow index of the suspensions, n , divided by the flow index of the liquid phase, n l , as a function of solid volume fraction (Fig. 7).
A changing interparticle friction would be able to account for both decreasing as well as increasing flow indices of concentrated non-Brownian suspensions. This variation is more difficult to explain from a pure hydrodynamic perspective. Although wall slip (Kroupa et al. 2017) and viscous heating (Mueller et al. 2010) can explain a decrease of the flow index, they cannot predict an increase of the flow index in the case of suspensions with strongly shear thinning liquid phases. It is unlikely that either slip or viscous heating become inversely proportional to the solid volume fraction, as they both depend on the addition of particles to become effective and are directly proportional by nature. Hidden non-Newtonian behaviour (Vázquez-Quesada et al. 2016) and shear-induced organization (Lin et al. 2015b) would be able to explain individual cases, but there is yet to be a physical explanation how they could account for the separation between liquid phases with varying degrees of shear thinning, as shown in Fig. 7. As the rheology of the liquid phase was only measured until macroscopic shear rates of ~ 500 s −1 , any extrapolation to larger values introduces uncertainty. However, considering the fact that the flow curves of the capillary rheometry data are well-described by Eq. (1) over the complete range of shear rates (Fig. 3), the extrapolation seems justified. The flow curves start at macroscopic shear rates ~ 10 s −1 , at which the microscopic shear rates of the liquid phase are still within the range of shear rates measured with rotational rheometry. If hidden non-Newtonian behaviour was present, an apparent change in the slope of the flow curves of the suspensions would be observed. Note that the results in Fig. 6 and Fig. 7 only suggest that, at high solid volume fractions, the relative importance of interparticle friction increases compared to hydrodynamic effects. This does not mean that factors with a hydrodynamic origin have no influence on the flow index of suspensions investigated in this study. In particular at low solid volume fractions ( < 0.40) for HPMC suspensions, it is possible that the initial decrease in flow index has a hydrodynamic origin. The assumption that the shear rate dependency of the hydrodynamic stress is equal to the flow index of the liquid over the entire range of solid volume fractions is likely a simplification and the existence of a more complex function is an interesting topic for future studies. These future studies should also assess the influence of air. Air might not only affect the consistency coefficient, as discussed in relation to Figs. 3, 4 and 5, but also the shear rate dependency. Ovarlez et al. (2015) studied non-colloidal spheres in a concentrated emulsion and suggested that their particles behaved as low friction particles. They hypothesized that the emulsion droplets hindered frictional contacts, which could also apply to air bubbles that are present in the suspension. As HPMC suspensions showed the highest air content (Table 1), they will be most affected, possibly explaining the limited change of the flow index.

Flow curve of the suspension
It is possible to predict the flow curves of the suspensions if one has a description of both the consistency coefficient and flow index as a function of solid volume fraction and liquid phase rheology. In the case of the consistency coefficient, this function is described by Eq. (4) (Fig. 5), but in the case of the flow index, knowledge is still lacking and only purely empirical models exist (Mueller et al. 2010). Mueller et al. (2010) proposed an exponential function that described the decrease of the flow index as a function of solid volume fraction, aspect ratio and maximum packing fraction for suspensions that are made with silicon oil. A similar fit is performed to the data in Fig. 7, using a simplified version; where b is a fit parameter. Equation (11) changes to 1 + b for suspensions where n h < n f , i.e. for strongly shear thinning liquid phases. Compared to Mueller et al. (2010), Eq. (11) describes the relative flow index but does not include an aspect ratio or maximum packing fraction. The absence of the maximum packing fraction agrees with the hypothesis of interparticle friction, as the flow index does not tend to 0 when → m . Equation (11) is purely empirical, with b = 1.83 (R 2 = 0.96) when fitted to sodium alginate and b = 2.45 (R 2 = 0.99) when fitted to PEG (Fig. 7). Equation (4) and Eq. (11) can be included into the Ostwald-de Waele function (Eq. (1)), to obtain; where the sign in the final term on the right side of the equation changes when n h < n f . Equation (12) describes the shear stress as a function of shear rate, solid volume fraction and liquid phase rheology. It includes two parameters (11) that describe the rheology of the liquid phase K l and n l and three fit parameters; m , [ ]and b . The parameters K l and n l are obtained from Table 2, m and [ ] from the fit to Eq. (4) (Fig. 5) and b from the fit to Eq. (11). The fits were only performed on the sodium alginate suspensions. The predicted and experimental values are compared for PEG, sodium alginate and HPMC suspensions at solid volume fractions of 0.40, i.e. measured using rotational rheometry, and 0.60, i.e. measured using capillary rheometry (Fig. 8).
The validity of Eq. (12) is confirmed by the minimal difference between predicted and experimental values for sodium alginate suspensions, being 8% on average and never exceeding 15%. The universal applicability then follows form the fits to PEG and HPMC suspensions. The smallest deviation between predicted and experimental values is found for the PEG suspension with a solid volume fraction of 0.40, having a deviation < 25% over the entire range of shear rates. However, for the other suspensions, the deviations are significantly larger. In the case of PEG with a solid volume fraction of 0.60, the values are underpredicted by a factor 2, i.e. 50% lower, and, in the case of HPMC with a solid volume fraction of 0.40 and 0.60, the values are overpredicting by a factor 1.75 and 2 respectively, i.e. 75% and 100% higher. Although these differences are too large to allow for application of the current model in practical situations, the fit is promising. Especially considering the large variation in rheological properties between the PEG and HPMC suspensions and the sodium alginate reference and the lack of knowledge concerning the changing flow index with increasing solid volume fraction.  The prediction of the flow curve can be improved by considering three sources of error; (i) the presence of air influences the complete rheological behaviour of the suspension, including the consistency coefficient, as it affects the particle packing, and the flow index, as it could reduce the friction between particles. (ii) In the case of strongly shear thinning liquid phases, the flow index is likely to follow a more complex function than Eq. (11), i.e. an initial decrease that is followed by an increase (Fig. 7). A minor change in the flow index can lead to a large change in the shear rate term (Eq. (12)), especially at high shear rates. (iii) In the case of a Newtonian liquid phase, the suspensions suffer from stability issues that lead to concentration fluctuations and larger entry pressure losses. As the consistency coefficient is an exponential function of solid volume fraction, a small change will lead to a large change in consistency coefficient when → m . By assessing these issues in future studies, the accuracy of the prediction will be improved.
In addition, several alterations can be made to account for specific rheological behaviour and/or decrease the experimental effort. If the studied material shows a high yield stress, Eq. (12) can be extended to include such yield stress, following Eq. (5). The yield stress itself will be a function of solid volume fraction, e.g. via Eq. (6), but most likely also a function of shear rate, which can be incorporated using structural kinetics models (Mujumdar et al. 2002). The accuracy of Eq. (12) relies on three fit parameters and the estimation of these parameters requires significant experimental work. A priori predictions of the fit parameters would increase the practical usefulness of the equation. The term describing the relative consistency coefficient (Eq. (4)) appears to be most suitable for theoretical predictions, considering the purely empirical nature of the term describing the flow index (Eq. (11)). Significant progress has been made regarding the prediction of the maximum packing fraction, as a function of polydispersity (Dörr et al. 2013) and particle shape (Gan et al. 2004). These models can serve as rough guidelines until they are translated to real suspensions that show a continuous distribution of particle properties. In the case of the intrinsic viscosity, individual liquid-particle interactions are fairly well understood, leading to good predictions for particles with varying shapes in the absence of multiparticle structures (Wierenga and Philipse 1998; Mueller et al. 2010). If multiparticle structures cannot be excluded, an intrinsic viscosity of 4.5 appears to be a reasonable estimate, following results reported by Pabst et al. (2006) and this study.

Conclusion
Modelling the macroscopic rheology of suspensions as a function of solid volume fraction is required to make a priori predictions regarding the influence of compositional variation on material processing. This modelling is generally performed using semi-empirical equations, such as the Krieger-Dougherty equation, relating the relative viscosity to the solid volume fraction. However, these equations can neither account for differences in the non-Newtonian behaviour of the liquid phase (Pal 2015), nor that of the suspension (Mueller et al. 2010). Theoretical approaches have been proposed to include non-linear behaviour (Tanner et al. 2010;Pal 2015), but the general applicability of these equations has not yet been verified.
In this contribution, a model is proposed for suspensions of which the macroscopic rheology is described by the Ostwald-de Waele relationship. It is experimentally validated and includes the non-linearity of the consistency coefficient, K , and flow index, n . Experimental data was obtained using rotational rheometry ( ≤ 0.40) and capillary rheometry (0.55 ≤ ≤ 0.60) at shear rates > 50 s −1 . The capillary rheometry is performed using custom-made dies with a rough internal surface, following the method of Wilms et al. (2021b). The consistency coefficient is described using a modified Krieger-Dougherty model that includes the non-linearity of the liquid phase, as proposed by Laven and Stein (1991) and reintroduced by Pal (2015). The model is fitted to suspensions in which the liquid phase is thickened using sodium alginate ( n l = 0.50) and the fit parameters are used to predict the relative consistency coefficient of suspensions with a large variation in their liquid phase rheology ( n l = 0.20-1.00). Predicted values were in good agreement with experimental data. The results verify that the relation that describes the relative consistency coefficient of suspensions in non-Newtonian liquid phases originates from a different local shear rate of the liquid phase. This difference can be corrected for by the proposed modification, even when → m .
The flow index of the suspensions was found to decrease with increasing solid volume fraction for Newtonian and mildly shear thinning liquid phases ( n l ≥0.50), whereas they were found to increase for strongly shear thinning phases ( n l ≤0.27). It is hypothesized that this is due to the contribution of interparticle contact and resulting friction forces to the viscosity of the suspension. The total stress that is experienced by the suspension is a combination of hydrodynamic stress and contact stress. Because the hydrodynamic-and contact stress have a different shear rate dependency, i.e. different flow index, the flow index of the suspension becomes intermediate between these values. As the frictional contribution increases with solid volume fraction and the friction coefficient is shown to decrease with increasing normal load (Chatté et al. 2018), the flow index of Newtonian and mildly shear thinning phases decreases with increasing solid volume fraction. In the case of strongly shear thinning liquid phases, the flow index of the hydrodynamic stress is lower than that of the contact stress, leading to an increase of the flow index with increasing solid volume fraction.
The relative change of the flow index was fitted to a purely empirical function and, together with the modified Krieger-Dougherty equation, incorporated into the Ostwald-de Waele function (Eq. (12)). The proposed equation describes the flow curve as a function of solid volume fraction and non-linear behaviour of both the liquid phase, as well as the resulting suspension. Although the difference between predicted and experimental values was large, i.e. showing under-and overpredictions by a factor 2, Eq. (12) can be considered to be a promising first step, especially considering the large variation in rheological behaviour of the liquid phases that are investigated. The fit can be improved in future studies by considering the influence of air, implementing a more complex function of the flow index of suspensions with a highly shear thinning liquid phase and reducing phase separation for suspensions with a Newtonian liquid phase.
The proposed model provides a basis for the prediction of the macroscopic rheology of suspensions with non-Newtonian liquid phases at high shear rates. It contributes towards the goal of anticipating on compositional variation and a priori prediction of changes in process parameters for numerous flow situations in industry and biological systems.