Introduction

Cellulose micro/nanofibrils (CMNFs), or micro/nanofibrillated cellulose, has been a topic of academic interest since the 1980’s due to their unique properties, such as mechanical robustness, barrier properties, high specific surface area, lightness, and complex rheology. More recently, CMNFs has also been perceived as a versatile, sustainable, and biodegradable material that enables developing eco-friendly all-cellulose products. The industrial interest in CMNFs has recently increased also due to the rising number of commercially available CMNF grades and due to innovations that have lowered the production costs of these materials (Klemm et al. 2011; Moon et al. 2011; Lavoine et al. 2012; Moon et al. 2016).

CMNF can be isolated from wood or plant cell walls with a purely mechanical or chemical treatment, or with a chemical or enzymatic pretreatment followed by a mechanical treatment (Nechyporchuk et al. 2016; Desmaisons et al. 2017). The fibrils may greatly differ in size and morphology depending on the fibrillation method used (see Fig. 1). Typically, the lateral dimension of fibrils is on the nanometer scale, and the length is up to several micrometers. The aspect ratio and the number of fibril-fibril contacts can, thus, be very high. The specific surface area (and, thus, hydroxyl group surface density) is also much higher than for regular cellulose fibers. For these reasons, CMNF suspensions can form yield-stress gels already at approx. 0.1–0.3% consistency (Varanasi et al. 2013; Raj et al. 2016; Li et al. 2016; Arola et al. 2018), and they have a strong tendency towards aggregation and flocculation (Karppinen et al. 2012; Pääkkönen et al. 2016; Hubbe et al. 2017; Raj et al. 2017; Koponen et al. 2018). Notice that highly crystalline and rigid nanoparticles having a length of a few hundred nanometers can be obtained from cellulose e.g. by acid hydrolysis (Peng et al. 2011; Hubbe et al. 2017). The resulting material, known as cellulose nanocrystals (CNC), is outside the scope of this article.

Fig. 1
figure 1

Examples of various types of CMNF fibrils: a mechanically disintegrated commercial CMNFs Celish (Turpeinen et al. 2019), b mechanically disintegrated CMNFs from birch (Turpeinen et al. 2019), c mechanically disintegrated CMNFs from softwood (Kumar et al. 2016), d purely chemically disintegrated CMNFs from sugar beet (Agoda-Tandjawa et al. 2010), e entzymatically disintegrated CMNFs from softwood (Pääkkö et al. 2007), and f TEMPO-oxidized CMNFs from eucalyptus (Honorato et al. 2015). Figures published with permission from the publishers

The rheological characteristics of various CMNF suspensions has become a widely discussed topic. Although knowledge on rheological behavior is naturally important in the use of CMNFs as a rheology modifier (Dimic-Misic et al. 2013; Shao et al. 2015; Li et al. 2015) and stabilizer (Andresen and Stenius 2007; Winuprasith and Suphantharika 2013), such information is also needed for CMNF production (Delisée et al. 2010; Pääkkönen et al. 2016; Colson et al. 2016) and for other CMNF-related processes (Saarikoski et al. 2015; Shao et al. 2015; Hoeng et al. 2017; Kumar et al. 2017). A recent thorough review of the rheology of CMNF suspensions can be found in (Hubbe et al. 2017).

As CMNF suspensions are non-Newtonian yield stress materials, the viscous (shear-thinning) behavior of aqueous CMNF suspensions is sometimes analyzed in the context of the Herschel-Bulkley equation (Mohtaschemi et al. 2014)

$$\tau = \tau_{y} + K\dot{\gamma }^{n} .$$
(1)

Above, \(\tau_{{}}\)is the shear stress, \(\tau_{y}\) is the yield stress, \(\dot{\gamma }\)is the shear rate, K is the consistency index, and n is the flow index. The use of Eq. (1), however, necessitates very accurate measurements, especially with small shear rates where, e.g., the effect of apparent slip flow on the geometry boundaries (Buscall, 2010; Saarinen et al. 2014; Haavisto et al. 2015; Martoia et al. 2015; Nazari et al. 2016) has been successfully eliminated. In practice, the rheological behavior of CMNF suspensions is usually described by the power law

$$\mu = K\dot{\gamma }^{n - 1}$$
(2)

(Lasseuguette et al. 2008; Moberg et al. 2014; Honorato et al. 2015). Note that it has been suggested that the consistency index K reflects individual fibril characteristics, whereas the flow index n reflects the structural property of the whole suspension (Tatsumi et al. 2002).

For suspensions with elongated particles shear thinning becomes more pronounced as the particle aspect ratio, particle flexibility, or consistency increases (Goto et al. 1986; Switzer and Klingenberg 2003; Bounoua et al. 2016). An explanation of shear thinning, which works well for polymers (Wagner et al. 2013), is the alignment of particles due to shear. Owing to the high entanglement of fibrils, this is unlikely a sufficient mechanism for most CMNF suspensions. The mechanism behind the shear-thinning behavior of CMNF suspensions, which is also seen with other fibrous materials (Mongruel and Cloitre 1999; Switzer and Klingenberg 2003; Derakhshandeh et al. 2011; Jiang et al. 2016), is likely to be caused by (adhesive) interactions between the fibrils. Hydrodynamic shear forces are more effective in breaking fibril-fibril contacts when the shear rate increases; this is reflected in a decrease in floc size and an increase in the orientation of the fibrils. As a result, the efficiency of momentum transport in the suspension declines and the viscosity decreases (Petrich et al. 2000; Iotti et al. 2011; Bounoua et al. 2016). This view is supported by the observation that in pulp suspensions, with a given shear rate, decreasing floc size decreases the viscosity of the suspension (Kerekes, 2006).

Factors that affect the rheological properties of standard CMNF suspensions and, thus, also the values of parameters K and n in Eqs. (1) and (2), fall into several categories. One is morphology, which includes fibril flexibility and shape, length and diameter distributions, aspect ratio, fibrillation, and network/floc structure. Another factor is the surface-chemical composition of the fibrils, which can affect surface charge and the colloidal interactions between them. Both morphology and surface composition can depend on the processes of treatments used to prepare the material. The most universal parameter that affects CMNF rheology is indisputably mass consistency. Dimensional analysis shows that for flows of spherical non-Brownian surface-chemically neutral particles, the viscosity of a suspension is a function of only two factors: consistency and the viscosity of the carrier fluid; particle size does not affect it (Mewis and Wagner 2012). Aspect ratio is another key parameter for elongated particles in this context. Two suspensions can, thus, have very similar viscous behavior even though their particle size differs considerably.

While the values of parameters K and n have been reported for CMNF suspensions in numerous studies, there are only a few studies where the dependence of the consistency index K and the flow index n on consistency have been studied systematically. In Schenker et al. (2018) and Turpeinen et al. (2019) these relations were found to be power laws, i.e. \(K \sim c^{{m_{K} }}\) and \(n \sim c^{{m_{n} }}\), where c is the mass consistency of the suspension. Note that similar scaling behavior is well known for the storage modulus G′ and the loss modulus G″ of CMNF suspensions (Tatsumi et al. 2002; Pääkkö et al. 2007; Naderi and Lindström 2014; Quennouz et al. 2016). A power law relation,\(\tau_{y} \sim c^{{m_{\tau } }}\) has been reported for yield stress in several studies (Lowys et al. 2001; Tatsumi et al. 2002; Varanasi et al. 2013; Kumar et al. 2016).

Figure 2 shows the viscosity of various CMNF suspensions as a function of consistency with a fixed shear rate of 1.0 1/s (Hubbe et al. 2017). The figure shows that there is great variability in the viscosity values caused, e.g., by different raw materials and production methods, which result in great versatility in fibril shape, flexibility, fibrillation, entanglements, and surface chemistry. (Hubbe et al. 2017) conclude that, due to this divergent nature of CMNF materials, “the relationship between measured viscosity and nanocellulose solids content is likely to be irregular”, and not necessarily conform to any of the many viscosity formulae they reviewed. Hubbe et al. (2017) also emphasizes the variability of exponent mK in the scaling law \(K \sim c^{{m_{K} }}\) in the literature. We believe that this view reflects well the current general opinion of the rheology of CMNF suspensions. However, while it is likely that no general formula can be found for predicting the viscosity of a given CMNF suspension a priori, it might still be possible to find useful general guidelines on its viscous behavior. Even in Fig. 2, which includes an extensive collection of different CMNF materials in various conditions, the general trend appears to be a power law with\(\mu \sim c^{2.3}\).

Fig. 2
figure 2

Viscosity of various aqueous CMNF suspensions as a function of consistency at the shear rate of 1.0 1/s. The data was collected from 40 studies (Hubbe et al. 2017). The aqueous solutions had water-like viscosity (not dominated by polyelectrolytes) and a wide range of ionic strengths, surface charges, and other details of experimentation. The solid line is a power law fit to the data. (To avoid the disturbing effect of outliers, which are common in the datasets analyzed here, we used the Matlab robustfit function throughout this study for fitting. The MATLAB version was 9.5.0.944444, R2018b)

In this review, we examined the dependence of shear viscosity and the yield stress of aqueous CMNF suspensions on consistency. Our focus was on the three scaling laws of\(K \sim c^{{m_{K} }}\), \(n \sim c^{{m_{n} }}\)and\(\tau_{y} \sim c^{{m_{\tau } }} .\) We included such studies where the rheology of at least three different CMNF consistencies were analyzed. Moreover, we included only studies where the measurements were performed on original CMNF materials. Studies where, e.g., ionic strength was modified or polymers were added, have been excluded. The observed similarities in the shear rheology of various CMNF suspensions suggest that the shear rheology of CMNF suspensions is more universal than has previously been realized.

Materials

Table 1 lists the included studies together with the CMNF raw material, CMNF manufacturing method, and in some cases the fibril dimensions. The aspect ratio was calculated from the middle values of the given length and width ranges. The fibril dimensions of Moberg et al. (2014) were obtained from Wågberg et al. (2008). The width of the Celish fibrils was taken from Kose et al. (2015). The fibril widths for Schenker et al. (2018) and Schenker et al. (2019) were obtained from the supplementary material of Schenker et al. (2019). Varanasi et al. (2013) used the finest 20% fraction of Celish. Notice that the rheological analysis of the measurement data presented in Kataja et al. (2017) can be found in Koponen et al. (2019).

Table 1 From left to right: reference, raw material, manufacturing method (purely mechanical/purely chemical/pretreatment before mechanical disintegration), is sonication used as post treatment for dispergation, fibril length, fibril width, and fibril aspect ratio

Table 1 shows that there is great variability in the fibril dimensions and aspect ratios among the different studies. This likely reflects both the great versatility of CMNF materials as well as experimental uncertainties related to the measurement of width and length distributions. Due to their small size and often very wide width and length distributions, measuring the average diameter, not to mention the average length, is very challenging for CMNF materials. For this reason, the fibril dimensions shown in Table 1 are most likely not commensurate.

Table 2 shows the consistency range, the number of measured consistencies, the shear rate region, the existence of yield stress analysis, and the rheometer geometry used. In some cases, the lowest consistencies are close to or even below the gel point of the CMNF suspension, see e.g. Geng et al. (2018) and Lasseuguette et al. (2008). Some studies have reported the flow index K and the power index n explicitly for every consistency. However, in many studies, this data was not given. When that occurred, a graphical interface (grabit function written by Jiro Doke for Matlab) for acquiring viscosity-shear rate data from the figures was used, and the power law in Eq. (2) was fitted to the data. As discussed in next paragraph, shear rate regions with suspicious behavior were sometimes excluded from this analysis. In such cases the original values of K and n reported by the authors were not used.

Table 2 From left to right: reference, consistency range, number of consistencies, shear rate region, is the yield stress analyzed, and the rheometer geometry used (cp—cone and plate, pp—plate plate, bc—bob and cup, r—roughened or serrated walls)

In e.g. Agoda-Tandjawa et al. (2010), the shear rate-viscosity rheograms have two shear rate regions with different slope and a plateauing region between them. In contrast, in e.g. Tatsumi et al. (2002), the shear stress-shear rate rheograms have regions where shear stress is constant. In both studies the shear stress was probably below the yield stress, and slip flow had probably taken place at the lowest shear rates (Vadodaria et al. 2018; Turpeinen et al. 2019). With higher shear rates; see e.g. (Charani et al. 2013); power law behavior was, in some cases, violated—probably due to some measurement-related instabilities. We have, whenever possible, excluded these shear rate regions from our analysis, and the shear rate region used for fitting the power law in Eq. (2) to the measurement data are then smaller than those shown in Table 2. (Nazari et al. 2016) have performed their measurements with exceptionally high shear rates. It seems that viscosity had saturated at the highest shear rates in their measurements. We have eliminated these measurement points from the analysis. In Mohtaschemi et al. (2014), Herschel-Bulkley Eq. (1) was used successfully in the rheological analysis. We have used these values for parameters K and n.

Finally note that, in a given consistency range, the values of the consistency index K can vary several orders of magnitude with different consistencies (see Fig. 3), while the variation of the flow index n is much smaller (see Fig. 5). For this reason, small discrepancies in the experimental setup (such as slip flow) do not have a big effect on the relative size of the values of K with different consistencies. The values of n, however, may, in such cases, be distorted much more, which is manifested, e.g., in the nonmonotous behavior of n as a function of consistency.

Fig. 3
figure 3

Consistency index K of Eq. (2) as a function of consistency obtained from different studies. The solid line is a power law fit to the data

Results

Shear viscosity

Figure 3 shows the consistency index K as a function of consistency for different studies. The figure shows that the consistency index usually follows a power law

$$K = K_{0} c^{{m_{K} }} ,$$
(3)

where K0 and mK are CMNF-dependent material parameters. Table 3 shows the values of K0 and mK for different CMNFs. We can see that K0 varies strongly; the difference between lowest and highest values being more than two orders of magnitude. The exponents mK, in contrast, are quite close to each other.

Table 3 Parameters K0 and mK of Eq. (3) and n0 and mn of Eq. (5) obtained with the fit to the data shown in Figs. 3 and 5. Note that, for (Lasseuguette et al. 2008) and (Geng et al. 2018), the obvious outliers (see Fig. 4) at 0.8% and 0.05%, respectively, have been omitted from the fits

A power law fit to all the data points in Fig. 3 gives \(K = 6.4 \times c^{2.75}\). However, the exponent is rather sensitive to the distribution of data points, and it does not necessarily accurately reflect the true average behavior of these materials. A better estimate is obtained by doing the analysis for the normalized consistency index \(\bar{K} = {K \mathord{\left/ {\vphantom {K {K_{0} }}} \right. \kern-0pt} {K_{0} }}\) (see Fig. 4). Then, the power law fit gives

Fig. 4
figure 4

Normalized consistency index \(\bar{K} = {K \mathord{\left/ {\vphantom {K {K_{0} }}} \right. \kern-0pt} {K_{0} }}\) as a function of consistency. The markers are the same as in Figure 3. A power law fit to the data gives \(K \sim c^{2.43}\)

$$\bar{K} = 0.98 \times c^{2.43} .$$
(4)

Figure 4 shows that while there are some outliers, all CMNFs generally follow very well the power law in Eq. (4).

Lowys et al. (2001), Lasseuguette et al. (2008) and Geng et al. (2018) have shown that below the gel point the slope of the viscosity versus consistency curve is much shallower and the flow is almost Newtonian. In Geng et al. e.g., the viscosity of the CMNF suspension scaled as \(\mu \sim c^{0.4}\) below the gel point of 0.2% (with shear rate 100 1/s). In Fig. 4, the deviation of the two low consistency data points of Geng et al. (2018) from the general trend is due to this behavior.

Figure 5 shows the flow index n as a function of consistency for different studies. We can see in the figure that there are several outliers, and in many cases, the values of n do not decrease monotonously with increasing consistency. Table 3 shows the values of n0 and mn of a power-law fit

$$n = n_{0} c^{{m_{n} }}$$
(5)

for different CMNFs. We can see that while n0 is of the same order in all cases, the variation of mn is rather high. It is likely that many of these variations are due to experimental uncertainties and not due to fundamental differences in the rheological behavior of these materials. This possibility is emphasized by the fact that the biggest deviations from the average were obtained with either smooth plate-plate or smooth cone-and-plate geometries (compare Table 2 and Table 3), which are rather sensitive geometries to slip flow. The solid line in Fig. 5 is a power law fit of Eq. (5) to all data points:

Fig. 5
figure 5

Flow index n of Eq. (2) as a function of consistency obtained from different studies. The solid line is a power law fit to the data

$$n = 0.30 \times c^{ - 0.43} .$$
(6)

This equation gives a reasonable approximation of the general rheological behavior of the flow index n of CMNF suspensions.

By combining Eqs. (4) and (6), we get the following general formula for the viscosity of CMNF suspensions:

$$\mu (c,\dot{\gamma }) = Mc^{2.43} \dot{\gamma }^{n - 1} ; \, n = 0.30c^{ - 0.43} .$$
(7)

Above M is a material parameter that depends on the CMNF grade. Table 4 shows the values of M for various CMNFs. We can see in the table that, in most cases, M is close to K0—in an ideal case they would be equal.

Table 4 Parameter M for various CMNF materials. Also shown are M/K0, and the mean and standard deviation of \(\delta\) defined in Eq. (8)

Figure 6 shows the viscosity given by Eq. (7) as a function of viscosity calculated from Eq. (2) for all the CMNFs shown in Table 4 for shear rates of 0.1, 1.0, 10, 100, and 1000 1/s. The correlation between these values is seen to be very high (R2 = 0.99). Table 4 shows the mean and standard deviation for the various CMNFs of

$$\delta \left( {\dot{\gamma },c} \right) = \frac{{\left| {\mu \left( {\dot{\gamma },c} \right)_{\text{model}} - \mu \left( {\dot{\gamma },c} \right)_{ \exp } } \right|}}{{\mu \left( {\dot{\gamma },c} \right)_{\exp } }} ,$$
(8)

where \(\mu \left( {\dot{\gamma },c} \right)_{\text{model}}\) and \(\mu \left( {\dot{\gamma },c} \right)_{ \exp }\)are the modelled and measured viscosity for shear rate \(\dot{\gamma }\) and consistency c, respectively. As we can see in Table 4, there is considerable variation in the values of mean(δ) and SD(δ). For most CMNFs, Eq. (7) gives viscosity with a good accuracy, while for some CMNFs, the accuracy is only decent. This difference can also be seen in Fig. 6, where blue and red circles show the best case (Moberg et al. 2014) and the worst case (Geng et al. 2018), respectively.

Fig. 6
figure 6

The viscosity given by Eq. (7) as a function of viscosity calculated from Eq. (2) for all the CMNFs shown in Table 4. The viscosities were calculated for shear rates of 0.1, 1.0, 10, 100, and 1000 1/s. The solid line shows a power law fit to the data. Blue and red circles show the best case (Moberg et al. 2014), R2 = 1.00, and the worst case (Geng et al. 2018), R2 = 0.84, respectively

In the literature, the consistency dependence of the viscosity of CMNF suspensions is often given in the simple form of \(\mu \sim c^{m}\) (Hubbe et al. 2017). This is, however, misleading, as parameter m is not unique but depends on the shear rate. Fig. 7 shows as an example the viscosity given by Eq. (7) for various shear rates with M  = 1. Also shown are power law fits for each shear rate. We can see in the figure that exponent m increases from 1.6 to 2.8 when the shear rate decreases from 1000 to 0.1 1/s.

Fig. 7
figure 7

Viscosity given by Eq. (7) with M=1 for various shear rates. The power law fit was performed for the consistency range of 0.3–3%

Yield stress

Figure 8 shows the yield stress \(\tau_{y}\) as a function of consistency for different studies. We can see that, just like the consistency index, the yield stress also follows a power law

$$\tau_{y} = \tau_{0} c^{{m_{\tau } }} ,$$
(9)

where τ0 and mτ are CMNF-dependent parameters. Table 5 shows the values of τ0 and mτ for different CMNFs. According to Bennington et al. (1990) τ0 depends on Young’s modulus and the aspect ratio of the fibrils. We can see in the figure that τ0 varies strongly; the difference between lowest and highest value is more than two orders of magnitude. The exponents mτ, on the other hand, are quite close to each other. Notice that there was no correlation between parameters K0 and τ0. Thus, the absolute levels of viscosity and yield stress appear to be uncorrelated. (This is also reflected in the high variability of parameter a in Table 6). It is possible that the experimental uncertainties discussed below are too high and the dataset is too small for this analysis. However, this surprising finding merits further investigation.

Fig. 8
figure 8

Yield stress as a function of consistency for different studies. The solid line shows a power law fit to the data

Table 5 Parameters τ0 and mτ of Eq. (9) obtained with a power law fit to the data shown in Fig. 8
Table 6 The fitting paremeters of Eq. (11) for the data shown in Fig. 10

A power law fit to all the data points in Fig. 8 gives \(\tau_{y} = 6.1 \times c^{2.35}\). A corresponding fit to the normalized yield stress \(\bar{\tau }_{y} = {{\tau_{y} } \mathord{\left/ {\vphantom {{\tau_{y} } {\tau_{0} }}} \right. \kern-0pt} {\tau_{0} }}\) gives (see Fig. 9)

Fig. 9
figure 9

Normalized yield stress \(\bar{\tau }_{y} = {{\tau_{y} } \mathord{\left/ {\vphantom {{\tau_{y} } {\tau_{0} }}} \right. \kern-0pt} {\tau_{0} }}\) as a function of consistency. The markers are the same as in Figure 8. A power law fit to the data gives \(\tau_{y} \sim c^{2.26}\)

$$\bar{\tau }_{y} = 1.00 \times c^{2.26} .$$
(10)

We can see in Fig. 9 that most CMNFs follow very well the power law in Eq. (10).

Relation between consistency index and yield stress

We can see from Eqs. (4) and (10) that the consistency dependence of K and \(\tau_{y}\)is almost identical. It is, thus, interesting to directly compare how these two quantities relate to each other. Figure 10 shows the consistency index as a function of the yield stress for studies where both quantities have been reported. Table 6 shows the values of parameters a and b for the power law

Fig. 10
figure 10

Consistency index K as a function of yield stress \(\tau_{y}\). The solid line shows a power law fit to the data

$$K = a\tau_{y}^{b} .$$
(11)

We can see in Table 6 that, while there is great variation in a, the values of exponent b are, in most cases, quite close to each other. Figure 11 shows the normalized consistency index \({{\bar{K}_{a} = K} \mathord{\left/ {\vphantom {{\bar{K}_{a} = K} a}} \right. \kern-0pt} a}\) as a function of yield stress. The solid line shows the power law fit

$$\bar{K}_{a} = 1.01 \times \tau_{y}^{1.12}$$
(12)

to the data. Most data points follow this formula very accurately. The relation between the consistency index and the yield stress is, consequently, generally almost linear.

Fig. 11
figure 11

Normalized consistency index \({{\bar{K}_{a} = K} \mathord{\left/ {\vphantom {{\bar{K}_{a} = K} a}} \right. \kern-0pt} a}\) as a function of yield stress. The markers are the same as in Figure 10. The solid line shows a power law fit to the data

Discussion

The rheological behavior of CMNF materials is often compared with pulp fiber suspensions. The yield stress of pulp suspensions also generally follows the power law in Eq. (9). In numerous studies, the values of mτ have varied in the range of 1.6–3.3 (Bennington et al. 1990; Huhtanen and Karvinen 2006; Derakhshandeh et al. 2010a; Derakhshandeh et al. 2011; Sumida and Fujimoto 2015). In Dalpke and Kerekes (2005), e.g., exponent mτ slowly approached the value of 2.4 (shorter fibers had a higher exponent) when the fiber length (or aspect ratio) increased, while in Bennington and Kerekes (1996), the energy dissipation needed for floc level fluidization of pulps scaled as ~ c2.5. Due to experimental difficulties, relatively few studies have been devoted to measuring the viscosity of pulp suspensions. Huhtanen and Karvinen (2006) used Eq. (3) and obtained mK=1.8, while in Bennington and Kerekes (1996), the viscosity of pulp scaled as ~ c3.1. In Derakhshandeh et al. (2010b), however, the relation between consistency index and consistency was approximately linear. Experiments with pulp suspensions have, thus, with some exceptions, given scaling exponents similar to CMNF suspensions with an even wider variation, at least partly due to higher experimental uncertainties. However, there is enough similarity between the rheological behaviour of pulp and CMNF suspensions to speculate that similar mechanisms might take place during the shearing of them. Thus, in their basic form (i.e. with low ionic strength and without extra polymers or surface modification), CMNFs might rheologically be somewhat surprisingly simply another group of fibrous materials with a high aspect ratio, not too far from other fiber suspensions.

A possible reason for the observed similarity of the rheological behavior of most CMNFs (and possibly also for the similar behavior of pulp fiber suspensions) may be the strong entanglement of the CMNF fibrils. Micrographs of CMNF products usually show highly complex structures that are, similar to pulp suspensions (Kerekes et al. 1985), better described as networks of flocs rather than individual fibrils. Typical lengths of the aggregated CMNF structures are in the range of 20–1000 μm (Hubbe et al. 2017). Consequently, although the fibrils that compose CMNFs can be clearly in the “nano” range, the gross structure is typically a lot larger. Björkman (2005) studied the dynamics of pulp suspension and concluded that many observed flow phenomena could be given more natural and comprehensible explanations by using the concept of floc instead of the traditional individual fibre view. Similarly, rather than having a suspension of individual fibrils, the rheology of a CMNF suspension might more appropriately be modelled as a suspension of flocs dispersed in a liquid phase or in a gel-like matrix (Saarikoski et al. 2012; Hubbe et al. 2017). As a consequence, the behavior of different CMNF suspensions is similar at mesoscopic and macroscopic scales. Figure 12 shows an example of the floc behavior of CMNFs during shearing. The floc structure is clearly seen in the figure; floc size decreases with increasing shear rate, but the suspension is not totally homogenized even at the highest shear rate of 500 1/s.

Fig. 12
figure 12

Mechanically disintegrated 1.0% CMNF suspension (See Fig. 1b) in a bob and cup geometry with a gap width of 1.0 mm and with shear rates of 1.0–500 1/s. The upper wall is stationary and the lower wall moves from left to rigtht. The dynamics of the system is obviously not that of individual fibrils but that of flocs. Floc size decreases with increasing shear rate. The white curves in a, b show the local axial velocity of the suspension. The suspension yields with shear rates higher than 20 1/s. So, in a, b, the shear stress is below the yield stress; there is slip flow on the walls and no shearing inside the supension. g) The local axial velocity of the suspension with the shear rates of 230, 340 and 500 1/s. The suspension is sheared but there is still strong apparent slip flow on the walls. Notice that the rheogram of this CMNF is shown in Fig. 13. Figures published with permission from the publisher (Haavisto et al. 2015)

Equation (3), with mK = 2.43, and Eq. (9), with mτ = 2.26, are a useful starting point when analyzing the measured rheological data. If the obtained fitting values of K00 and τ00 differ significantly from the values K0 and τ0 (obtained with a fit that also has mK and mτ as fitting parameters), one should consider if the difference is real or due to experimental uncertainties. For the studies reviewed in the present study, K00 and τ00 were, in most cases, close to K0 and τ0. The two exceptions were (Nazari et al. 2016) where K00/K0 = 0.84 and τ000 = 1.73, and (Quennouz et al. 2016), where K00/K0 = 0.51 and τ000 = 1.19. As there is a clear difference only in either K00 or τ00, it is likely that the discrepancy is due to experimental uncertainty and not due to exceptional rheology. Also, when comparing the general rheological behavior of different CMNFs, it might be useful to compare the values of K00 and τ00 instead of K0 and τ0. If one wants to include shear thinning behavior into the analysis, Eq. (7), or its generalized version

$$\mu (c,\dot{\gamma }) = Wc^{m} \dot{\gamma }^{n - 1} ; \, n = qc^{ - p}$$
(13)

can be used (W, m, q and p above are fitting parameters). Due to high degrees of freedom, Eq. (13) usually fits very well (R2 close to one) to the viscosity-shear rate data. Eq. (13) is, consequently, preferred if a good quantitative description of the data is required. However, the use of Eq. (7) should always be considered as it may give a better description of the real rheological behavior, especially outside the measurement range. Moreover, the single free fitting parameter M can then be used for comparing the general rheological behavior of different CMNFs.

Due to flocculation (see Fig. 12a–f) and apparent slip flow (see Fig. 12a, b, g), the shear viscosity of CMNF suspensions is exceptionally difficult to measure quantitatively with good accuracy (Saarinen et al. 2014; Haavisto et al. 2015; Martoia et al. 2015; Nazari et al. 2016). Apparent slip flow may have influenced viscosity data considerably also when it is not apparent in the rheogram. Figure 13 shows an example of two rheograms obtained with a bob and cup rheometer for the mechanically disintegrated CMNFs seen in Figs. 1b and in Fig. 12. We can see that even though the curves look reasonable, there are only a few legitimate measurement points that reflect the true viscous behavior of the suspensions with reasonable accuracy. When the shear rate is below 20 1/s, shear stress is below the yield stress of the suspensions. The apparent flow is in this case only due to slip flow at the solid walls. Moreover, even above the yield stress, the slip flow introduces high errors to the measured viscosity.

Fig. 13
figure 13

Rheograms for a mechanically disintegrated CMNF (See also Figs. 1b and 12) for consistencies of 0.5% (spheres) and 1.0% (triangles) obtained with a bob and cup rheometer with a 1 mm gap. Solid lines are fits of Eq. (2) to the measurement points (obtained fitting values of n are also shown). There are only a couple of measurement points that reflect the true rheological behavior of the CMNF suspensions with reasonable accuracy. Most points are below the yield stress, and above the yield stress errors in viscosity values are in most cases high. Notice that Fig. 12 shows the floc structure of this CMNF in the rheometer. Figure published with permission from the publisher. (Turpeinen et al. 2019)

Schenker et al. (2018) have studied shear viscosity and yield stress with different rheometer geometries. The difference between the measured viscosity and yield stress was found to be, in the worst cases, an order of magnitude. According to Vadodaria et al. (2018), scientific literature on the rheological properties of CMNFs is due to slip replete with such erroneous viscosity data. Note that for yield stress, a unique definition does not even exist, and different measurement techniques may give different results even though the measurements have been performed rigorously. Nazari et al. (2016) have studied the yield stress of a CMNF suspension with various methods on rotational rheometers (in the present article, we used their values obtained with a vane geometry) and found that the obtained yield stress values could vary as much as a factor of 4. Table 7 shows the fitting parameters of Eqs. (3), (5), and (9) for different rheometer geometries for the viscosity and yield stress data of Schenker et al. (2018). We can see in the table that, although the parameters K0 and τ0 vary significantly between the different rheometer geometries, parameters mK and mτ are rather close to each other. The variation of parameters n0 and mn is also rather small. Also notice that in Fig. 13, despite of relatively high quantitative errors in viscosity, the values of exponent n of power law Eq. (2) are with higher shear rates well in line with those shown in Fig. 5. While the reason behind such behavior warrants further study, it explains why the scaling behavior of shear viscosity and yield stress on consistency is often rather similar among different studies even though the measurement data may have in some cases been quantitatively wrong due to slip flow.

Table 7 Fitting parameters of Eqs. (3), (5), and (9) for different rheometer geometries for mechanically disintegrated CMNFs in the consistency range of 0.52.0% (bcbob and cup, r—roughened or serrated walls, the gap size is in millimeters). The obtained power law parameters were found to grow/decrease rather systematically with improved rheometer geometry. Red and blue arrows show the direction of the increasing and decreasing tendency, respectively (Schenker et al. 2018)

Rough or serrated walls are often used to eliminate the slip at walls, but they may not solve the problem entirely for basic rotational geometries [see Table 7 and (Nechyporchuk et al. 2014)]. However, vane-in-cup geometry, especially with a wide gap, has been found to decrease slip effects for both wood fiber suspensions (Mosse et al. 2012) and CMNF suspensions (Mohtaschemi et al. 2014). A wide gap, however, is not always problem-free. It may introduce heterogeneous flow in the vane geometry, which must be properly addressed in order to obtain correct results. A wide gap may also cause the system to be susceptible to secondary flows (Mohtaschemi et al. 2014). Finally notice that the effect of slip flow on the rheological analysis can be totally eliminated by combining the rheometer with a velocity profiling technique, such as Ultrasound Velocity Profiling (UVP), Nuclear Magnetic Resonance Imaging (MRI), or Optical Coherence Tomography (OCT) (Haavisto et al. 2017; Koponen and Haavisto 2018). Of the studies reviewed here, this approach was used in Kataja et al. (2017) and Turpeinen et al. (2019). Notice that the scaling parameters obtained in these two studies are very similar to those obtained in the other studies (see Table 3). So, while slip flow is likely involved in many of the rheological studies reviewed here, it hasn’t dominated the scaling behaviour of the CMNF suspensions.

Conclusions

Although it is obvious that CMNF raw material and production method strongly affect the absolute values of yield stress and the viscosity of CMNF suspensions, it seems that the scaling laws on consistency are similar for a wider group of CMNF materials than previously believed. The consistency index and yield stress of CMNF suspensions reviewed in this study generally followed the scaling laws \(K \sim c^{2.43}\) and\(\tau_{y} \sim c^{2.26}\), and the relation between the consistency index and the yield stress was almost linear. While the measured values of the flow index n varied—possibly mostly due to experimental uncertainties, a decent approximation for n was given by the formula\(n = 0.30 \times c^{ - 0.43}\).

The variability of the scaling law parameters of CMNF suspensions found in the literature (see Tables 3 and 5) is probably not so much due to real differences in the physical behavior of the suspensions rather than due to experimental uncertainties and to general difficulties in measuring the rheological behavior of these suspensions rigorously. Here the biggest source of uncertainly is apparent slip flow that may introduce high quantitative errors in the measured viscosities. However, it seems somewhat surprisingly that the general scaling behaviour of these materials is rather intact to slip flow. There is no obvious explanation for this observed behaviour—this problem definitely warrants further study.

The reason behind the universal rheological behavior of CMNF at mesoscopic and macroscopic scales might be the strong entanglement of fibrils; the flow dynamics of typical CMNF suspensions is dominated by interactions between fibril flocs and not so much by interactions between individual fibrils.

The obtained scaling laws were used to form a general formula, Eq. (7), for CMNF viscosity as a function of shear rate and consistency. Although this formula has only one fitting parameter, it worked quite well for most CMNFs reviewed in this study. This formula can be a useful framework when interpreting and analyzing measured CMNF rheological data.

In the future, a corresponding analysis should be performed in the LVE-region for parameters of oscillatory rheology, i.e. for the storage modulus G’ and the loss modulus G″, and their relation to the consistency index K and the yield stress τy.