Chapter 11 Complex Fluids and Non-Newtonian Rheology

To begin with, let’s recap what we know for Newtonian fluids. The fluid is Newtonian if the relationship between the rate of deformation (${\gamma }_{ij}=\partial u_i/\partial x_j+\partial u_j/\partial x_i$) and the stress $𝑻$ is local, linear, instantaneous and isotropic. If the fluid is also incompressible, the Newtonian constitutive law is

$𝑻=-p𝑰+\mu \left(\mathbf{\nabla }𝒖+\mathbf{\nabla }{𝒖}^T\right)\equiv -p𝑰+\underset{𝝉}{\underbrace{\mu \dot{𝜸}}},$

(11.2)

where $\mu$ is constant and $𝝉$ is the deviatoric stress tensor and represents the traceless ($\mathbf{\nabla }\cdot 𝒖=0$) part of the stress tensor, while the pressure is the isotropic part.

A generalised Newtonian fluid assumes that the applied flow only changes the energy dissipation in the fluid, represented by the viscosity, but does not change the tensorial structure of the Newtonian model. The generalised constitutive law is given by

$𝑻=-p𝑰+\mu (\dot{𝜸})\dot{𝜸},$

(11.3)

where $\mu (\dot{𝜸})$ is the viscosity now dependent on the fluid flow. To ensure that the viscosity does not change under a coordinate transformation, which would be unphysical, the viscosity can only depend on invariants of the strain-rate tensor $\dot{𝜸}$. The second tensorial invariant is the lowest non-trivial invariant, so we may write

$𝑻=-p𝑰+\mu (\dot{\gamma })\dot{𝜸},\dot{\gamma }=\sqrt{\frac{1}{2}{\mathrm{\Sigma }}_{jk}{\dot{\gamma }}_{jk}^2}.$

(11.4)

For example, a Bingham fluid has constitutive law

$𝝉=\left({\displaystyle \frac{{\tau }_Y}{\dot{\gamma }}}+\mu \right)\dot{𝜸},\tau =\sqrt{\frac{1}{2}{\mathrm{\Sigma }}_{jk}{\tau }_{jk}^2}>{\tau }_Y,$

(11.7)

and $\dot{𝜸}=\mathrm{𝟎}$ otherwise, where $\mu$ is a constant viscosity. Note, $\tau$ is the second invariant of the deviatoric stress tensor $𝝉$. This is chosen to decide whether the fluid is yielded or not as it does not change under a coordinate transformation and the first invariant is zero.

A Herschel-Bulkley fluid has constitutive law

$𝝉=\left({\displaystyle \frac{{\tau }_Y}{\dot{\gamma }}}+K{|\dot{\gamma }|}^{n-1}\right)\dot{𝜸},\tau =\sqrt{\frac{1}{2}{\mathrm{\Sigma }}_{jk}{\tau }_{jk}^2}>{\tau }_Y,$

(11.8)

and $\dot{𝜸}=\mathrm{𝟎}$ otherwise, where $K$ is the consistency with units $P{T}^n$ and $n$ is the power-law index.

Let’s first look at a simple model of a viscoelastic fluid called a Maxwell fluid. Consider a shear deformation where adjacent layers of the material are shifted in the same direction along the planes relative to each other. The deformation can be described by the strain $\gamma$, which is approximately the total relative shift divided by the distance between the two layers assuming the displacement is small.

The shear stress $\tau$ for an elastic solid is given by Hooke’s law

${\tau }_e=G{\gamma }_e,$

(11.39)

where $G$ is the shear modulus of the material. For a Newtonian fluid the constitutive law is

${\tau }_v={\mu }_p{\dot{\gamma }}_v,$

(11.40)

where ${\mu }_p$ is the constant Newtonian viscosity. Graphically we could depict these two elements as a spring and a dampner. Depending on how we connect them, we form a viscoelastic fluid or viscoelastic solid.

The 1D Maxwell model $\tau +\lambda \dot{\tau }=\mu \dot{\gamma }$, where $\lambda =\mu /G$, can be solved to give

${\displaystyle \frac{d}{dt}}\left(e^{t/\lambda }\tau \right)={\displaystyle \frac{\mu \dot{\gamma }}{\lambda }}{e}^{t/\lambda }\Rightarrow \tau (t)={\displaystyle \frac{\mu }{\lambda }}{\displaystyle {\int }_{-\mathrm{\infty }}^t}\dot{\gamma }(t^{\prime })e^{-(t-t^{\prime })/\lambda }d{t}^{\prime }.$

(11.43)

The parameter $\lambda$ is the Maxwell relaxation time. Looking at the expression for the stress response, we see that the stress created by a stepwise deformation relaxes exponentially on a time scale $\lambda$ (viscous-fluid-like property). Whereas, on short times $\tau (t)\sim \mu \gamma (t)/\lambda$ and the material behaves like a solid. For example, consider the stress response to a periodic deformation with frequency $\omega$, where $\gamma ={\gamma }_0\sin (\omega t),\dot{\gamma }={\gamma }_0\omega \cos (\omega t)$. Substituting into the integral expression gives

	$\tau (t)$	$={\displaystyle \frac{\mu {\gamma }_0\omega }{\lambda }}{\displaystyle {\int }_{-\mathrm{\infty }}^t}\cos (\omega t^{\prime })e^{-(t-t^{\prime })/\lambda }d{t}^{\prime }=\mu {\gamma }_0\omega {\displaystyle \frac{\cos (\omega t)+\lambda \omega \sin (\omega t)}{1+{(\lambda \omega )}^2}}$		(11.44)
		$={\displaystyle \frac{\mu \dot{\gamma }}{1+{(\lambda \omega )}^2}}+{\displaystyle \frac{G{(\lambda \omega )}^2\gamma }{1+{(\lambda \omega )}^2}}=\tilde{\mu }(\omega )\dot{\gamma }+\tilde{G}(\omega )\gamma ,$		(11.45)

where $\tilde{\mu }(\omega ),\tilde{G}(\omega )$ are the frequency dependent viscosity and shear modulus. At short times ($\lambda \omega \gg 1$), $\tilde{G}\simeq G$ and the material behaves like a solid. Whereas, at long times ($\lambda \omega \ll 1$), $\tilde{\mu }\simeq \mu$ and the material behaves as a viscous fluid. The crossover between these two behaves occurs when the time scale of deformation is the same as the time scale of relaxation, ${\omega }^{-1}\sim \lambda$.

One version of a frame invariant derivative is the upper convected derivative defined as

${\stackrel{▽}{𝝉}}_p={\displaystyle \frac{\partial {𝝉}_p}{\partial t}}+𝒖\cdot \mathbf{\nabla }{𝝉}_p-{(\mathbf{\nabla }𝒖)}^T\cdot {𝝉}_p-{𝝉}_p\cdot (\mathbf{\nabla }𝒖),$

(11.48)

which transforms as a covariant tensor. N.B. there are other choices for a frame invariant derivative but we will only consider the upper convected derivative during this course. Rewriting the relationship between stress and strain-rate (11.46) using the upper convected derivative gives

${𝝉}_p+{\displaystyle \frac{{\mu }_p}{G}}{\stackrel{▽}{𝝉}}_p={\mu }_p\dot{𝜸}$

(11.49)

to describe the polymeric stress which incorporates both viscous and elastic properties. The relaxation time of the fluid is defined as $\lambda ={\mu }_p/G$. For a polymeric fluid, there is an additional Newtonian contribution from the solvent viscosity ${\mu }_s$. Hence, the total deviatoric stress tensor is written as

$𝝉={\mu }_s\dot{𝜸}+{𝝉}_p.$

(11.50)

This model is known as the Oldroyd B model. In the limit of vanishing deformation rates it describes a Newtonian fluid with viscosity $\mu ={\mu }_s+{\mu }_p$.

 Oldroyd B Model:

	$\mathbf{\nabla }\cdot 𝒖=0$		(11.51)
	$\rho \left({\displaystyle \frac{\partial 𝒖}{\partial t}}+𝒖\cdot \mathbf{\nabla }𝒖\right)=\mathbf{\nabla }\cdot 𝑻$		(11.52)
	$𝑻=-p𝑰+𝝉=-p^{*}𝑰+{\mu }_s\dot{𝜸}+G𝑨$		(11.53)
	$\stackrel{▽}{𝑨}={\displaystyle \frac{\partial 𝑨}{\partial t}}+𝒖\cdot \mathbf{\nabla }𝑨-{(\mathbf{\nabla }𝒖)}^T\cdot 𝑨-𝑨\cdot (\mathbf{\nabla }𝒖)=-{\displaystyle \frac{1}{\lambda }}(𝑨-𝑰)$		(11.54)

where

$p^{*}=p+{\displaystyle \frac{{\mu }_p}{\lambda }},G={\displaystyle \frac{{\mu }_p}{\lambda }}\mathrm{and}𝑨={\displaystyle \frac{1}{G}}\left[𝝉-{\mu }_s\dot{𝜸}+G𝑰\right].$

(11.55)

 Here, we have rewritten the Oldroyd B model in terms of variable $𝑨$; a symmetric second rank tensor. This form of the model shows that $𝑨$ tends towards its relaxed state $𝑨=𝑰$ at long times governed by long timescale $\lambda$ at which point the polymer stress ${𝝉}_p$ vanishes. To understand the properties of Oldroyd B fluids let’s see how it behaves in some simple fluid flows.

First consider a shear flow of the form

$𝒖=(\dot{\gamma }y,0,0),$

(11.56)

which only depends on $t$ ($\dot{\gamma }=\dot{\gamma }(t)$) and $y$ so all variables must also only be functions of $t$ and $y$. By construction this satisfies mass conservation. The momentum equation reduces to

$\rho {\displaystyle \frac{\partial 𝒖}{\partial t}}=\mathbf{\nabla }\cdot 𝑻,$

(11.57)

where

(11.64)

and

(11.71)

Now let’s consider a two dimensional extensional flow of the form

$𝒖=(\dot{\epsilon }x,-\dot{\epsilon }y),$

(11.83)

where $\dot{\epsilon }$ is constant. Again, this automatically satisfies mass conservation. The strain-rate and stress tensors are

	$\dot{𝜸}=$			(11.86)
	$𝑻=$			(11.93)

with the components $A_{ij}$ satisfying

(11.100)

Rearranging, we find

$A_{xx}={\displaystyle \frac{1}{1-2\lambda \dot{\epsilon }}},A_{xy}=0,A_{yy}={\displaystyle \frac{1}{1+2\lambda \dot{\epsilon }}}.$

(11.101)

The total stress is then

(11.104)

where

${\mu }_{ext}={\mu }_s+{\displaystyle \frac{\lambda G}{1-4{\lambda }^2{\dot{\epsilon }}^2}}={\mu }_s+{\displaystyle \frac{{\mu }_p}{1-4{\lambda }^2{\dot{\epsilon }}^2}}$

(11.105)

If we plot this viscosity as a function of the strain rate $\dot{\epsilon }$ (figure 11.6) we see that initially the viscosity increases (hence thickens, or strain hardens). But when $\dot{\epsilon }=1/2\lambda$ the viscosity diverges and then is negative for larger values. This result is clearly unphysical. This happens because the Oldroyd B model is derived using Hooke’s law which is infinitely extensible. To get around this a modification, a nonlinear spring law must be used.

A weakly non-linear flow is a situation wherein the flow changes on time scales much longer than the relaxation time, $\lambda$ , and therefore $\lambda$ is in some sense small and can be used as an expansion parameter. It is generally a bad practice to perform a Taylor expansion in a dimensional variable; a better expansion parameter might be $\lambda /T_0$, where $T_0$ is the typical time scale set by the flow.

From the Oldroyd B model we have

$𝝉+\lambda \stackrel{▽}{𝝉}=\mu \left(\dot{𝜸}+\lambda {\displaystyle \frac{{\mu }_s}{\mu }}\stackrel{▽}{\dot{𝜸}}\right),$

(11.106)

where $\mu ={\mu }_s+{\mu }_p$. Let $P$ be a characteristic pressure scale, $L$ a characteristic length scale, $L/T_0$ a characteristic velocity scale. The non-dimensional version of equation (11.106) is

$𝝉+ϵ\stackrel{▽}{𝝉}={\displaystyle \frac{\mu }{P{T}_0}}\left(\dot{𝜸}+ϵ{\displaystyle \frac{{\mu }_s}{\mu }}\stackrel{▽}{\dot{𝜸}}\right),$

(11.107)

where $ϵ=\lambda /T_0$. Expanding in this parameter, we may write

$𝝉={𝝉}^{(0)}+ϵ{𝝉}^{(1)}+O({ϵ}^2).$

(11.108)

Substituting this into the Oldroyd B model, we recover a Newtonian constitutive law to leading order,

${𝝉}^{(0)}={\displaystyle \frac{\mu }{P{T}_0}}\dot{𝜸},\mathrm{and}{𝝉}^{(1)}=-{\stackrel{▽}{𝝉}}^{(0)}+{\displaystyle \frac{{\mu }_s}{P{T}_0}}\stackrel{▽}{\dot{𝜸}}=-{\displaystyle \frac{{\mu }_p}{P{T}_0}}\stackrel{▽}{\dot{𝜸}}.$

(11.109)

Putting the dimensions back in, we have

$𝝉=\mu \dot{𝜸}-\lambda {\mu }_p\stackrel{▽}{\dot{𝜸}}.$

(11.110)

This is an example of a second-order fluid.

 Second-order fluid: The simplest class of equations for viscoelastic solutions involves the expression of the stress tensor as a sum of all admissible combinations of the velocity gradient tensor. Depending on the highest algebraic power of the velocity gradient tensor involved, they are called the second-order fluid, third-order fluid, etc. For example, the deviatoric stress of the second-order fluid is given by

$𝝉=\mu \dot{𝜸}+b_2\stackrel{▽}{\dot{𝜸}}+b_{11}\dot{𝜸}\cdot \dot{𝜸},$

(11.111)

where $b_2,b_{11}$ are material constants.  

Returning back to our Oldroyd B model, the weakly non-linear expansion is a second-order fluid with material parameters $b_2=-\lambda {\mu }_p,b_{11}=0$.

For an Oldroyd B fluid we showed that there is a normal stress difference in a shear flow. In this section we will determine the steady free surface of a weakly non-linear second-order fluid caused by a rod of radius $r=a$ rotating at a frequency $\mathrm{\Omega }$ in the fluid. Flow is purely in the azimuthal direction so we write $𝐮=(0,v(r),0)$ in cylindrical polars.

The flow is weakly non-linear so to leading order, the fluid behaves like a Newtonian fluid where ${\dot{\gamma }}_{r\theta }$ is the only non-zero strain-rate component. In terms of the deviatoric stress components, taking $P$ to be a characteristic stress scale, ${\tau }_{r\theta }(r)\sim O(P)$ and ${\tau }_{rr}(r),{\tau }_{\theta \theta }(r)\sim O(ϵP)$ with pressure $p(r,z)\sim O(P)$. The leading order momentum equations give

$0={\displaystyle \frac{\partial }{\partial r}}\left(T_{r\theta }\right)+{\displaystyle \frac{2}{r}}{T}_{r\theta }={\displaystyle \frac{1}{r^2}}{\displaystyle \frac{\partial }{\partial r}}\left(r^2{T}_{r\theta }\right),0={\displaystyle \frac{\partial T_{zz}}{\partial z}}-\rho g.$

(11.112)

We can use these to solve for the azimuthal component of the velocity. The leading order r-momentum equation gives $\partial p/\partial r=0$. The r-momentum equation at order $O(\lambda /T_0)$ can be used to calculate the higher order components of the stress

$0={\displaystyle \frac{\partial T_{rr}}{\partial r}}+{\displaystyle \frac{T_{rr}-T_{\theta \theta }}{r}}.$

(11.113)

For the leading order balance, the fluid has a Newtonian constitutive law:

${\dot{\gamma }}_{r\theta }=r{\displaystyle \frac{d}{dr}}\left({\displaystyle \frac{v}{r}}\right),T_{r\theta }=\mu {\dot{\gamma }}_{r\theta }\Rightarrow T_{r\theta }={\displaystyle \frac{A}{r^2}}={\displaystyle \frac{d}{dr}}\left({\displaystyle \frac{v}{r}}\right)$

(11.114)

Integrating and applying boundary conditions $v(a)=\mathrm{\Omega }a,v\to 0$ as $r\to \mathrm{\infty }$, gives

$v={\displaystyle \frac{\mathrm{\Omega }{a}^2}{r}},{\dot{\gamma }}_{r\theta }=r{\displaystyle \frac{d}{dr}}\left({\displaystyle \frac{v}{r}}\right)=-{\displaystyle \frac{2\mathrm{\Omega }{a}^2}{r^2}}.$

(11.115)

The vertical component of the Stokes equation gives a hydrostatic balance

$0={\displaystyle \frac{\partial T_{zz}}{\partial z}}-\rho g\Rightarrow p=\rho g(h(r)-z).$

(11.116)

We now consider the second-order contribution to the stress. Using the second-order fluid model, we have

$T_{r\theta }=\mu {\dot{\gamma }}_{r\theta },T_{rr}=-p+b_{11}{\dot{\gamma }}_{r\theta }^2,T_{\theta \theta }=-p-2b_2{\dot{\gamma }}_{r\theta }^2+b_{11}{\dot{\gamma }}_{r\theta }^2,T_{zz}=-p.$

(11.117)

since ${(\stackrel{▽}{\dot{𝜸}})}_{rr}=0,{(\stackrel{▽}{\dot{𝜸}})}_{\theta \theta }=-2{\dot{\gamma }}_{r\theta }^2$. Substituting into the r-momentum balance gives

	${\displaystyle \frac{\partial T_{rr}}{\partial r}}+{\displaystyle \frac{2b_2{\dot{\gamma }}_{r\theta }^2}{r}}=0\Rightarrow -p+b_{11}{\dot{\gamma }}_{r\theta }^2={\displaystyle \frac{b_2}{2}}{\dot{\gamma }}_{r\theta }^2+f(z)$		(11.118)
	$-\rho g(h(r)-z)+(b_{11}-b_2/2){\dot{\gamma }}_{r\theta }^2=f(z).$		(11.119)

Choosing $f(z)=\rho gz$, we arrive at

$h(r)={\displaystyle \frac{1}{\rho g}}\left(b_{11}-{\displaystyle \frac{b_2}{2}}\right){\dot{\gamma }}_{r\theta }^2.$

(11.120)

Referring back to the material parameters for the Oldroyd B model, we see the surface has height

$h(r)={\displaystyle \frac{2\lambda {\mu }_p{\mathrm{\Omega }}^2{a}^4}{\rho g{r}^4}}.$

(11.121)

Hence, for an Oldroyd B model, rod climbing occurs.

1. 1.

   The linear Maxwell model is not objective
   The model is not frame invariant so none of the conclusions drawn from studying the model can be guaranteed to be physical.

2. 2.

   Time-dependent flows in weakly non-linear viscoelastic fluids are unstable
   Consider the second-order fluid model

   $𝝉=\mu \dot{𝜸}+b_2\stackrel{▽}{\dot{𝜸}}+b_{11}\dot{𝜸}\cdot \dot{𝜸}$

   (11.122)

   in a two-dimensional shear flow $𝐮=(u(y,t),0)$ in a channel $0\le y\le h$. From momentum balance we have

   $\rho {\displaystyle \frac{\partial u}{\partial t}}={\displaystyle \frac{\partial T_{xy}}{\partial y}}\equiv {\displaystyle \frac{\partial {\tau }_{xy}}{\partial y}}.$

   (11.123)

   From the constitutive law, the deviatoric stress ${\tau }_{xy}$ is

   ${\tau }_{xy}=\mu {\displaystyle \frac{\partial u}{\partial y}}+b_2{\displaystyle \frac{\partial }{\partial t}}{\displaystyle \frac{\partial u}{\partial y}}.$

   (11.124)

   These two can be combined to give

   $\rho {\displaystyle \frac{\partial u}{\partial t}}=\mu {\displaystyle \frac{{\partial }^{2}u}{\partial y^2}}+b_2{\displaystyle \frac{\partial }{\partial t}}{\displaystyle \frac{{\partial }^{2}u}{\partial y^2}}.$

   (11.125)

   Taking no-slip boundary conditions at the walls of a channel located at $y=0$ and $y=h$, we write the flow velocity as a Fourier sine series:

   $u(y,t)={\mathrm{\Sigma }}_{m=1}^{\mathrm{\infty }}{u}_m{e}^{{\alpha }_{m}t}\sin \left({\displaystyle \frac{m\pi y}{h}}\right).$

   (11.126)

   Substituting this into the above equation and rearranging, we find growth rates ${\alpha }_m$

   ${\alpha }_m=-{\displaystyle \frac{\mu }{b_2+\frac{\rho h^2}{{(m\pi )}^2}}}.$

   (11.127)

   For polymer solutions $b_2$ is typically negative (recall its value based on the Oldroyd B model, $b_2=-\lambda {\mu }_p$), so that ${\alpha }_m$ is positive for sufficiently large $m$. This implies that a steady shear flow of a second-order fluid is unstable to short-wavelength perturbations and cannot be realized. In the case of negligible inertia, achieved in the above by setting $\rho =0$, this predicts that all Fourier modes are unstable. The implications of this result are profound: it shows that an approximation of slow flows, or, in other words, Taylor expansions of the stress in terms of the relaxation time cannot be used in time-dependent flows where any shear component would result locally in a linear instability and exponential growth of the stress. Hence, the second-order fluid and similar approximations should generally not be used to study time-dependent flows!