Incorporating Rheological Nonlinearity into Fractional Calculus Descriptions of Fractal Matter and Multi-Scale Complex Fluids

Abstract

In this paper, we use ideas from fractional calculus to study the rheological response of soft materials under steady-shearing flow conditions. The linear viscoelastic properties of many multi-scale complex fluids exhibit a power-law behavior that spans over many orders of magnitude in time or frequency, and we can accurately describe this linear viscoelastic rheology using fractional constitutive models. By measuring the non-linear response during large step strain deformations, we also demonstrate that this class of soft materials often follows a time-strain separability principle, which enables us to characterize their nonlinear response through an experimentally determined damping function. To model the nonlinear response of these materials, we incorporate the damping function with the integral formulation of a fractional viscoelastic constitutive model and develop an analytical framework that enables the calculation of material properties such as the rate-dependent shear viscosity measured in steady-state shearing flows. We focus on a general subclass of fractional constitutive equations, known as the Fractional Maxwell Model, and consider several different analytical forms for the damping function. Through analytical and computational evaluations of the shear viscosity, we show that for sufficiently strong damping functions, for example, an exponential decay of fluid memory with strain, the observed shear-thinning behavior follows a power-law response with exponents that are set by the power-law indices of the linear fractional model. For weak damping functions, however, the power-law index of the high shear rate viscosity is set by the terminal behavior of the damping function itself at large strains. In the limit of a very weak damping function, the theoretical formulation predicts an unbounded growth of the shear stress with time and a continuously growing transient viscosity function that does not converge to a meaningful steady-state value. By determining the leading terms in our analytical solution for the viscosity at both low and high shear rates, we construct an approximate analytic expression for the rate-dependent viscosity. An error analysis shows that, for each of the damping functions considered, this closed-form expression is accurate over a wide range of shear rates.

1. Introduction

A wide range of soft solids and complex fluids from consumer products to fracking fluids exhibit viscoelastic properties characterized by a broad spectrum of relaxation time scales. This can make the quantitative constitutive modeling of the material response very challenging. At its simplest level, the incorporation of both solid-like and liquid-like responses for a viscoelastic material can be achieved by using a simple mechanical combination of a linear elastic spring and a constant viscosity dashpot [1]. The serial combination of a dashpot and a spring is commonly referred to as the Maxwell model [1,2] and, in this case, the constitutive relationship between the shear stress and strain can be written as:

$$\sigma +\frac{{\eta }_0}{G_0}\frac{d\sigma }{dt}={\eta }_0\frac{d\gamma }{dt},$$

(1)

where $\sigma$ is the shear stress, $\gamma$ denotes the strain, ${\eta }_0$ is the viscosity and $G_0$ is the elastic modulus of the material. The time-dependent relaxation modulus $G\left(t\right)$ for the linear Maxwell model can be derived by finding the stress response to an imposed step strain function:

$$G\left(t\right)=G_{0}exp(-t/\tau ),$$

(2)

where $\tau ={\eta }_0/G_0$ is the single characteristic relaxation time of the Maxwell model. The complex modulus $G^{*}\left(\omega \right)=G^{\prime }\left(\omega \right)+i{G}^{″}\left(\omega \right)$ characterizing the response of the Maxwell model to a small amplitude oscillatory shear (SAOS) deformation can be obtained by taking a Fourier transform of Equation (1):

$$G^{\prime }\left(\omega \right)=G_0\frac{{\left(\omega \tau \right)}^2}{1+{\left(\omega \tau \right)}^2}$$

(3)

$$G^{″}\left(\omega \right)=G_0\frac{\left(\omega \tau \right)}{1+{\left(\omega \tau \right)}^2}$$

(4)

where the dimensionless product $\omega \tau$ is often referred to as the Deborah number. Several classes of complex materials have been shown to closely approximate Maxwell behavior with a single relaxation time scale, at least in the limit of small deformations. In particular, complex fluids constructed from ideal reversible networks or self assembling worm-like surfactants often exhibit a single dominant relaxation time scale and show Maxwell-like behavior for small deformation amplitudes [3,4,5].

However, most complex fluids, such as soft glassy materials, colloids, polymer melts, biopolymers, hydrogels, and microgel dispersions exhibit a wide range of relaxation timescales due to the broad range of length scales that characterize their microstructure [6,7,8]. This results in very broad relaxation spectra that cannot be described well using single relaxation models such as the classical Maxwell or Kelvin–Voigt representation [9]. Such multiscale complex materials are often observed to exhibit power-law–like responses in standard rheological experiments. To describe this rheological behavior, an additive combination of several mechanical elements can be used, resulting in multiple relaxation timescales; however, the number of discrete modes required for a quantitative description of the material response often becomes excessive [1,6]. Thus, Gemant [10,11] and Scott Blair [12] proposed to write the stress response to an imposed deformation in terms of a fractional derivative of the strain to accurately and compactly capture the power-law response observed in real industrial materials. The resulting material properties characterizing the strength of the material response (denoted as $\mathbb{G}$ and/or $\mathbb{V}$ in the present work) are, thus, quasi-properties with non-integer physical dimensions [13,14]. This class of material response can be represented in terms of the response of a ‘spring-pot’ mechanical element also commonly referred to as the Scott Blair element (shown in Figure 1b) [14]. As shown in Figure 1a, the Scott Blair element compactly describes the bulk rheological behavior of a class of viscoelastic materials known as critical gels [15], whose linear viscoelastic properties are characterized by power-law responses of the form $G\left(t\right)\sim t^{-\alpha }$ and $G^{*}\left(\omega \right)\sim {\omega }^{\alpha }$.

Researchers have additively combined spring-pot elements in series and parallel in order to construct fractional constitutive models to model more complex rheological response [14,16,17]. One such combination that is widely used in the literature is the Fractional Maxwell Model (FMM), shown in Figure 1c, which consists of a linear combination of two spring-pots in series [18]. This fractional framework is broadly applicable, since most complex materials exhibit distinct power-law behaviors at short and long time scales, which can be well described by the FMM [6,18,19,20,21,22,23,24,25,26,27].

Although these linear fractional models can accurately describe the linear viscoelastic response of a wide range of materials, the results are independent of the applied deformation amplitude and, in consequence, fail to capture the non-linear response that arises at large deformations [28]. Real soft solids and complex fluids typically undergo softening or stiffening at large deformations which affects the relaxation time scales and the magnitudes of rheological properties such as the relaxation modulus [1,29,30,31]. Thus, for example, the rate-dependent steady shear viscosity of a real viscoelastic fluid cannot be described by typical linear constitutive models formulated as fractional differential equations of any order in time.

There are many integral constitutive models that have been developed over the years in order to characterize the non-linear response of complex fluids [32,33,34,35]. One such integral formulation that has been used ubiquitously to describe non-linear viscoelastic properties is the Wagner model [30,36]. It is a special limit of the more general Rivlin–Sawyers constitutive equation [37] but ignores any possible contribution of the Cauchy–Green strain tensor to the strain history and, consequently, predicts that the second normal stress difference is zero [1]. This constitutive model has been carefully and extensively investigated and it has been shown to describe the rate-dependent shear viscosity, large amplitude oscillatory shear response, orthogonal superposition modulus as well as other viscoelastic properties of complex fluids that change with deformation or flow [21,38,39].

In the absence of non-linear effects, any general linear viscoelastic constitutive model can be written in an integral form:

$$\sigma \left(t\right)={\int }_{-\infty }^{t}G(t-t^{\prime })\dot{\gamma }\left(t^{\prime }\right)d{t}^{\prime }$$

(5)

By a change of variables and integration by parts, we can rewrite Equation (5) into:

$$\sigma \left(t\right)=-{\int }_{-\infty }^{t}M(t-t^{\prime })\gamma (t,t^{\prime })d{t}^{\prime }$$

(6)

where $M(t-t^{\prime })\equiv \partial G(t-t^{\prime })/\partial t^{\prime }$ is known as the memory function and $\gamma (t,t^{\prime })$ is the strain between time t and $t^{\prime }$. The integral model in Equation (6), when extended to non-linear deformations by combining the memory function (with a form that remains to be specified) together with a strain-dependent damping function in its integral to incorporate the softening or stiffening observed in a real viscoelastic material at large strains, is known as the Wagner integral model [40]. There have been several mathematical expressions developed for damping functions to model the material softening that arises from large deformations. Jaishankar and McKinley [6,21] selected a memory kernel consistent with the fractional integral formulation of the FMM and combined it with a simple damping function in order to numerically compute the Wagner integral and predict the rate-dependent steady shear viscosity for multiscale materials. This framework of using a fractional relaxation kernel combined with a damping function to capture the linear and non-linear response, respectively, through a Wagner model or other Rivlin–Sawyers- type formulation, is becoming increasingly common [21,29,41]. However, evaluating the rate-dependent viscosity requires a numerical computation of the Mittag–Leffler function (MLF) and careful integration over the entire strain history and, consequently, there is no closed-form analytical expression for the steady shear viscosity as a function of shear rate in the literature. In addition, there is little systematic understanding on how the steady shear viscosity for fractional models depends on the fractional model parameters that describe the linear properties or on the damping function parameters that describe the softening (or stiffening) of the material at large non-linear deformations.

In this paper, we begin by briefly summarizing the Fractional Maxwell model and highlighting two important distinguishing limits: the Fractional Maxwell Liquid (FML) and the Fractional Maxwell Gel (FMG). In Figure 1c, we provide a graphical representation for these two models. Using experimental data from three biopolymer systems (xanthan gum, casein, and alginate), we show how well these simple fractional models capture the linear viscoelastic response over a wide range of time and frequency scales. We introduce the Wagner integral model and motivate the functional form of the damping functions through measurements of the strain-dependent relaxation modulus that is obtained from step-strain experiments with xanthan gum. In [21], Jaishankar and McKinley utilize a Doi–Edwards-type damping form [42] with a linear FMM kernel and obtain the high shear-rate asymptotes for the steady shear viscosity through the Wagner integral model. Here, we derive analytical expressions for the rate-dependent viscosity predicted by the Wagner model with a fractional Maxwell kernel and different forms of the damping function: the exponential damping function [30] and the Soskey–Winter damping function [43], which is a more generalized version of the Doi–Edwards damping function considered by [21]. Through the analytical expressions we discuss the dependence of the rate-dependent viscosity on the parameters characterizing the linear fractional kernel and the nonlinear strain-dependent damping function. We observe that, for a sufficiently strong strain damping, the power-law exponents of the steady shear viscosity response are set by the power-law exponents of the linear rheological response (e.g., by the complex modulus measured in Small Amplitude Oscillatory Strain (SAOS), or the relaxation modulus determined in step strain). The material parameters in the damping function simply affect the numerical values of the pre-factors. However, for sufficiently weak strain damping, we find that the power-law exponents characterizing the rate-dependent viscosity are set by the functional form of the damping function at large shear strains. Finally, we propose a simple compact analytic expression for the shear-thinning viscosity function of a complex fluid based on the asymptotic expressions that we obtain for a generalized damping function under weak and strong shear flow conditions.

Figure 1. (a) Relaxation modulus of a 3.5 wt.% alginate gel and a casein milk protein gel (4 wt.% casein gel acidified with 1 wt.% glucono-$\delta$-lactone (GDL)). Power-law fits (Equation (11)) to the experimental data show that the materials closely obey the constitutive response predicted by a single spring-pot element or Scott Blair element, Equation (7). (b) Mechanical representation of the classical linear Maxwell model and the spring-pot element of Koeller [14] as a compact representation of a hierarchical ladder model constructed from recursive Maxwell elements [44]. (c) Mechanical representation of (i) Fractional Maxwell Gel (FMG), (ii) Fractional Maxwell Model (FMM), and (iii) Fractional Maxwell Liquid (FML). (d) The general form of the FMM simplifies to a single spring-pot when either of the quasi-properties of one mechanical element ($\mathbb{V}$ or $\mathbb{G}$) diverge to infinity.

Figure 1. (a) Relaxation modulus of a 3.5 wt.% alginate gel and a casein milk protein gel (4 wt.% casein gel acidified with 1 wt.% glucono-$\delta$-lactone (GDL)). Power-law fits (Equation (11)) to the experimental data show that the materials closely obey the constitutive response predicted by a single spring-pot element or Scott Blair element, Equation (7). (b) Mechanical representation of the classical linear Maxwell model and the spring-pot element of Koeller [14] as a compact representation of a hierarchical ladder model constructed from recursive Maxwell elements [44]. (c) Mechanical representation of (i) Fractional Maxwell Gel (FMG), (ii) Fractional Maxwell Model (FMM), and (iii) Fractional Maxwell Liquid (FML). (d) The general form of the FMM simplifies to a single spring-pot when either of the quasi-properties of one mechanical element ($\mathbb{V}$ or $\mathbb{G}$) diverge to infinity.

2. Fractional Maxwell Model (FMM) for Linear Rheology

The Fractional Maxwell Model consists of two spring-pots, each characterized by a pair of material parameters (here, denoted by $\mathbb{G},\beta$ and $\mathbb{V},\alpha$, respectively), arranged in series as illustrated in Figure 1c (ii). It can be shown that the spring-pot is a compact representation of a recursive ladder model constructed from Maxwell elements [26,44] as indicated schematically in Figure 1b: (ii) and (iii). The shear stress response of a single spring-pot is characterized by a quasi-property $\mathbb{G}$ (with units of $\mathrm{Pa}$$·s^{\beta }$). Schiessel [26] writes Equation (7) in terms of three parameters; a modulus E, a time scale $\tau$, and a fractional exponent $0\le \beta \le 1$. The stress is, thus, written as $\sigma =E{\tau }^{\beta }\left(d^{\beta }\gamma /d{t}^{\beta }\right)$; however, only the combination of ($E{\tau }^{\beta }$) and $\beta$ can actually be measured experimentally. This response can be written in terms of a fractional derivative as [14,45]:

$$\sigma \left(t\right)=\mathbb{G}\frac{d^{\beta }\gamma \left(t\right)}{d{t}^{\beta }}$$

(7)

where the fractional derivative we denoted $d^{\beta }\gamma /d{t}^{\beta }$ represents the compact form of a Caputo derivative when $\gamma \left(t\right)=0$ for $t\le 0$, and is defined as [6,46,47]:

$$\frac{d^{\beta }\gamma \left(t\right)}{d{t}^{\beta }}=\frac{1}{\mathsf{\Gamma }(1-\beta )}{\int }_0^t{(t-t^{\prime })}^{-\beta }\dot{\gamma }\left(t^{\prime }\right)d{t}^{\prime }$$

(8)

where $0\le \beta \le 1$, and $\mathsf{\Gamma }\left(a\right)={\int }_0^{\infty }{x}^{a-1}{e}^{-x}dx$ is the complete Gamma function. The overdot indicates the conventional first order derivative with time, so that $\dot{\gamma }\left(t^{\prime }\right)=d\gamma \left(t^{\prime }\right)/d{t}^{\prime }$ is the shear rate. The fractional constitutive equation in Equation (7) can, thus, also be written in an integral form:

$$\sigma \left(t\right)=\mathbb{G}{\int }_0^t{(t-t^{\prime })}^{-\beta }\dot{\gamma }\left(t^{\prime }\right)d{t}^{\prime }$$

(9)

for initial value problems in which $\dot{\gamma }\left(t^{\prime }\right)=0$$\forall t^{\prime }<0$. We also note that Equation (9) can be integrated by parts to represent the stress at time t in terms of the history of the strain field. Defining $\gamma (t,t^{\prime })={\int }_t^{t^{\prime }}\dot{\gamma }\left(t^{{}^{″}}\right)d{t}^{{}^{″}}$, we obtain:

$$\sigma \left(t\right)=-\mathbb{G}\beta {\int }_0^t{(t-t^{\prime })}^{-\beta -1}\gamma (t,t^{\prime })d{t}^{\prime }$$

(10)

All three expressions given by Equations (8)–(10) are equivalent forms of the linear viscoelastic response predicted by a single spring-pot. The relaxation modulus of the spring-pot element following a step strain $\gamma \left(t\right)={\gamma }_0\mathcal{H}\left(t\right)$ results in a decaying power-law response with time:

$$\frac{\sigma \left(t\right)}{{\gamma }_0}\equiv G\left(t\right)=\frac{\mathbb{G}}{\mathsf{\Gamma }(1-\beta )}{t}^{-\beta }$$

(11)

where $\mathbb{G}/\mathsf{\Gamma }(1-\beta )$ is known as the strength of the critical gel (often denoted by S with units of $\mathrm{Pa}·s^{\beta }$) and $\beta$ is the relaxation exponent. Boltzmann considered a power-law decaying kernel of the form of Equation (10), as he developed his original theories of viscoelasticity (see Markovitz [48]) This is exactly the form of the material response observed in a particular class of soft matter known as a critical gel [15]. As we showed in Figure 1, the simple two parameter model provides a quantitative description of the relaxation modulus of two different biopolymer gels over many decades of elapsed time. The complex moduli of a single spring-pot in a small amplitude oscillatory deformation can be computed by Fourier transformation of Equation (11) and both the elastic storage modulus and viscous loss modulus show power-law dependencies with frequency:

$$G^{\prime }\left(\omega \right)=\mathbb{G}cos\left(\frac{\pi \beta }{2}\right){\omega }^{\beta }$$

(12)

$$G^{″}\left(\omega \right)=\mathbb{G}sin\left(\frac{\pi \beta }{2}\right){\omega }^{\beta }$$

(13)

A much broader class of soft solids and complex fluids are not quite critical gels with a pure power-law response, but can be described by more complex fractional constitutive models constructed by linear combinations of the spring-pot element (see, for example, [49]). The constitutive equation relating the stress and strain for the FMM, which comprises of a series of two spring-pots arranged in series as shown in Figure 1c (ii), is described by two quasi-properties ($\mathbb{G}$, $\mathbb{V}$) and two power-law exponents ($\alpha$, $\beta$) [6,21]. This model can be written in the form:

$$\sigma \left(t\right)+\frac{\mathbb{V}}{\mathbb{G}}\frac{d^{\alpha -\beta }\sigma \left(t\right)}{d{t}^{\alpha -\beta }}=\mathbb{V}\frac{d^{\alpha }\gamma \left(t\right)}{d{t}^{\alpha }}$$

(14)

The parameters $\alpha$ and $\beta$ need to obey the inequality $0\le \beta \le \alpha \le 1$ to be thermodynamically consistent [50]. The relaxation modulus of the FMM can be obtained by calculating the response of Equation (14) for a step strain input ${\gamma }_0\mathcal{H}\left(t\right)$ [18]:

$$\frac{\sigma \left(t\right)}{{\gamma }_0}\equiv G\left(t\right)=\mathbb{G}{t}^{-\beta }{E}_{\alpha -\beta ,1-\beta }\left(-\frac{\mathbb{G}}{\mathbb{V}}{t}^{\alpha -\beta }\right)$$

(15)

where $E_{a,b}$ is the Mittag–Leffler function (MLF) [51] which is defined as:

$$E_{a,b}\left(z\right)=\sum _{k=0}^{\infty }\frac{z^k}{\mathsf{\Gamma }(ak+b)}$$

(16)

A careful consideration of the physical dimensions of Equations (14) and (15) shows that the characteristic relaxation time of the FMM can be defined as [21]:

$${\tau }_c={\left(\frac{\mathbb{V}}{\mathbb{G}}\right)}^{1/(\alpha -\beta )}$$

(17)

This characteristic time corresponds to the intersection of the two asymptotic responses of the relaxation modulus in Equation (15) at long and short times. The relaxation modulus in Equation (15) can be non-dimensionalized by defining a characteristic modulus $G_c$ in terms of the quasi-property $\mathbb{V}$ and characteristic relaxation time ${\tau }_c$ as $G_c=\mathbb{V}{\tau }_c^{-\alpha }\equiv {({\mathbb{G}}^{\alpha }/{\mathbb{V}}^{\beta })}^{1/(\alpha -\beta )}$. Consequently, the relaxation modulus can be written in a dimensionless form as:

$$\frac{G}{G_c}={\left(t/{\tau }_c\right)}^{-\beta }{E}_{\alpha -\beta ,1-\beta }\left(-{\left(t/{\tau }_c\right)}^{\alpha -\beta }\right)$$

(18)

The Fourier transform may be taken of Equation (14) to obtain the complex modulus:

$$G^{*}\left(\omega \right)=\frac{\mathbb{G}{\left(i\omega \right)}^{\beta }·\mathbb{V}{\left(i\omega \right)}^{\alpha }}{\mathbb{G}{\left(i\omega \right)}^{\beta }+\mathbb{V}{\left(i\omega \right)}^{\alpha }}$$

(19)

Similar to the relaxation modulus in Equation (18), the complex modulus can be non-dimensionalized using $G_c$ and ${\tau }_c$ and thus rewritten in a simpler, cleaner form:

$$\frac{G^{*}\left(\omega \right)}{G_c}=\frac{{\left(i\omega {\tau }_c\right)}^{\alpha }}{1+{\left(i\omega {\tau }_c\right)}^{\alpha -\beta }}$$

(20)

The real and imaginary part of the complex modulus provides us with the storage and loss modulus of the complex fluid, respectively:

$$\frac{G^{\prime }\left(\omega \right)}{G_c}=\frac{{\left(\omega {\tau }_c\right)}^{\alpha }cos(\pi \alpha /2)+{\left(\omega {\tau }_c\right)}^{2\alpha -\beta }cos(\pi \beta /2)}{1+{\left(\omega {\tau }_c\right)}^{\alpha -\beta }cos\left(\pi (\alpha -\beta )/2\right)+{\left(\omega {\tau }_c\right)}^{2(\alpha -\beta )}}$$

(21)

$$\frac{G^{″}\left(\omega \right)}{G_c}=\frac{{\left(\omega {\tau }_c\right)}^{\alpha }sin(\pi \alpha /2)+{\left(\omega {\tau }_c\right)}^{2\alpha -\beta }sin(\pi \beta /2)}{1+{\left(\omega {\tau }_c\right)}^{\alpha -\beta }cos\left(\pi (\alpha -\beta )/2\right)+{\left(\omega {\tau }_c\right)}^{2(\alpha -\beta )}}$$

(22)

The FMM has two important limiting cases: the Fractional Maxwell Liquid (FML) and the Fractional Maxwell Gel (FMG). In the FML limit, exponent $\alpha$ is set to unity ($\alpha =1$). By considering the low frequency limit of Equation (22) and the definition $G^{″}\equiv {\eta }^{\prime }\omega$, Jaishankar and McKinley [21] show that a bounded steady shear viscosity of the general Fractional Maxwell Model only exists in the FML limit (i.e., $\alpha =1$). Consequently, the three-parameter FML model ($\mathbb{G}$, $\beta$ and $\mathbb{V}$) has been used extensively to describe complex fluids in the pre-gel state [24]. It can also be shown that the constant shear viscosity of the FML model corresponds to the quasi-property $\mathbb{V}$ in the FML limit corresponding to $\alpha =1$ ($\mathbb{V}={\eta }_0$). The mechanical representation of the FML is illustrated in Figure 1c (iii), and comprises of a series arrangement of a dashpot and a spring-pot.

Secondly, in the Fractional Maxwell Gel (FMG) limit, exponent $\beta$ is constrained to be zero ($\beta =0$). This FMG model accurately captures the predominantly elastic behavior of viscoelastic materials beyond the gel point [52], and the quasi-property $\mathbb{G}$ characterizes the plateau modulus of the gel ($G^{\prime }(\omega \to \infty )=\mathbb{G}$ (Pa)). The three-parameter FMG model ($\mathbb{G}$, $\alpha$ and $\mathbb{V}$) mechanically translates to a series combination of a spring and a spring-pot as shown in Figure 1c (i). Finally, for completeness, we note that the special case of $\alpha =1$ and $\beta =0$ leads to the simple Maxwell model introduced in Equation (1) and the quasi-properties $\mathbb{V}$ and $\mathbb{G}$ simplify to the viscosity ${\eta }_0$ and $G_0$, respectively. As indicated schematically in Figure 1, the FMM also reduces to a simple single spring-pot element or Scott Blair element when either one of the quasi-properties diverges to infinity ($\mathbb{V}\to \infty$ or $\mathbb{G}\to \infty$). The asymptotes of the complex moduli for these three canonical models can be obtained from Equations (21) and (22) are tabulated in

Table 1. The asymptotes of the storage ($G^{\prime }$) and loss modulus ($G^{″}$) for the Fractional Maxwell Model (FMM) ($0<\beta <\alpha <1$), the Fractional Maxwell Liquid (FML) ($\alpha =1)$, and the Fractional Maxwell Gel (FMG) ($\beta =0$).

		$\mathit{\omega }{\mathit{\tau }}_{\mathit{c}}\ll 1$	$\mathit{\omega }{\mathit{\tau }}_{\mathit{c}}\gg 1$
FMM	$G^{\prime }/G_c$	${\left(\omega {\tau }_c\right)}^{\alpha }cos\left(\pi \alpha /2\right)$	${\left(\omega {\tau }_c\right)}^{\beta }cos\left(\pi \beta /2\right)$
	$G^{″}/G_c$	${\left(\omega {\tau }_c\right)}^{\alpha }sin\left(\pi \alpha /2\right)$	${\left(\omega {\tau }_c\right)}^{\beta }sin\left(\pi \beta /2\right)$
FML	$G^{\prime }/G_c$	${\left(\omega {\tau }_c\right)}^{2-\beta }cos\left(\pi \beta /2\right)$	${\left(\omega {\tau }_c\right)}^{\beta }cos\left(\pi \beta /2\right)$
	$G^{″}/G_c$	$\omega {\tau }_c$	${\left(\omega {\tau }_c\right)}^{\beta }sin\left(\pi \beta /2\right)$
FMG	$G^{\prime }/G_c$	${\left(\omega {\tau }_c\right)}^{\alpha }cos\left(\pi \alpha /2\right)$	1
	$G^{″}/G_c$	${\left(\omega {\tau }_c\right)}^{\alpha }sin\left(\pi \alpha /2\right)$	${\left(\omega {\tau }_c\right)}^{-\alpha }sin\left(\pi \alpha /2\right)$

.

Using the expressions given above in Equations (21) and (22), the magnitude of the complex viscosity defined as $|{\eta }^{*}\left(\omega \right)|=\sqrt{G^{\prime }{\left(\omega \right)}^2+G^{″}{\left(\omega \right)}^2}/\omega$ can be written in terms of the four FMM parameters as:

$$\frac{|{\eta }^{*}\left(\omega \right)|}{\mathbb{V}{\tau }_c^{1-\alpha }}=\frac{1}{\sqrt{{\left(\omega {\tau }_c\right)}^{2-2\alpha }+2{\left(\omega {\tau }_c\right)}^{2-\alpha -\beta }cos\left(\pi (\alpha -\beta )/2\right)+{\left(\omega {\tau }_c\right)}^{2-2\beta }}}$$

(23)

The dependence of the relaxation modulus on dimensionless time and the variation of the complex moduli and complex viscosity with the Deborah number ($\omega {\tau }_c$) for the FMM, FML, and FMG models are illustrated graphically in Figure 2.

3. Wagner Formulation for Steady Shear Viscosity

The fractional models outlined above can quantitatively describe the linear rheological response of a range of complex materials (see examples in [25]). However, they fail to capture the non-linear rheological responses measured at large deformations, such as the strain-dependent relaxation modulus or the rate-dependent viscosity [28]. To describe the non-linear response of complex materials, we make use of the Wagner integral formulation [53], in which the tensorial stress can be expressed in terms of an integral over the strain history experienced by the material, which is described by the Finger strain tensor ${\mathit{C}}^{-1}$ [30,36,40,54]:

$$\sigma \left(t\right)={\int }_{-\infty }^{t}M(t-t^{\prime })h(I_1,I_2)\left[{\mathit{C}}^{-1}(t,t^{\prime })-\mathit{I}\right]d{t}^{\prime }$$

(24)

where $I_1$ and $I_2$ are the first and second invariants of the Finger strain tensor (see [1] and [21] for more details). Here, we focus only on shear flows for which $I_1=I_2={\gamma }^2+3$, and the shear stress can, thus, be written in the form:

$$\sigma \left(t\right)=-{\int }_{-\infty }^{t}M(t-t^{\prime })h\left(\gamma \right)\gamma (t,t^{\prime })d{t}^{\prime }$$

(25)

where the minus sign in front of the integral appears because we define the finite strain as $\gamma (t,t^{\prime })=\gamma \left(t^{\prime }\right)-\gamma \left(t\right)$, and $h\left(\gamma \right)$ is the damping function (discussed below in detail). The function, $M(t-t^{\prime })$ is a memory kernel that is defined by the rate of change in the relaxation modulus (with physical dimensions of Pa/s):

$$M(t-t^{\prime })=\frac{\partial G(t-t^{\prime })}{\partial t^{\prime }}$$

(26)

Equation (25) is, evidently, a non-linear generalization of Equation (10) and we can use this for constructing the non-linear viscoelastic constitutive equations of fractional order by a suitable choice of the memory function. Thus, for a Fractional Maxwell Model with a relaxation modulus given by Equation (15), the memory function is written as:

$$M(t-t^{\prime })=-\mathbb{V}{\tau }_c^{-(1+\alpha )}{\left((t-t^{\prime })/{\tau }_c\right)}^{-1-\beta }{E}_{\alpha -\beta ,-\beta }\left(-{\left((t-t^{\prime })/{\tau }_c\right)}^{\alpha -\beta }\right)$$

(27)

In polymeric systems, the damping function $h\left(\gamma \right)$ that enters Equation (25) is related to the survival probability of network chains after the application of a large deformation [30,36]. The non-linear response of a given material can be calculated from measurements of the strain-dependent relaxation modulus $G(t,\gamma )$. Values of the strain-dependent relaxation moduli for a 0.5 wt.% solution of xanthan gum in water are shown in Figure 3a. For strain amplitudes up to $\gamma =0.4$, each measurement of the relaxation modulus overlapped on top of one another, corresponding to the linear viscoelastic response. This response is well described by the non-linear regression of the data to the FMM model with $\mathbb{G}=11.6$ Pa·${\mathrm{s}}^{\beta }$, $\beta =0.2$, $\mathbb{V}=458$ Pa·s, and $\alpha =1$. The characteristic relaxation time (corresponding to an appropriate moment of the relaxation spectrum) is thus ${\tau }_c={\left(\mathbb{V}/\mathbb{G}\right)}^{1/(\alpha -\beta )}$ = 98.9 s. However, for strain amplitudes $\gamma >0.4$, there is a clear softening of the material with large deformations reducing the strain-dependent relaxation modulus. In fact, at the largest strain ($\gamma =6$), the relaxation modulus of this complex fluid at any time t decreases by a factor of fifty. It is also clear that the functional form of each measured response in Figure 3a is self-similar and the data show that the material follow a principle of time-strain separability [54], so that $G(t,\gamma )\cong G\left(t\right)h\left(\gamma \right)$ is consistent with Equation (25). The damping function for this complex fluid can, thus, be found directly from experimental measurements by collapsing the strain-dependent relaxation modulus onto the linear relaxation modulus as shown in Figure 3b using the relation:

$$h\left(\gamma \right)=\frac{G(t,\gamma )}{G\left(t\right)}$$

(28)

The corresponding values of $h\left(\gamma \right)$ obtained at each imposed strain are plotted in Figure 3c. We compared the measurements with the predictions of three different damping functions proposed in the literature to describe the strain softening: exponential damping $h\left(\gamma \right)=exp(-\gamma /{\gamma }^{*})$ [31], Doi–Edwards damping $h\left(\gamma \right)=1/\left(1+{\left(\gamma /{\gamma }^{*}\right)}^2\right)$ [42,55], and the discontinuous Tanner–Simmons damping function $h\left(\gamma \right)=\mathcal{H}(1-\gamma /{\gamma }^{*})$ [56], where $\mathcal{H}$ is the Heaviside step function. The value of the critical strain ${\gamma }^{*}$ used in drawing Figure 3c was evaluated from fitting the Doi–Edwards damping function to the measurements of $h\left(\gamma \right)$, and was found to be ${\gamma }^{*}=0.89$. For comparative purposes, rather than refitting ${\gamma }^{*}$ for each model, we use the same value of ${\gamma }^{*}=0.89$ also for plotting the other damping functions. We note that both the exponential damping function and the Doi–Edwards damping function describe the measured values of the damping function much better than the simple Tanner–Simmons damping function and, hence, they have been used extensively in recent years [30,34] to capture the strain softening observed in a wide range of complex fluids such as polymer melts, liquid foods, and biopolymer solutions.

For a steady shear flow with a constant shear rate $\dot{\gamma }$, Equation (25) can be rewritten using the change of variables $u=t-t^{\prime }$ as:

$$\sigma \left(\dot{\gamma }\right)={\int }_0^{\infty }M\left(u\right)h\left(\dot{\gamma }u\right)\dot{\gamma }udu$$

(29)

The steady shear viscosity can, thus, be written as:

$$\eta \left(\dot{\gamma }\right)\equiv \frac{\sigma \left(\dot{\gamma }\right)}{\dot{\gamma }}={\int }_0^{\infty }M\left(u\right)h\left(\dot{\gamma }u\right)udu$$

(30)

As the shear rate is progressively incremented, it is clear that the impact of the strong decrease in $h\left(\gamma \right)$ with accumulated strain is to reduce the area integral embodied by the right hand side of Equation (30) and reduce the total shear viscosity; thus, giving rise to the shear thinning observed experimentally. Substituting the linear fractional kernel given by Equation (27) in Equation (30) and replacing $t-t^{\prime }=u$, the resulting fractional Wagner model predicts that the shear viscosity is given by:

$$\eta \left(\dot{\gamma }\right)=-\mathbb{V}{\tau }_c^{\beta -\alpha }{\int }_0^{\infty }{u}^{-\beta }{E}_{\alpha -\beta ,-\beta }\left(-{\left(u/{\tau }_c\right)}^{\alpha -\beta }\right)h\left(\dot{\gamma }u\right)du$$

(31)

which can be evaluated when the appropriate form of the damping function $h\left(\dot{\gamma }u\right)$ is specified. Evaluating the above integral results in a rate-dependent steady shear viscosity that can not be predicted by linear fractional constitutive models such as the Fractional Maxwell Model (FMM). In the subsequent sections, we find analytic results for the rate-dependent viscosity for the case of an exponential damping function and, then, extend our analysis for a more generalized form of damping function.

3.1. Exponential Damping

For an exponential damping function, substituting $h\left(\gamma \right)=exp(-\gamma /{\gamma }^{*})$ and recalling that in a steady shear flow the accumulated strain is $\gamma (t,t^{\prime })=\dot{\gamma }\times (t^{\prime }-t)=-\dot{\gamma }u$, in Equation (31), the rate-dependent viscosity can be written in a dimensionless form as:

$$\frac{\eta \left(\dot{\gamma }\right)}{\mathbb{V}{\tau }_c^{1-\alpha }}={\tau }_c^{\beta -1}{\int }_0^{\infty }-u^{-\beta }{E}_{\alpha -\beta ,-\beta }\left(-{\left(u/{\tau }_c\right)}^{\alpha -\beta }\right)exp(-\dot{\gamma }u/{\gamma }^{*})du$$

(32)

Using the following recurrence property of the MLF,

$$E_{a,b}\left(z\right)=\frac{b}{\mathsf{\Gamma }(1+b)}+z{E}_{a,b+a}\left(z\right).$$

(33)

we simplify Equation (32) to find a new expression for the steady shear viscosity:

$$\begin{array}{cc}\hfill \eta \left(\dot{\gamma }\right)=& \mathbb{V}{\tau }_c^{\beta -\alpha }\{\frac{\beta }{\mathsf{\Gamma }(1-\beta )}{\int }_0^{\infty }{u}^{-\beta }exp(-\dot{\gamma }u/{\gamma }^{*})du+\hfill \\ & {\left(\frac{1}{{\tau }_c}\right)}^{\alpha -\beta }{\int }_0^{\infty }{u}^{\alpha -2\beta }{E}_{\alpha -\beta ,\alpha -2\beta }\left(-{\left(u/{\tau }_c\right)}^{\alpha -\beta }\right)exp(-\dot{\gamma }u/{\gamma }^{*})du\}\hfill \end{array}$$

(34)

The integral in the first term of the brackets can be identified as a complete Gamma function with $x=\dot{\gamma }u/{\gamma }^{*}=\gamma /{\gamma }^{*}$, and can be expressed as:

$$\begin{array}{cc}\hfill \eta \left(\dot{\gamma }\right)=& \mathbb{V}{\tau }_c^{\beta -\alpha }\{\frac{\beta }{\mathsf{\Gamma }(1-\beta )}{\int }_0^{\infty }{x}^{-\beta }exp(-x)dx+\hfill \\ & {\left(\frac{1}{{\tau }_c}\right)}^{\alpha -\beta }{\int }_0^{\infty }{u}^{\alpha -2\beta }{E}_{\alpha -\beta ,\alpha -2\beta }\left(-{\left(u/{\tau }_c\right)}^{\alpha -\beta }\right)exp(-\dot{\gamma }u/{\gamma }^{*})du\}\hfill \end{array}$$

(35)

Thus, Equation (34) can be further simplified to:

$$\eta \left(\dot{\gamma }\right)=\mathbb{V}{\tau }_c^{\beta -\alpha }\left\{\beta {\left(\frac{\dot{\gamma }}{{\gamma }^{*}}\right)}^{\beta -1}+{\left(\frac{1}{{\tau }_c}\right)}^{\alpha -\beta }{\int }_0^{\infty }{u}^{\alpha -2\beta }{E}_{\alpha -\beta ,\alpha -2\beta }\left(-{\left(u/{\tau }_c\right)}^{\alpha -\beta }\right)exp(-\dot{\gamma }u/{\gamma }^{*})du\right\}$$

(36)

Using this recursive relationship, one can repeat this process and replace $E_{\alpha -\beta ,\alpha -2\beta }$ with its corresponding recursive representation in Equation (33). Hence, we write Equation (36) in the form of an infinite series:

$$\eta \left(\dot{\gamma }\right)=\mathbb{V}{\tau }_c^{\beta -\alpha }{\left(\dot{\gamma }/{\gamma }^{*}\right)}^{\beta -1}\sum _{k=1}^{\infty }{(-1)}^{k+1}\left(\alpha -k(\alpha -\beta )\right){\left(\dot{\gamma }\tau /{\gamma }^{*}\right)}^{-(k-1)(\alpha -\beta )}$$

(37)

The above equation is valid if the kth term of the infinite series approaches 0 as $k\to \infty$ [59]. Thus, the infinite series in Equation (37) is convergent and equivalent to Equation (32) only when the argument $\left(\dot{\gamma }{\tau }_c/{\gamma }^{*}\right)>1$. The dimensionless product $\dot{\gamma }{\tau }_c$ is known as the Weissenberg number, Wi (see Dealy [60] for additional details), which is a parameter that measures the strength of a given non-linear flow. It was clear that the relevant measure in our study emerges as a comparison between the strength of the flow and the value of the strain at which substantial damping in the relaxation modulus is encountered. Accordingly, the steady shear rate-dependent viscosity can be expressed in the following form:

$$\frac{\eta \left(\dot{\gamma }\right)}{\mathbb{V}{\tau }_c^{1-\alpha }}={\left(\frac{\mathrm{Wi}}{{\gamma }^{*}}\right)}^{\beta -1}\sum _{k=1}^{\infty }{(-1)}^{k+1}\left(\alpha -k(\alpha -\beta )\right){\left(\frac{\mathrm{Wi}}{{\gamma }^{*}}\right)}^{-(k-1)(\alpha -\beta )}$$

(38)

In order to find an analogous analytical expression between the flow viscosity and the applied shear rate for weak flows in which $\mathrm{Wi}/{\gamma }^{*}<1$, we used the following recurrence property for the MLF [51]:

$$E_{a,b}\left(z\right)=\frac{1}{z}\left[E_{a,b-a}\left(z\right)-\frac{1}{\mathsf{\Gamma }(b-a)}\right]$$

(39)

Using Equation (39) and a similar recursive procedure to that outlined above, we obtain an expression for the steady shear viscosity, when $\mathrm{Wi}/{\gamma }^{*}<1$:

$$\frac{\eta \left(\dot{\gamma }\right)}{\mathbb{V}{\tau }_c^{1-\alpha }}={\left(\frac{\mathrm{Wi}}{{\gamma }^{*}}\right)}^{\alpha -1}\sum _{k=1}^{\infty }{(-1)}^{k+1}\left(\beta +k(\alpha -\beta )\right){\left(\frac{\mathrm{Wi}}{{\gamma }^{*}}\right)}^{(k-1)(\alpha -\beta )}$$

(40)

This infinite series is, again, a valid representation of Equation (32) if the kth term of the infinite series Equation (40) approaches 0 as $k\to \infty$[59]. This requirement is satisfied only for $\mathrm{Wi}/{\gamma }^{*}<1$.

We can combine Equation (40) and Equation (38) using a Heaviside function $\mathcal{H}$ to write an expression for the steady shear viscosity that is valid for all Wi as:

$$\begin{array}{cc}\hfill \frac{\eta \left(\dot{\gamma }\right)}{\mathbb{V}{\tau }_c^{1-\alpha }}=& \mathcal{H}\left(1-\mathrm{Wi}/{\gamma }^{*}\right){\left(\frac{\mathrm{Wi}}{{\gamma }^{*}}\right)}^{\alpha -1}\sum _{k=1}^{\infty }{(-1)}^{k+1}\left(\beta +k(\alpha -\beta )\right){\left(\frac{\mathrm{Wi}}{{\gamma }^{*}}\right)}^{(k-1)(\alpha -\beta )}\hfill \\ & +\mathcal{H}\left(\mathrm{Wi}/{\gamma }^{*}-1\right){\left(\frac{\mathrm{Wi}}{{\gamma }^{*}}\right)}^{\beta -1}\sum _{k=1}^{\infty }{(-1)}^{k+1}\left(\alpha -k(\alpha -\beta )\right){\left(\frac{\mathrm{Wi}}{{\gamma }^{*}}\right)}^{-(k-1)(\alpha -\beta )}\hfill \end{array}$$

(41)

This series solution for the rate-dependent viscosity of a fractional Maxwell memory function with an exponential damping function is illustrated in Figure 4a for the FML and Figure 4b for the FMM, respectively, for various numbers of summation terms. The parameters used in Figure 4a,b were $\beta =0.3$, $\alpha =1$, ${\tau }_c=1$ s, ${\gamma }^{*}=1$ and $\beta =0.3$, $\alpha =0.8$, ${\tau }_c=1$ s, ${\gamma }^{*}=1$, respectively. The steady shear viscosity decreases by over a factor of 1000, with increasing Weissenberg number. At high shear rates, the viscosity approaches a power-law shear thinning response with $\eta \left(\dot{\gamma }\right)\sim {\dot{\gamma }}^{\beta -1}$. The series solution for the rate-dependent viscosity approaches the accurate numerical calculation of the integral in Equation (32) as the number of summed terms (K) increased. We observe that $K=O\left(10\right)$ gives an extremely good approximation of the rate-dependent viscosity. A similar series expansion can be implemented for the Tanner–Simmons damping function as well (see Appendix A for more details).

The exponential damping expression captures the experimentally observed levels of damping quite well (especially if the value of ${\gamma }^{*}$ was to be further adjusted to minimize the mean square error between the function $exp(-\gamma /{\gamma }^{*})$ and the data in Figure 3c), and it is widely used in the literature since it provides the possibility of obtaining solutions with different relaxation kernels in the integral expression given by Equation (31). Vermant et al. in [31] showed that the damping function of a 4 wt. % poly-isobutylene in a decalin solution follows exponential damping at large strains and, also, measured the corresponding steady shear viscosity. However, a more general damping function is needed to correctly capture the rate-dependent response of most materials. In particular, we note that for the exponential damping function $\partial h\left(\gamma \right)/\partial \gamma \ne 0$ as $\gamma \to 0$ so the range of linear viscoelastic behavior is in fact vanishingly small, whereas most real materials exhibited a clearly defined (and measurable) linear viscoelastic behavior for applied strains less than some characteristic value ${\gamma }^{*}$.

3.2. Generalized Damping Function

We introduced a generalized damping function (also known as the Soskey–Winter damping) which can capture the functional form of the strain softening observed in most real materials at large deformations [43]:

$$h\left(\gamma \right)=\frac{1}{1+{\left(\gamma /{\gamma }^{*}\right)}^m}$$

(42)

where the damping exponent $m>0$ and ${\gamma }^{*}$ is the critical strain of the material. The form of damping function in Equation (42) has been shown to accurately describe the damping effects or strain softening due to large strains for polymer melts, polymers, suspensions, magnetorheological fluids, polymer blends, and also in food rheology [30,43,57,61,62,63,64,65,66,67,68]. In [21,66], both the damping of the form given in Equation (42) and the resulting steady shear viscosity for xanthan gum solutions as well as for a lubricating grease are illustrated, respectively. The well-known Doi–Edwards constitutive model for entangled polymers predicts strain softening with a damping function that is a special case of Equation (42) with $m=2$ [42]. Using the recursive relationship procedure outlined in Section 3.1, the rate-dependent viscosity for a damping function $h\left(\gamma \right)$ of the form in Equation (42) can be derived to be:

$$\begin{array}{cc}\hfill \frac{\eta \left(\dot{\gamma }\right)}{\mathbb{V}{\tau }_c^{1-\alpha }}=& \mathcal{H}\left(1-\mathrm{Wi}/{\gamma }^{*}\right){\left(\frac{\mathrm{Wi}}{{\gamma }^{*}}\right)}^{\alpha -1}\sum _{k=1}^{\infty }{(-1)}^{k+1}p{\left(\frac{\mathrm{Wi}}{{\gamma }^{*}}\right)}^{(k-1)(\alpha -\beta )}\frac{1}{\mathsf{\Gamma }(1-p)}{\int }_0^{\infty }\frac{x^{-p}}{1+x^m}dx\hfill \\ & +\mathcal{H}\left(\mathrm{Wi}/{\gamma }^{*}-1\right){\left(\frac{\mathrm{Wi}}{{\gamma }^{*}}\right)}^{\beta -1}\sum _{k=1}^{\infty }{(-1)}^{k+1}q{\left(\frac{\mathrm{Wi}}{{\gamma }^{*}}\right)}^{-(k-1)(\alpha -\beta )}\frac{1}{\mathsf{\Gamma }(1-q)}{\int }_0^{\infty }\frac{x^{-q}}{1+x^m}dx\hfill \end{array}$$

(43)

where $p=\beta +k(\alpha -\beta )$, $q=\alpha -k(\alpha -\beta )$ and $x=\gamma /{\gamma }^{*}$. Using the Euler reflection formula:

$${\int }_0^{\infty }\frac{x^a}{1+x^m}dx=\mathsf{\Gamma }\left(\frac{1+a}{m}\right)\mathsf{\Gamma }\left(1-\frac{1+a}{m}\right)=\frac{\pi }{msin\left[\pi \left(\frac{1+a}{m}\right)\right]}\mathrm{for}0<\frac{1+a}{m}<1$$

(44)

and identifying $a=-p$ and $a=-q$ in the first and second term of Equation (43), respectively, we can simplify the rate-dependent shear viscosity as follows:

$$\begin{array}{cc}\hfill \frac{\eta \left(\dot{\gamma }\right)}{\mathbb{V}{\tau }_c^{1-\alpha }}=& \mathcal{H}\left(1-\mathrm{Wi}/{\gamma }^{*}\right){\left(\frac{\mathrm{Wi}}{{\gamma }^{*}}\right)}^{\alpha -1}\sum _{k=1}^{\infty }\left[{(-1)}^{k+1}\frac{p\pi }{\mathsf{\Gamma }(1-p)sin\left(\pi \frac{(1-p)}{m}\right)}{\left(\frac{\mathrm{Wi}}{{\gamma }^{*}}\right)}^{(k-1)(\alpha -\beta )}\right]\hfill \\ & +\mathcal{H}\left(\mathrm{Wi}/{\gamma }^{*}-1\right){\left(\frac{\mathrm{Wi}}{{\gamma }^{*}}\right)}^{\beta -1}\sum _{k=1}^{\infty }\left[{(-1)}^{k+1}\frac{q\pi }{\mathsf{\Gamma }(1-q)sin\left(\pi \frac{(1-q)}{m}\right)}{\left(\frac{\mathrm{Wi}}{{\gamma }^{*}}\right)}^{-(k-1)(\alpha -\beta )}\right]\hfill \end{array}$$

(45)

Importantly, we note that the precise form of the damping function does not affect the powers of the polynomial terms in Equation (45) just the magnitude of the coefficients in the summation. It is noteworthy to mention that the first term in the brackets of the series expansion given above in Equation (45) is valid only in the case of:

$$0<(1-p)/m<1$$

(46)

for the first series sum, and the second series term is valid when:

$$0<(1-q)/m<1$$

(47)

The asymptotes at low and high Wi are dominated by the $(k=1)$ terms in Equation (45). The asymptotes can be directly calculated from Equation (45), and for cases when the inequalities in Equation (46) or Equation (47) are not satisfied (i.e., when $(1-\alpha )<m<(1-\beta ))$ or $m<1-\alpha <1-\beta$), we evaluated the asymptotes from Equation (31) (see Appendix A for more details).

For completeness, we provide the asymptotic behavior of the rate-dependent shear viscosity for both the FMM and FML model with non-linear damping and different values of the damping function parameter m in

Table 2. Asymptotic forms of the steady shear viscosity predicted by the FMM ($0<\beta <\alpha <1$) for various values of the damping exponent (m) relative to the FML parameters $\alpha$ and $\beta$, that characterize the linear viscoelastic response of complex multiscale materials. The function $f\left(m\right)$ in the table is given by: $f\left(m\right)={\int }_0^{\infty }-x^{-\beta -m}{E}_{\alpha -\beta ,-\beta }\left(-x^{\alpha -\beta }\right)dx$.

$\mathit{\eta }/\left({\mathit{G}}_{\mathit{c}}{\mathit{\tau }}_{\mathit{c}}\right)$	Wi$/{\mathit{\gamma }}^{*}\ll 1$	Wi$/{\mathit{\gamma }}^{*}\gg 1$
$1-\alpha <1-\beta <m$ (strong damping)	$\frac{{\displaystyle \pi \alpha }}{{\displaystyle \mathsf{\Gamma }(1-\alpha )sin\left[\frac{\pi (1-\alpha )}{m}\right]}}{\left(\frac{{\displaystyle \mathrm{Wi}}}{{\displaystyle {\gamma }^{*}}}\right)}^{\alpha -1}$	$\frac{{\displaystyle \pi \beta }}{{\displaystyle \mathsf{\Gamma }(1-\beta )sin\left[\frac{\pi (1-\beta )}{m}\right]}}{\left(\frac{{\displaystyle \mathrm{Wi}}}{{\displaystyle {\gamma }^{*}}}\right)}^{\beta -1}$
$1-\alpha <m\le 1-\beta$ (weak damping)	$\frac{{\displaystyle \pi \alpha }}{{\displaystyle \mathsf{\Gamma }(1-\alpha )sin\left[\frac{\pi (1-\alpha )}{m}\right]}}{\left(\frac{{\displaystyle \mathrm{Wi}}}{{\displaystyle {\gamma }^{*}}}\right)}^{\alpha -1}$	$f\left(m\right){\left(\frac{{\displaystyle \mathrm{Wi}}}{{\displaystyle {\gamma }^{*}}}\right)}^{-m}$
$1-\alpha <m\le 1-\beta$ (v. weak damping)	Unbounded stress growth	Unbounded stress growth

and

Table 3. Asymptotic forms of the dimensionless shear viscosity $\eta \left(\dot{\gamma }\right)/\mathbb{V}$ for the FML model ($\alpha =1$) for various values of the damping exponent (m) relative to the parameter $\beta$ that characterizes the high frequency linear viscoelastic moduli. The function $f\left(m\right)$ in the table is given by: $f\left(m\right)={\int }_0^{\infty }-x^{-\beta -m}{E}_{\alpha -\beta ,-\beta }\left(-x^{\alpha -\beta }\right)dx$.

$\mathit{\eta }/\left({\mathit{G}}_{\mathit{c}}{\mathit{\tau }}_{\mathit{c}}\right)$	Wi$/{\mathit{\gamma }}^{*}\ll 1$	Wi$/{\mathit{\gamma }}^{*}\gg 1$
$1-\beta <m$ (strong damping)	1	$\frac{{\displaystyle \pi \beta }}{{\displaystyle \mathsf{\Gamma }(1-\beta )sin\left[\frac{\pi (1-\beta )}{m}\right]}}{\left(\frac{{\displaystyle \mathrm{Wi}}}{{\displaystyle {\gamma }^{*}}}\right)}^{\beta -1}$
$m<1-\beta$ (weak damping)	1	$f\left(m\right){\left(\frac{{\displaystyle \mathrm{Wi}}}{{\displaystyle {\gamma }^{*}}}\right)}^{-m}$

, respectively.

From the asymptotes, it is clear that at low shear rates, the exponent of the power-law behavior of the steady shear viscosity is set by the FMM parameter $\alpha$, i.e., the power-law exponent that characterizes the linear viscoelastic properties of the fluid at low $\mathrm{De}=\omega {\tau }_c\ll 1$ (cf. Figure 2). At high shear rates $\mathrm{Wi}/{\gamma }^{*}\gg 1$, the power-law behavior of the viscosity depends on the strength of damping. For a sufficiently strong damping function (i.e., $m>1-\beta$), the FMM parameter $\beta$ characterizing the linear viscoelastic response at $\mathrm{De}=\omega {\tau }_c\gg 1$ sets the rate of shear thinning. However, for a weaker choice of the damping function (i.e., $m<1-\beta$), the rate-thinning in the shear viscosity is in fact determined by the damping parameter m. For the case of very weak (insufficient) damping (i.e., $m<1-\alpha$ ), we can not obtain a meaningful value for the steady shear viscosity (i.e., unbounded stress growth in time is expected during start up of a steady shear flow, and no steady shear viscosity is reached).

Similar arguments can be extended to the FML case (corresponding to $\alpha =1$), and the results are tabulated in Table 3. For $\alpha =1$, one of the spring-pots becomes a purely viscous dashpot (i.e., $\mathbb{V}\to {\eta }_0)$, and we observe a plateau in the shear viscosity equal to ${\eta }_0$ at low shear rates (irrespective of the form of the damping function). At high shear rates, the shear thinning index is set by either the FML parameter $\beta$ (for the case of strong damping with $(m>1-\beta$) or by the exponent m for weak damping with $(m<1-\beta )$.

This asymptotic behavior of the steady shear viscosity and its dependence on the strength of the non-linear damping is illustrated visually for the FMM and FML models in Figure 5. The rate of shear thinning at high shear rates is independent of m for all cases when $m>1-\beta$. It is also very close to the result obtained for exponential damping (shown by the black dashed line). This suggests that as long as the non-linear strain damping described by Equation (42) is sufficently strong ( $m>1-\beta$), the rate of shear thinning in the steady shear viscosity is set by the linear viscoelastic properties ($\alpha ,\beta$) of the material. This analysis provides a rationalization for why a variety of different functional forms for the damping functions used in the literature are all able to provide a good prediction of the rate-dependent viscosity.

Our approach of incorporating non-linearities into a fractional description of viscoelastic behavior using an FMM linear kernel in the Wagner constitutive equation can be further extended to an even wider class of materials such as non-soft materials or viscoelastic solids by replacing the FMM kernel with a Fractional Zener kernel, Fractional Kelvin–Voigt [26], or the more generalized Modified Fractional Maxwell (MFM) introduced by Xu et al. [69,70]. The MFM has been shown to accurately describe the linear viscoelastic properties of complex liquids, rubber-like materials, and glassy viscoelastic solids across a wide range of frequencies.

4. Approximate Form of the Steady Shear Viscosity

The zero shear viscosity of the FMM diverges when $\alpha \ne 1$[21] as seen from Figure 5b. However, most complex fluids, such as biopolymer solutions or polymer melts, exhibit a distinct, measurable finite value of the shear viscosity in the limit of small shear-rates and are, thus, best described by fractional models in the FML limit ($\alpha =1$). The full viscosity for the non-linear formulation of the FML curve can be determined from the integral expression (Equation (32)) or from the series expansions given in Equations (41) and (45). However, for regression to data and the fitting of material constants, it is helpful to have a simpler closed-form expression for the flow curve. Thus, we propose a compact analytical approximation for the rate-dependent viscosity based on the results we obtained in Section 3 with an FML kernel. We express the rate-dependent viscosity in terms of a four parameter function:

$$\frac{\eta }{{\eta }_0}=\frac{1}{1+(1/B){\left(\mathrm{Wi}/{\gamma }^{*}\right)}^s}$$

(48)

where $\mathrm{Wi}={\tau }_c\dot{\gamma }$ and ${\tau }_c={\left(\mathbb{V}/\mathbb{G}\right)}^{1/(1-\beta )}$, and the constants B and s are power-law parameters that remain to be determined. This general functional form is motivated by commonly used empirical expressions such as the Carreau model or Cross model [1,71,72,73] used in the rheology literature for describing experimentally measured flow curves of $\eta \left(\dot{\gamma }\right)$. For strong damping with $m>1-\beta$, the two adjustable parameters in Equation (48) are given by:

$$\begin{array}{cc}\hfill & s=1-\beta \hfill \end{array}$$

(49a)

$$\begin{array}{cc}\hfill & B=\frac{\pi \beta }{\mathsf{\Gamma }(1-\beta )sin\left[\pi \left(\frac{1-\beta }{m}\right)\right]}\hfill \end{array}$$

(49b)

Whereas, for weak damping with $m<1-\beta$, the parameters were:

$$\begin{array}{cc}\hfill & s=m\hfill \end{array}$$

(50a)

$$\begin{array}{cc}\hfill & B=f\left(m\right)={\int }_0^{\infty }-x^{-\beta -m}{E}_{1-\beta ,-\beta }\left(-x^{1-\beta }\right)dx\hfill \end{array}$$

(50b)

To evaluate the rate-dependent viscosity from Equation (48), we require the values of the three parameters in the FML model that describe the linear rheological properties of the complex fluid ($\mathbb{G},\beta ,\mathbb{V}\equiv {\eta }_0$) and the damping exponent m for describing the strain dependence of the relaxation modulus $G(t,\gamma )$. In Figure 6a, we compare the compact expression for the rate-dependent viscosity given in Equation (48) with a numerical evaluation of Equation (31) for two different values of m: a strong damping case, $m=2$; and a weak damping case, $m=0.2(<(1-\beta \left)\right)$. It can be seen that Equations (48)–(50) capture the correct asymptotic limits at high and low shear rates, but do not precisely capture the gradual transition in the shear viscosity around $\mathrm{Wi}/{\gamma }^{*}\sim O\left(1\right)$. A contour plot comparing the proposed viscosity approximation given in Equation (48) with the numerical integration of Equation (31) using the generalized damping function of Equation (42) is presented in Figure 6b. The error is only substantial for the unlikely condition where $m=1-\beta$, i.e., when the damping parameter and the FML power-law exponents ($1-\beta$) are equal to each other. By contrast, for strong damping cases with $m\ge 2$, the error incurred in using the viscosity form given by Equations (48) and (49) is always less than 7.5% across the whole range of shear rates.

5. Conclusions

In this contribution, we have sought to extend fractional calculus descriptions of soft matter to the non-linear domain that is important for modeling the real flows of soft shear-sensitive complex fluids. To do this, we made use of the close connection between the Caputo definition of the fractional derivative and the general Rivlin–Sawyers hereditary integral form of a viscoelastic material. The steps of our analysis may be summarized as follows:

• We considered a general model known as the Fractional Maxwell Model (FMM) as well as two important limits studied in the fractional viscoelasticity literature (the Fractional Maxwell Liquid and the Fractional Maxwell Gel). The limiting frequency responses of the linear viscoelastic moduli for these models are summarized in Table 1.
• By considering measurements of the strain-dependent relaxation moduli for a 0.5 wt.% aqueous solution of xanthan gum, we illustrated the time-strain separability in $G(t,\gamma )$ observed in many complex fluids. This motivated the use of a damping function $h\left(\gamma \right)$ and the time-strain separable Wagner model in conjunction with a fractional relaxation modulus described by the Mittag-Leffler Function (MLF) to capture the strain softening that arises at large deformations in complex fluids.
• Using the FMM as the linear kernel together with the Wagner integral framework, we derived analytical expressions for the rate-dependent viscosity in a steady simple shear flow for an exponential damping function given by Equation (41), as well as a more general damping function given by Equation (45). From this rate-dependent viscosity expression, we concluded that for weak shear flows (i.e., low shear rates with $\mathrm{Wi}\ll 1$), the rate-dependence of the viscosity is always set by the FMM parameter $\alpha$ and approaches a constant plateau value $\mathbb{V}\to {\eta }_0$ for the FML model when $\alpha =1$.
• For strong shear flows (i.e., high shear rates with $\mathrm{Wi}\gg 1$), the power-law exponent of the rate-dependent viscosity is set by the other fractional Maxwell exponent $\beta$ and scales as ${\dot{\gamma }}^{\beta -1}$ for strong damping ($m>1-\beta$), whereas it depends on the damping exponent m and scales as ${\dot{\gamma }}^{-m}$ for weak damping ($m<1-\beta$). For sufficiently weak strain damping ($m<1-\alpha$), when combined with the FMM kernel, the shear stress following the inception of steady shear flow grows unbounded with time and a bounded shear stress is never attained.
• Finally, we developed an approximate, closed-form expression given by Equations (48)–(50) for the steady shear viscosity that incorporates the three parameters of the linear viscoelastic Fractional Maxwell Liquid model (${\eta }_0=\mathbb{V},\mathbb{G},\beta$) with a non-linear damping function specified by two non-linear parameters (m and ${\gamma }^{*}$).

Our analysis, thus, provides a compact analytical formulation that enabled the versatility provided by fractional calculus descriptions of soft matter and fractal media to be more broadly adopted by rheologists in the industry and beyond.