Recent Advances in Complex Fluids Modeling

Abstract

In this chapter, we present a brief description of existing viscoelastic models, starting with the classical differential and integral models, and then focusing our attention on new models that take advantage of the enhanced properties of the Mittag-Leffler function (a generalization of the exponential function). The generalized models considered in this work are the fractional Kaye-Bernstein, Kearsley, Zapas (K-BKZ) integral model and the differential generalized exponential Phan-Thien and Tanner (PTT) model recently proposed by our research group. The integral model makes use of the relaxation function obtained from a step-strain applied to the fractional Maxwell model, and the differential model generalizes the familiar exponential Phan-Thien and Tanner constitutive equation by substituting the exponential function of the trace of the stress tensor by the Mittag-Leffler function. Since the differential model is based on local operators, it reduces the computational time needed to predict the flow behavior, and, it also allows a simpler description of complex fluids. Therefore, we explore the rheometric properties of this model and its ability (or limitations) in describing complex flows.

1. Introduction

Since viscoelastic materials are abundant in nature and present in our daily lives (examples are paints, blood, polymers, biomaterials, etc.), it is important to study and understand viscoelastic behavior. Therefore, in this chapter, we further develop the modeling of viscoelasticity making use of fractional calculus tools.

We start this section with some basic concepts that are needed to derive and understand classical and fractional viscoelastic models. These are trivial concepts such as force, stress, viscosity, Hooke’s law of elasticity and also Newton’s law of viscosity. Later, we evolve to more complex concepts of viscoelasticity that involve the knowledge of fractional calculus, integral and differential models.

It is well known that a force is any interaction that when unopposed will change the motion of an object/body. Stress is an internal resistance provided by the body itself whenever it is under deformation. Stress is defined as the intensity of internal forces developed in the material. The intensity of any quantity is defined as the ratio of the quantity to the area on which it is acting, leading to: Average Stress = Force/Area. If we want to know the stress in one material point, then we must take the limit of the area to zero. A good example on how stress works is given by imagining a person lying on top of thin layer of ice. When the person is lying down on the ice, the force (weight) divided by the area of the surface of the person in contact with the ice is smaller, when compared to the case when someone is standing up (the weight is the same, but the area in contact with the ice is smaller). Therefore, eventually, the ice will break due to the high internal stresses when the person is standing. Finally, we refer to elasticity as the ability of a body to resist a distorting influence and to return to its original size and shape when that influence or force is removed. See for example Figure 1 where three springs are stretched. If we remove the weights attached to the springs, the spring would ideally return to its initial/natural position.

Figure 1.

Schematic of an experiment to obtain the relationship between force and deformation: (a) Experimental setup where three springs are stretched with the use of weights; (b) Graph showing the experimental results obtained from stretching three springs (the force is proportional to the deformation).

Figure 1(b) also shows an experiment where we observe that the force (mass times gravity) applied to the spring (increasing weight) is proportional to the displacement. This is known as Hooke’s law (the force $F$ needed to extend or compress a spring by some distance ${\mathrm{\gamma }}_e=\left(x-x_0\right)$ is proportional to that distance, $F=k{\gamma }_e$ ). Note that if we continue to increase the weight, eventually the spring will break. Therefore, Hooke’s law is a very good linear approximation of what happens in the real world.

We will now explore the concept of viscosity in fluids. The viscosity of a fluid is a measure of the internal resistance to the rate of deformation:

As an example, imagine that we have a thin film of fluid in between two parallel plates, as shown in Figure 2. The fluid is at rest, and suddenly the upper plate starts moving with constant velocity U. This velocity will be felt at the bottom layer due to diffusion of momentum, and to keep the bottom wall fixed, we must exert a restraining force, that is measured with a force gage or dynamometer attached to that wall. Note that if we take the view of this portion of fluid as infinitesimally thin layers, we observe that each layer will drag the underlying layer due to the action of viscosity (internal resistance). The higher the viscosity, the more force will be required to deform the fluid at a given speed U.

Figure 2.

Schematic of an experiment to verify Newton’s law of viscosity: (a) Liquid at rest between parallel plates; (b) The top wall is pulled with velocity U and a force meter is used to measure the force exerted on the bottom wall; (c) Experimental results.

Since the velocity of the thin layer adjacent to the top wall is U and the velocity of the bottom layer is 0, the velocity of each layer (for a Newtonian fluid) is given by u(y) = Uy/h, with y the coordinate shown in Figure 2(a). Figure 2(b) shows the experimental forces measured for different ratios of U/h. We observe that the force is proportional to U/h and U/h = du(y)/dy; therefore, we conclude the following (Newton’s law of viscosity):

$$\frac{\mathit{\text{Force}}}{\mathit{\text{Area}}}=\sigma =\eta \frac{U}{h}=\eta \frac{\mathit{du}\left(y\right)}{\mathit{dy}}\Rightarrow \sigma =\eta \frac{\mathit{du}\left(y\right)}{\mathit{dy}}$$ E1

with $\sigma$ the unidirectional stress and $\eta$ as the constant of proportionality, known as the Newtonian shear viscosity. Note that du/dy is known as the rate of shear deformation, usually denoted by $\dot{\mathrm{\gamma }}$ .

A good example of something we may see every day and something that verifies Newton’s law of viscosity is a dashpot. It is used for example as a door closer to prevent it from slamming shut.

1.1. Viscoelastic models

The simplest model that considers both viscous and elastic behavior is the linear Maxwell model [1] and can be obtained from a combination in series of a dashpot, $\sigma =\mathit{\eta d}{\mathrm{\gamma }}_f\left(t\right)/\mathit{dt}$ , and a spring, $\sigma =G{\gamma }_e\left(t\right)$ (with the subscripts f and e standing for Newtonian fluid and Hookean elastic solid, respectively), as shown in Figure 3.

Figure 3.

Maxwell model.

The total deformation $\gamma$ is the sum of the deformation obtained from the spring ${\gamma }_e$ and the dashpot ${\gamma }_f$ , and the rate of deformation is given by:

$$\begin{array}{l}\frac{\mathit{d\gamma }\left(t\right)}{\mathit{dt}}=\frac{d{\gamma }_f\left(t\right)}{\mathit{dt}}+\frac{d{\gamma }_e\left(t\right)}{\mathit{dt}}\\ \iff \frac{\mathit{d\gamma }\left(t\right)}{\mathit{dt}}=\frac{\sigma }{\eta }+\frac{1}{G}\frac{\mathit{d\sigma }}{\mathit{dt}}\\ \iff \underset{\mathit{\text{Maxwell\hspace{0.17em}\hspace{0.17em}Model}}}{\underbrace{\sigma +\lambda \frac{\mathit{d\sigma }}{\mathit{dt}}=\eta \frac{\mathit{d\gamma }\left(t\right)}{\mathit{dt}}}},\lambda =\frac{\eta }{G}\end{array}$$ E2

The three-dimensional version of this model can be easily obtained by considering appropriate tensors instead of the scalar properties of stress and deformation, leading to the following model:

$$\mathbf{\sigma }+\lambda \frac{d\mathbf{\sigma }}{\mathit{dt}}=\eta \frac{d\mathbf{\gamma }\left(\mathit{t}\right)}{\mathit{dt}}$$ E3

with $\mathbf{\sigma }$ the stress tensor, $\dot{\mathbf{\gamma }}=\left(\mathbf{\nabla }\mathit{u}+{\left(\mathit{\nabla }\mathit{u}\right)}^{\mathbf{T}}\right)$ the rate of deformation tensor, $\mathbf{u}$ the velocity vector, $\lambda$ the relaxation time of the fluid and $\eta$ the zero shear rate viscosity. This model can be equivalently written in integral form as

$$\mathbf{\sigma }\left(t\right)=\underset{0}{\overset{t}{\mathit{\int }}}{\mathit{Ge}}^{-\left(t-t^{\prime }\right)/\lambda }\frac{d\mathbf{\gamma }}{d{t}^{\prime }}d{t}^{\prime },$$ E4

where $G=\eta /\lambda$ and it was assumed that the fluid is at rest for $t<0$ .

The Maxwell model is not observer independent (frame invariant) and, therefore, the results obtained with this model may not be correct if large deformations are considered (e.g., we may obtain a viscosity that depends directly on the velocity rather than the velocity gradient, which is not correct, and is unphysical). To solve this problem, new models were proposed in the literature that can deal with this non-invariance problem.

Two well-known examples of frame invariant models are the upper-convected Maxwell (UCM) model given by $\mathbf{\sigma }+\lambda \widehat{\mathbf{\sigma }}=\eta \dot{\mathbf{\gamma }}$ (with $\widehat{\mathbf{\sigma }}=\partial \mathbf{\sigma }/\partial t+\mathbf{u}\cdot \nabla \mathbf{\sigma }-{\left(\nabla \mathbf{u}\right)}^T\cdot \mathbf{\sigma }-\mathbf{\sigma }\cdot \nabla \mathbf{u}$ the upper-convected derivative) that can also be written in integral form as

$$\mathbf{\sigma }\left(t\right)=\underset{0}{\overset{t}{\int }}\frac{\eta }{{\lambda }^2}{e}^{-\left(t-t^{\prime }\right)/\lambda }\left({\mathbf{C}}_{t^{\prime }}^{-1}-\mathbf{I}\right)d{t}^{\prime }$$ E5

where ${\mathbf{C}}^{-1}$ is the Finger strain tensor (a frame-invariant measure of deformation) [1]. The term $m\left(t-t^{\prime }\right)=\eta /{\lambda }^2{e}^{-\left(t-t^{\prime }\right)/\lambda }$ is known as the memory function (the derivative of the relaxation modulus $G\left(t-t^{\prime }\right)$ ). Note that the relaxation modulus can be easily obtained by imposing a step strain (constant deformation), as shown in Figure 4, resulting in $G\left(t\right)=\sigma /{\mathrm{\gamma }}_0=G{e}^{-\left(t-t^{\prime }\right)/\lambda }$ .

Figure 4.

Step strain of a Maxwell model. The step strain is given by $\gamma ={\gamma }_{0}H\left(t-t_0\right)$ with $H\left(t\right)$ the Heaviside function, and the stress relaxation is the solution of $\sigma +\lambda \mathit{d\sigma }/\mathit{dt}=\mathit{\eta d}{\gamma }_{0}H\left(t-t_0\right)/\mathit{dt}$ with $\sigma \left(t_0\right)={\sigma }_0$ , given by $\sigma ={\sigma }_0{e}^{-\left(t-t_0\right)/\lambda }$ ( ${\sigma }_0=G{\gamma }_0$ ).

Other well-known example of a frame-invariant but now nonlinear viscoelastic model is the variation of the K-BKZ [2] model proposed by Wagner, Raible and Meissner [3, 4],

$$\mathbf{\sigma }\left(t\right)=\underset{0}{\overset{t}{\int }}m\left(t-t^{\prime }\right)h\left(I_1,I_2\right)\left({\mathbf{C}}_{t^{\prime }}^{-1}-\mathbf{I}\right)d{t}^{\prime },$$ E6

where ${\mathbf{C}}^{-1}$ is the Finger tensor [1], $I_1,I_2$ are the traces of ${\mathbf{C}}^{-1}$ and $\mathbf{C}$ , respectively, and $h\left(I_1,I_2\right)$ is termed the damping function [5] (note that it is again assumed that the fluid is at rest for $t<0$ ). A large number of damping functions can be found in the literature (see [5]). The term $m\left(t-t^{\prime }\right)$ was proposed to be of the form:

$$m\left(t-t^{\prime }\right)=\frac{a}{\lambda }{e}^{-\left(t-t^{\prime }\right)/\lambda },$$ E7

where $a$ and $\lambda$ are model parameters. Note that the relaxation modulus is the response of the stress to a step in deformation (see Figure 4). It should be remarked that when $a=\eta /\lambda$ and $h\left(I_1,I_2\right)=1$ , we recover the integral version of the UCM model.

Different differential models were proposed in the literature along the years, with the aim of improving the modeling of complex viscoelastic materials, and with the aim of achieving the same modeling quality of integral models (by only using differential operators). Note that integral models are non-local (in time) operators that take into account all the past deformation of the fluid while differential models ones describe the material response in terms of the rate of change of stress to the local deformation, thus influencing the fitting quality of the model and the computational effort to numerically solve them (when performing numerical simulations).

More recently, new models have been proposed in the literature that basically take advantage of the generalization of the exponential function appearing in Eqs. (4), (5), and (7), thus allowing a more broad and accurate description of the relaxation of complex fluids (while the commonly used continuum approach describes the fluid as a whole, with only one relaxation, unless a Prony series is considered, that is, considering a series of the form $\sum _i{a}_i{e}^{-\left(t-t^{\prime }\right)/{\lambda }_i}$ ). This generalized function is the Mittag-Leffler function that naturally arises when solving problems involving fractional derivatives (more precisely, derivatives of non-integer order). This function will be introduced later in Section 3.

2. Fractional derivatives

To understand the need and the concept of a fractional derivative and its importance in the context of modeling physical processes, let us start with a simple example (Figure 5).

Figure 5.

Material formed by two regions where the same physical process occurs at different rates.

Imagine a portion of material that is principally formed of two different regions. In these regions, two similar physical processes ${\phi }_1$ and ${\phi }_2$ occur (for the time being it does not matter what is the process under study), but, at different rates, $d{\phi }_1/\mathit{dt}=0.1$ for Region I and $d{\phi }_2/\mathit{dt}=1$ for Region II. If we look at the portion of material as a whole, one would naturally choose the rate of $1$ as representative of the material’s behavior, because this region is bigger (when compared to Region I). However, this clearly neglects entirely the local variation in the deformation associated with the neighboring Region I. With the help of fractional calculus, we may define derivatives/rates of non-integer order, and we may have (for example) a rate given by $d^{\beta }\phi /{\mathit{dt}}^{\beta }$ with $\beta =0.9$ (possibly better representing the material behavior as a whole, by providing intermediate rates).

Although we have not defined yet what a fractional derivative is, the fact of having the possibility of non-integer derivatives seems quite attractive, allowing the creation of a continuous path between integer-order derivatives that may lead to a better description of the different rates of a certain physical process occurring in the same material. This means that fractional derivatives can transport more and more precise local information from the microscopic world to the continuum description.

2.1. Riemann-Liouville and Caputo fractional derivatives

Now, to understand a fractional derivative, we start by acknowledging that the n-fold integral of a generic function $f\left(t\right)$ is given by the formula

$$\underset{n\mathit{\text{times}}}{\underbrace{{\int }_a^t{\int }_a^t\cdots {\int }_a^{t}f\left(t\right)\mathit{\text{dtdt}}\dots \mathit{dt}=}}{J}_a^{n}f\left(t\right)=\frac{1}{\left(n-1\right)!}\underset{a}{\overset{t}{\int }}{\left(t-t^{\prime }\right)}^{n-1}f\left(t^{\prime }\right)d{t}^{\prime }.$$ E8

A generalization to non-integer values of $n$ can be performed using the Euler Gamma function $\mathrm{\Gamma }\left(x\right)$ , leading to the Riemann-Liouville fractional integral

$$J_a^{n}f\left(t\right)=\frac{1}{\mathrm{\Gamma }\left(\alpha \right)}{\int }_a^t{\left(t-t^{\prime }\right)}^{\alpha -1}f\left(t^{\prime }\right)d{t}^{\prime },$$ E9

where we have used $\alpha$ to represent the generalization of $n$ to non-integer values. A fractional derivative of any order can then be obtained by manipulating the number of integrations and differentiations of the function $f\left(t\right)$ . By performing the $m-\alpha$ -fold integration of the $m^{\mathit{th}}$ derivative of $f\left(t\right)$ , $J_a^{m-\alpha }{D}^{m}f\left(t\right)$ with $m=⌈\alpha ⌉$ , we arrive at the generalized derivative formula (Caputo fractional derivative [6]) of order $m-1<\alpha <m$ ,

$$\frac{d^{\alpha }f\left(t\right)}{{\mathit{dt}}^{\alpha }}=\frac{1}{\mathrm{\Gamma }\left(m-\alpha \right)}\underset{a}{\overset{t}{\int }}{\left(t-t^{\prime }\right)}^{-\alpha +m-1}\frac{d^{m}f\left(t^{\prime }\right)}{d{t^{\prime }}^m}d{t}^{\prime },m-1<\alpha <m,$$ E10

This last fractional derivative is the one chosen to deal with physical processes due to the ease in handling initial and boundary conditions [7].

Next, we present two models that rely on the Mittag-Leffler function (a function closely related to fractional calculus) to improve their modeling and fitting capabilities when describing the behavior of viscoelastic materials. These are the fractional K-BKZ (integral) and the generalized Phan-Thien and Tanner (differential) models.

3. Viscoelastic models based on the Mittag-Leffler function

3.1. The fractional K-BKZ model

We first note that the Maxwell-Debye relaxation of stress (exponential decay—see Eqs. (4) and (5)) is quite common, but there are many real materials showing different types of fading memory, such as a power law decay $G\left(t\right)\sim t^{-\alpha },0<\alpha <1$ [8]. For example, the critical gel model investigated by Winter and Chambon is written $G\left(t\right)={\mathit{St}}^{-\alpha }$ . If we assume the relaxation modulus for an arbitrary loading history in such materials is given by $G\left(t-t^{\prime }\right)=\mathbb{V}{\left(\mathrm{\Gamma }\left(1-\alpha \right)\right)}^{-1}{\left(t-t^{\prime }\right)}^{-\alpha }$ ( $\mathbb{V}$ is known as a quasi-property [9] and is connected to the critical gel strength by $S=\mathbb{V}/\mathrm{\Gamma }\left(1-\alpha \right)$ ), then we have that:

$$\mathbf{\sigma }\left(t\right)=\frac{1}{\mathrm{\Gamma }\left(1-\alpha \right)}\underset{0}{\overset{t}{\int }}\mathbb{V}{\left(t-t^{\prime }\right)}^{-\alpha }\frac{d\mathbf{\gamma }}{d{t}^{\prime }}d{t}^{\prime }.$$ E11

By recognizing that the Caputo fractional derivative of a general function $\mathbf{\gamma }\left(t\right)$ (in our case $\mathbf{\gamma }\left(t\right)$ is the deformation) is defined as [10]:

$$\frac{d^{\alpha }\mathbf{\gamma }\left(t\right)}{{\mathit{dt}}^{\alpha }}=\frac{1}{\mathrm{\Gamma }\left(1-\alpha \right)}\underset{0}{\overset{t}{\int }}{\left(t-t^{\prime }\right)}^{-\alpha }\frac{d\mathbf{\gamma }}{d{t}^{\prime }}d{t}^{\prime },$$ E12

we obtain a generalized viscoelastic model [10, 11], that can be written in the simple compact form:

$$\mathbf{\sigma }=\mathbb{V}\frac{d^{\alpha }\mathbf{\gamma }\left(t\right)}{{\mathit{dt}}^{\alpha }},0<\alpha <1,$$ E13

This model provides a generalized viscoelastic response, in the sense that when $\alpha =1$ we obtain a Newtonian fluid, and when $\alpha =0$ we obtain a Hookean elastic solid. The corresponding mechanical element is intermediate to the spring and dashpot shown in Figure 3 and is thus known as a spring-pot [11, 12]. Note that care must be taken when $\alpha \to 1$ because of the singularity in $\mathrm{\Gamma }\left(1-\alpha \right)$ [12].

We can define the fractional Maxwell model (FMM) as a combination of two linear fractional elements (spring-pots) in series. In a series configuration, the stress felt by each spring-pot is the same, that is, $\mathbf{\sigma }=\mathbb{V}{d}^{\alpha }{\mathbf{\gamma }}_1\left(t\right)/{\mathit{dt}}^{\alpha }=\mathbb{G}{d}^{\beta }{\mathbf{\gamma }}_2\left(t\right)/{\mathit{dt}}^{\beta },0<\alpha ,\beta <1$ , and the total deformation is given by the sum of the deformation obtained for each spring-pot, $\mathbf{\gamma }\left(t\right)={\mathbf{\gamma }}_1\left(t\right)+{\mathbf{\gamma }}_2\left(t\right)$ . The FMM can then be written as

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

This model allows a much better fit of rheological data, as shown in [12] but it is not frame invariant. However, following the same procedure employed with the Maxwell and K-BKZ model, that is, using the derivative of the relaxation function obtained for the Maxwell model as the memory function of the K-BKZ model, one can also use the derivative of the relaxation function of the FMM and insert it in the K-BKZ model, thus, obtaining a frame-invariant constitutive model, that retains all the good fitting properties of the FMM.

The relaxation function of the FMM can be obtained by solving the fractional differential Eq. (14) considering a constant deformation $\gamma ={\gamma }_{0}H\left(t\right)$ ( $H\left(t\right)$ is the Heaviside function) together with $\sigma \left(t_0\right)={\sigma }_0$ , leading to the relaxation modulus $G\left(t\right)=\sigma \left(t\right)/{\mathrm{\gamma }}_0$ given by:

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

where $E_{a,b}\left(z\right)$ is the generalized Mittag-Leffler function [7],

$$E_{\alpha ,\beta }\left(z\right)=\sum _{k=0}^{\infty }\frac{z^k}{\mathrm{\Gamma }\left(\mathit{\alpha k}+\beta \right)},$$ E16

and a characteristic measure of the relaxation spectrum described by the two spring-pots in series is $\lambda ={\left(\mathbb{V}/\mathbb{G}\right)}^{1/\left(\alpha -\beta \right)}$ .

This leads to the fractional K-BKZ model proposed by Jaishankar and Mckinley [12, 13], with $m\left(t-t^{\prime }\right)$ the memory function [2] in Eq. (6) now given by.

$$m\left(t-t^{\prime }\right)=\frac{\mathit{dG}\left(t-t^{\prime }\right)}{d{t}^{\prime }}=-\mathbb{G}{\left(t-t^{\prime }\right)}^{-1-\beta }{E}_{\alpha -\beta ,-\beta }\left(-\frac{\mathbb{G}}{\mathbb{V}}{\left(t-t^{\prime }\right)}^{\alpha -\beta }\right),$$ E17

Note that here the relaxation modulus $G\left(t-t^{\prime }\right)$ is the one obtained for the FMM. Please see [11, 12, 13] for more details. It should be remarked that the Mittag-Leffler function was used in the past by Guy Berry to describe polymeric materials exhibiting Andrade creep [14].

The fractional K-BKZ model is therefore given by:

$$\mathbf{\sigma }\left(t\right)=-\mathbb{G}\underset{0}{\overset{t}{\int }}{\left(t-t^{\prime }\right)}^{-1-\beta }{E}_{\alpha -\beta ,-\beta }\left(-\frac{\mathbb{G}}{\mathbb{V}}{\left(t-t^{\prime }\right)}^{\alpha -\beta }\right)h\left(I_1,I_2\right)\left({\mathbf{C}}_{t^{\prime }}^{-1}-\mathbf{I}\right)d{t}^{\prime },$$ E18

and we need to ensure that the integral converges (see the Foundations of Linear Viscoelasticity by Coleman and Noll [15]). The main problem seems to be the term ${\left(t-t^{\prime }\right)}^{-1-\beta }$ that diverges as $t^{\prime }\to t$ , and ${\int }_0^t{\left(t-t^{\prime }\right)}^{-1-\beta }d{t}^{\prime }$ diverges. Also, the term $E_{\alpha -\beta ,-\beta }\left(\dots \right)h\left(I_1,I_2\right)$ is finite $\forall t^{\prime },t^{\prime }\le t$ . Therefore, we must have $\left({\mathbf{C}}_{t^{\prime }}^{-1}-\mathbf{I}\right)=O\left({\left(t-t^{\prime }\right)}^m\right),m\ge 1$ as $t^{\prime }\to t$ so that ${\left(t-t^{\prime }\right)}^{-1-\beta }\left({\mathbf{C}}_{t^{\prime }}^{-1}-\mathbf{I}\right)=O\left({\left(t-t^{\prime }\right)}^n\right),n\le 1$ and therefore the integral converges.

It can be easily shown [1] that a Taylor series expansion of ${\mathbf{C}}_{t^{\prime }}^{-1}-\mathbf{I}$ about $t^{\prime }=t$ leads to.

$$\left({\mathbf{C}}_{t^{\prime }}^{-1}-\mathbf{I}\right)=-\sum _{k=1}^{\infty }\frac{{\left(t^{\prime }-t\right)}^k}{k!}{A}_k\left(t\right),$$ E19

with $A_k\left(t\right)$ the Rivlin-Ericksen tensors. Note that these tensors can be obtained directly from the velocity field without having to find the strain tensor [16]. We may therefore conclude that the integral is convergent, assuming a smooth velocity field is provided/obtained. Note that this does not mean that convergence problems will not arise during numerical calculations.

In Refs. [11, 12, 17], the beneficial fitting qualities of this constitutive model framework are discussed in detail. Here, we are interested in determining to what extent the properties of the Mittag-Leffler function can be used to improve the fitting quality of differential models, and this will be discussed in the next subsection.

3.2. Generalized Phan-Thien and Tanner model

The previous integral model given by Eq. (18) allows a good fit to experimental rheological data, in flows with defined kinematics where ${\mathbf{C}}^{-1}$ can be computed explicitly; but, it would be desirable to obtain also an improved frame-invariant differential model, that is easier to handle both mathematically and numerically, when compared to integral models, for solving complex flow problems with spatially varying kinematics. The model to be presented was recently proposed by our research group [18], and basically takes advantage of the flexible functional form of the Mittag-Leffler function by inserting this function into the already well-known Phan-Thien and Tanner (PTT) model, replacing the classical linear and exponential functions of the trace of the stress tensor.

The original exponential PTT model [19, 20] is given by.

$$\text{exp}\left(\frac{\mathit{\epsilon \lambda }}{\eta }{\mathbf{\sigma }}_{\mathit{kk}}\right)\mathbf{\sigma }+\lambda {\displaystyle \stackrel{\square }{\mathbf{\sigma }}}=\eta \dot{\mathbf{\gamma }},$$ E20

with ${\displaystyle \stackrel{\square }{\mathbf{\sigma }}}=\partial \mathbf{\sigma }/\partial t+\mathbf{u}\cdot \nabla \mathbf{\sigma }-{\left(\nabla \mathbf{u}\right)}^T\cdot \mathbf{\sigma }-\mathbf{\sigma }\cdot \nabla \mathbf{u}+\xi \left(\mathbf{D}\cdot \mathbf{\sigma }-\mathbf{\sigma }\cdot \mathbf{D}\right)$ being the Gordon-Schowalter derivative and ${\mathbf{\sigma }}_{\mathit{kk}}$ the trace of the stress tensor. Here, the parameter ξ accounts for slip between the molecular network and points in the continuous medium. The model was derived from a Lodge-Yamamoto type of network theory for polymeric fluids, in which the network junctions are not assumed to move strictly as points of the continuum but instead they are allowed a certain effective slip as well as a rate of destruction that depends on the state of stress in the network. Phan-Thien proposed that an exponential function form would be quite adequate to represent the rate of destruction of junctions and in [17] it was shown that the Mittag-Leffler function could improve the quality of model fits to real data by allowing different forms for the rates of destruction.

The model is then given by

$$\mathrm{\Gamma }\left(\beta \right)E_{\alpha ,\beta }\left(\frac{\mathit{\epsilon \lambda }}{\eta }{\mathbf{\sigma }}_{\mathit{kk}}\right)\mathbf{\sigma }+\lambda {\displaystyle \stackrel{\square }{\mathbf{\sigma }}}=\eta \dot{\mathbf{\gamma }},$$ E21

where the factor Γ(β) is used to ensure that $\mathrm{\Gamma }\left(\beta \right)E_{\alpha ,\beta }\left(0\right)=1$ .

This new model can further improve the accuracy of the description of real data obtained with the original exponential function of the trace of the stress tensor, as shown in [18].

4. Parametric study of the GPTT model

We will now present a detailed parametric study on the influence of the new parameters $\alpha ,\beta$ (arising from the Mittag-Leffler function) on the rheological behavior of the generalized exponential PTT model.

4.1. Steady-state shear flows

As shown in [18], the steady shear viscosity is given by $\eta \left(\dot{\gamma }\right)={\sigma }_{\mathit{xy}}\left(\dot{\gamma }\right)/\dot{\gamma }$ with

$${\sigma }_{\mathit{xy}}\left(\dot{\gamma }\right)=\frac{\eta \dot{\gamma }-{\sigma }_{\mathit{xx}}\mathit{Wi\xi }}{\mathrm{\Gamma }\left(\beta \right)E_{\alpha ,\beta }\left(\frac{\mathit{\epsilon \lambda }}{\eta }\left(\frac{2-2\xi }{2-\xi }\right){\sigma }_{\mathit{xx}}\right)},$$ E22

and ${\sigma }_{\mathit{xx}}$ is given by the solution of

$${\mathrm{\Gamma }}^2\left(\beta \right)E_{\alpha ,\beta }{\left(\frac{\mathit{\epsilon \lambda }}{\eta }\left(\frac{2-2\xi }{2-\xi }\right){\sigma }_{\mathit{xx}}\right)}^2{\sigma }_{\mathit{xx}}=\left(2-\xi \right){\left(\lambda \dot{\gamma }\right)}^2\left[\frac{\eta }{\lambda }-{\sigma }_{\mathit{xx}}\xi \right].$$ E23

Here $\mathit{Wi}=\lambda \dot{\gamma }$ is the dimensionless strength of the shear flow and $\eta ,\lambda ,\epsilon ,\xi ,\alpha ,\beta$ are the constitutive parameters of the generalized PTT (or GPTT) model.

Since we consider a simple plane shear flow aligned with the x-axis, we have that ${\sigma }_{\mathit{yy}}\left(\dot{\gamma }\right)={\sigma }_{\mathit{xx}}\xi /\left(2-\xi \right)$ and ${\sigma }_{\mathit{zz}}\left(\dot{\gamma }\right)=0$ (see [18] for more details).

Eqs. (22) and (23) can readily be solved using the Newton-Raphson method (solving first Eq. (23) and then substituting the numerical values obtained for ${\sigma }_{\mathit{xx}}$ into Eq. (22)).

Figure 6 shows the dimensionless steady shear viscosity obtained for the different parameters of the Mittag-Leffler function, $\alpha ,\beta$ . On the left, we show the influence of $\alpha$ by keeping constant all other parameters. On the right, we show the influence of $\beta$ (it should be remarked that when $\alpha ,\beta =1$ the exponential PTT model is recovered). We observe that when compared to the classical exponential PTT model, when $\alpha ,\beta <1$ , shear-thinning occurs for lower dimensionless shear rates and when $\alpha ,\beta >1$ there is a delay in the shear-thinning effect. For $\alpha ,\beta >1$ the shear viscosity increases, especially for high shear rates. Also, when we increase $\alpha$ , the slope of the shear viscosity curve for high dimensionless shear rates decreases (observed in Figure 6(a)), while varying $\beta$ , the slope seems to be the same, but a higher viscosity is obtained (observed in Figure 6(b)).

Figure 6.

Dimensionless shear viscosity obtained for the different parameters of the Mittag-Leffler function; (a) varying $\alpha$ holding the other five parameters constant and (b) varying $\beta$ .

Figure 7(a) shows the dimensionless steady shear viscosity, now obtained for different values of $\alpha ,\beta$ and $\epsilon$ . These plots allow one to see that the $\epsilon$ parameter may not be compared directly to the value used in the classical models (featuring linear and exponential functions of the trace of the stress tensor). For comparison purposes, we plot again the curve obtained for the exponential PTT model with $\epsilon =0.25$ ( $\alpha =\beta =1$ ) by the dash-dot lines.

Note that (see Figure 7(b)) small variations of the parameter $\epsilon$ allows one to control the rate of transition to the shear-thinning at high $\mathit{Wi}$ while maintaining a similar shear thinning set point.

Figure 7.

Dimensionless shear viscosity obtained for the different parameters of the Mittag-Leffler function varying: (a) Constant $\epsilon$ ; (b) Varying $\epsilon$ .

Figure 7 shows that by setting different combinations of $\alpha ,\beta$ we may obtain different slopes at higher dimensional shear rates. For low $\beta$ and high $\alpha$ , we obtain a lower slope but a premature shear viscosity thinning, while for high $\beta$ and low $\alpha$ , we obtain a higher slope but a delayed shear-thinning.

4.2. Steady-state elongational flows

The steady unidirectional extensional viscosity is defined as ${\eta }_E=\left({\sigma }_{\mathit{xx}}-{\sigma }_{\mathit{yy}}\right)/\dot{\varsigma }$ , where $\dot{\varsigma }$ is the imposed elongation rate [18], and can be obtained by solving the system of equations (for a simpler technique that does not involve an iterative procedure, please consult [18])

$${\sigma }_{\mathit{xx}}\left(\mathrm{\Gamma }\left(\beta \right)E_{\alpha ,\beta }\left(\frac{\mathit{\epsilon \lambda }}{\eta }{\sigma }_{\mathit{kk}}\right)-2\lambda \dot{\varsigma }\left(1-\xi \right)\right)=2\eta \dot{\varsigma },$$ E24

$${\sigma }_{\mathit{xx}}\left(\mathrm{\Gamma }\left(\beta \right)E_{\alpha ,\beta }\left(\frac{\mathit{\epsilon \lambda }}{\eta }{\sigma }_{\mathit{kk}}\right)+\lambda \dot{\varsigma }\left(1-\xi \right)\right)=-\eta \dot{\varsigma },$$ E25

with ${\sigma }_{\mathit{kk}}={\sigma }_{\mathit{xx}}+2{\sigma }_{\mathit{yy}}$ .

Figure 8 shows the dimensionless steady elongational viscosity obtained for different parameters of the Mittag-Leffler function. In Figure 8(a), we show the influence of $\alpha$ by keeping constant all other parameters. In Figure 8(b), we show the influence of $\beta$ . Note that we have used the same parameters as in the shear viscosity case.

Figure 8.

Dimensionless elongational viscosity obtained for different parameters of the Mittag-Leffler function: (a) Varying α; (b) Varying β.

Note that when we increase $\alpha ,\beta$ , we observe an increase of the elongational viscosity, with the maximum value being reached for higher dimensionless extensional rates. Again, we observe different asymptotic slopes for high extension rates (when varying $\alpha$ ). Note that there is no overshoot for low values of $\beta$ .

We may conclude that by varying $\alpha ,\beta$ , we change both the shear and elongational viscosities, and therefore the fit to experimental data should be performed with care, taking into account this dependence.

Figure 9 shows the effect of the parameters used in Figure 7, for the case of elongational viscosity. The results are qualitatively similar to the ones obtained in Figure 7, that is, in terms of changes to the asymptotic slopes at high deformation rates and premature/delayed thinning. It can be observed that the elongational viscosity is more sensitive to changes in the parameters $\alpha ,\beta$ and $\epsilon$ . This result is to be expected since this is a strong flow, and, the exponential PTT model was originally proposed to be able to describe the response of complex fluids in strong flows. Figure 9(a) shows that the overshoot can be suppressed using a low $\beta$ and high $\alpha$ values. Note that the maximum extensional viscosity is obtained for the exponential PTT model, and that the values of $\alpha ,\beta$ have a strong influence on the asymptotic slope of the curve for high extensional rates. Figure 9(b) shows three different curves for different combinations of $\alpha ,\beta$ and $\epsilon$ . Note that for $\alpha =0.1,\beta =0.1$ and $\epsilon =10$ we can also suppress the overshoot in the extensional viscosity, and for $\alpha =2,\beta =0.1$ and $\epsilon =0.05$ we can decrease the curvature of the overshoot, and at the same time decrease the slope of curve.

Figure 9.

Dimensionless elongational viscosity obtained for the different parameters of the Mittag-Leffler function: (a) Constant $\epsilon$ ; (b) Varying $\epsilon$ .

4.3. Steady-state shear and elongational flows

Until now, we have explored generally the influence of the different model parameters on the behavior of the GPTT model for steady flows, but, a more quantitative side-by-side comparison between the shear and elongational flow curves was not performed, and the limited flexibility of the classical exponential PTT model for fitting experimental data (when compared to the GPTT) was not explored. In Figure 10, we try to illustrate the advantages of using the Mittag-Leffler function instead of the classical exponential one. To this end, we present the viscometric predictions obtained for both shear and elongational flows for both models (GPTT and exponential PTT).

Figure 10.

Comparison of the dimensionless elongational and shear viscosity obtained for different parameters of the Mittag-Leffler function, varying $\epsilon$ , and the classical exponential PTT model (α = 1, β = 0).

Figure 10 illustrates the additional flexibility of using the Mittag-Leffler function, by showing that we can manipulate the magnitude of the increase in the elongational viscosity and at the same time only slightly change the shear viscosity. This allows better fits to rheological data when using the Mittag-Leffler function [18]. Note that in the exponential PTT model, when we increase the $\epsilon$ parameter, both the shear and elongational viscosities increase concomitantly.

5. Conclusions

In this chapter, we have presented a brief introduction to the world of viscoelastic models capable of describing the rheology of complex fluids, and we have summarized some of the well-known classical differential and integral models.

With incorporation of ideas from fractional calculus, most of these models can be further improved, either by changing classical local operators for improved (non-local) fractional versions, or, either using new analytic functions that arise in the realm of fractional differential equations, such as the Mittag-Leffler function.

As an example, we present the fractional K-BKZ model and the recently proposed generalized PTT model. The fractional K-BKZ model allows a better description of fluid flow behavior (when compared to the generalized PTT model), but, increases the need for high computational power. Therefore, the novelty of the present work is our detailed study on the influence of the Mittag-Leffler function in shear and elongational flows of a generalized PTT model.

Acknowledgments

L.L. Ferrás and A.M. Afonso acknowledge the Project PTDC/EMS-ENE/3362/2014-POCI-010145-FEDER-016665—funded by FEDER funds through COMPETE2020—Programa Operacional Competitividade e Internacionalização (POCI) and by national funds through FCT—Fundação para a Ciência e a Tecnologia; I.P. L.L. Ferrás would also like to thank the funding by FCT through the scholarship SFRH/BPD/100353/2014. M.L. Morgado would like to thank the funding by FCT through Project UID/MULTI/04621/2013 and M. Rebelo would also like to thank the funding by FCT through Project UID/MAT/00297/2013 (Centro de Matemática e Aplicações).