Open access peer-reviewed chapter

# Recent Advances in Complex Fluids Modeling

By Luís L. Ferrás, Maria L. Morgado, Magda Rebelo, Rosalía T. Leiva, António Castelo, Gareth H. McKinley and Alexandre M. Afonso

Submitted: March 1st 2018Reviewed: November 23rd 2018Published: May 29th 2019

DOI: 10.5772/intechopen.82689

## 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.

### Keywords

• Mittag-Leffler
• viscoelastic
• memory function
• fractional calculus
• rheology

## 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(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 Fneeded to extend or compress a spring by some distance γe=xx0is proportional to that distance, F=kγ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.

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):

ForceArea=σ=ηUh=ηduydyσ=ηduydyE1

with σthe unidirectional stress and η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 γ̇.

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, σ=ηdγft/dt, and a spring, σ=Gγet(with the subscripts f and e standing for Newtonian fluid and Hookean elastic solid, respectively), as shown in Figure 3.

The total deformation γis the sum of the deformation obtained from the spring γeand the dashpot γf, and the rate of deformation is given by:

tdt=dγftdt+dγetdttdt=ση+1Gdtσ+λdt=ηtdtMaxwell  Model,λ=ηGE2

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:

σ+λdσdt=ηdγtdtE3

with σthe stress tensor, γ̇=u+uTthe rate of deformation tensor, uthe velocity vector, λthe relaxation time of the fluid and ηthe zero shear rate viscosity. This model can be equivalently written in integral form as

σt=0tGett/λdγdtdt,E4

where G=η/λ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 σ+λσ̂=ηγ̇(with σ̂=σ/t+uσuTσσuthe upper-convected derivative) that can also be written in integral form as

σt=0tηλ2ett/λCt1IdtE5

where C1is the Finger strain tensor (a frame-invariant measure of deformation) [1]. The term mtt=η/λ2ett/λis known as the memory function (the derivative of the relaxation modulus Gtt). Note that the relaxation modulus can be easily obtained by imposing a step strain (constant deformation), as shown in Figure 4, resulting in Gt=σ/γ0=Gett/λ.

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],

σt=0tmtthI1I2Ct1Idt,E6

where C1is the Finger tensor [1], I1,I2are the traces of C1and C, respectively, and hI1I2is 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 mttwas proposed to be of the form:

mtt=aλett/λ,E7

where aand λ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=η/λand hI1I2=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 iaiett/λ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).

Imagine a portion of material that is principally formed of two different regions. In these regions, two similar physical processes φ1and φ2occur (for the time being it does not matter what is the process under study), but, at different rates, dφ1/dt=0.1for Region I and dφ2/dt=1for Region II. If we look at the portion of material as a whole, one would naturally choose the rate of 1as 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βφ/dtβwith β=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 ftis given by the formula

atatatftdtdtdt=ntimesJanft=1n1!atttn1ftdt.E8

A generalization to non-integer values of ncan be performed using the Euler Gamma function Γx, leading to the Riemann-Liouville fractional integral

Janft=1Γαatttα1ftdt,E9

where we have used αto represent the generalization of nto non-integer values. A fractional derivative of any order can then be obtained by manipulating the number of integrations and differentiations of the function ft. By performing the mα-fold integration of the mthderivative of ft, JamαDmftwith m=α, we arrive at the generalized derivative formula (Caputo fractional derivative [6]) of order m1<α<m,

dαftdtα=1Γmαatttα+m1dmftdtmdt,m1<α<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 Gttα,0<α<1[8]. For example, the critical gel model investigated by Winter and Chambon is written Gt=Stα. If we assume the relaxation modulus for an arbitrary loading history in such materials is given by Gtt=VΓ1α1ttα(Vis known as a quasi-property [9] and is connected to the critical gel strength by S=V/Γ1α), then we have that:

σt=1Γ1α0tVttαdγdtdt.E11

By recognizing that the Caputo fractional derivative of a general function γt(in our case γtis the deformation) is defined as [10]:

dαγtdtα=1Γ1α0tttαdγdtdt,E12

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

σ=Vdαγtdtα,0<α<1,E13

This model provides a generalized viscoelastic response, in the sense that when α=1we obtain a Newtonian fluid, and when α=0we 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 α1because of the singularity in Γ1α[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, σ=Vdαγ1t/dtα=Gdβγ2t/dtβ,0<α,β<1, and the total deformation is given by the sum of the deformation obtained for each spring-pot, γt=γ1t+γ2t. The FMM can then be written as

σt+VGdαβσtdtαβ=Vdαγtdtα,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 γ=γ0Ht(Htis the Heaviside function) together with σt0=σ0, leading to the relaxation modulus Gt=σt/γ0given by:

Gt=GtβEαβ,1βGVtαβ,E15

where Ea,bzis the generalized Mittag-Leffler function [7],

Eα,βz=k=0zkΓαk+β,E16

and a characteristic measure of the relaxation spectrum described by the two spring-pots in series is λ=V/G1/αβ.

This leads to the fractional K-BKZ model proposed by Jaishankar and Mckinley [12, 13], with mttthe memory function [2] in Eq. (6) now given by.

mtt=dGttdt=Gtt1βEαβ,βGVttαβ,E17

Note that here the relaxation modulus Gttis 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:

σt=G0ttt1βEαβ,βGVttαβhI1I2Ct1Idt,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 tt1βthat diverges as tt, and 0ttt1βdtdiverges. Also, the termEαβ,βhI1I2is finite t,tt. Therefore, we must have Ct1I=Ottm,m1as ttso that tt1βCt1I=Ottn,n1and therefore the integral converges.

It can be easily shown [1] that a Taylor series expansion of Ct1Iabout t=tleads to.

Ct1I=k=1ttkk!Akt,E19

with Aktthe 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 C1can 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.

expελησkkσ+λσ=ηγ̇,E20

with σ=σ/t+uσuTσσu+ξDσσDbeing the Gordon-Schowalter derivative and σkkthe 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

ΓβEα,βελησkkσ+λσ=ηγ̇,E21

where the factor Γ(β) is used to ensure that ΓβEα,β0=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 α,β(arising from the Mittag-Leffler function) on the rheological behavior of the generalized exponential PTT model.

As shown in [18], the steady shear viscosity is given by ηγ̇=σxyγ̇/γ̇with

σxyγ̇=ηγ̇σxxWiξΓβEα,βελη22ξ2ξσxx,E22

and σxxis given by the solution of

Γ2βEα,βελη22ξ2ξσxx2σxx=2ξλγ̇2ηλσxxξ.E23

Here Wi=λγ̇is the dimensionless strength of the shear flow and η,λ,ε,ξ,α,β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 σyyγ̇=σxxξ/2ξand σzzγ̇=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 σxxinto Eq. (22)).

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

Figure 7(a) shows the dimensionless steady shear viscosity, now obtained for different values of α,βand ε. These plots allow one to see that the ε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 ε=0.25(α=β=1) by the dash-dot lines.

Note that (see Figure 7(b)) small variations of the parameter εallows one to control the rate of transition to the shear-thinning at high Wiwhile maintaining a similar shear thinning set point.

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

The steady unidirectional extensional viscosity is defined as ηE=σxxσyy/ς̇, where ς̇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])

σxxΓβEα,βελησkk2λς̇1ξ=2ης̇,E24
σxxΓβEα,βελησkk+λς̇1ξ=ης̇,E25

with σkk=σxx+2σ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 αby keeping constant all other parameters. In Figure 8(b), we show the influence of β. Note that we have used the same parameters as in the shear viscosity case.

Note that when we increase α,β, 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 α). Note that there is no overshoot for low values of β.

We may conclude that by varying α,β, 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 α,βand ε. 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 βand high αvalues. Note that the maximum extensional viscosity is obtained for the exponential PTT model, and that the values of α,β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 α,βand ε. Note that for α=0.1,β=0.1and ε=10we can also suppress the overshoot in the extensional viscosity, and for α=2,β=0.1and ε=0.05we can decrease the curvature of the overshoot, and at the same time decrease the slope of curve.

### 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 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 ε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).

## How to cite and reference

### Cite this chapter Copy to clipboard

Luís L. Ferrás, Maria L. Morgado, Magda Rebelo, Rosalía T. Leiva, António Castelo, Gareth H. McKinley and Alexandre M. Afonso (May 29th 2019). Recent Advances in Complex Fluids Modeling, Fluid Flow Problems, Farhad Ali, IntechOpen, DOI: 10.5772/intechopen.82689. Available from:

### Related Content

Next chapter

#### Super-Speeding Jets in MHD Couette Flow

By Krzysztof Mizerski

#### Pattern Formation and Stability in Magnetic Colloids

Edited by Nicolás O. Rojas

First chapter

#### Introductory Chapter: An Example in Superparamagnetic Colloids

By Nicolás O. Rojas

We are IntechOpen, the world's leading publisher of Open Access books. Built by scientists, for scientists. Our readership spans scientists, professors, researchers, librarians, and students, as well as business professionals. We share our knowledge and peer-reveiwed research papers with libraries, scientific and engineering societies, and also work with corporate R&D departments and government entities.

View all Books