Non Linear Viscosity 2D Stokes#

(Lecture 6)

Although fracture is important in shallow crustal rock at low temperatures and pressures, there are many circumstances in which rock behaves as a ductile material. In determining the transition from brittle to ductile behavior, pressure, temperature, and strain rate are important.

At temperatures that are a signiﬁcant fraction of the melt temperature the atoms and dislocations in a crystalline solid become suﬃciently mobile to result in creep when the solid is subjected to deviatoric stresses. At very low stresses diﬀusion processes dominate, and the crystalline solid behaves as a Newtonian ﬂuid with a viscosity that depends exponentially on pressure and the inverse absolute temperature.

At higher stresses the motion of dislocations becomes the dominant creep process resulting in a non-Newtonian or nonlinear ﬂuid behavior that also has an exponential pressure and inverse absolute temperature dependence.

Figure 7.1: Shear strain in the fluid is measured by the displacement $dx$ between two nearby layers with a distance of $dy$.

(189)#\[ \epsilon = \frac{dx}{dy}\]

The rate of strain measures the speed of deformation, and is defined by

(190)#\[ \dot{\epsilon} = \frac{dx}{dy} \div t = \frac{du}{dy}\]

Recall that the cause for strain or deformation is stress, which is defined by

(191)#\[ \tau = \sigma = \mu \frac{du}{dy}\]

where $\tau$ and $\sigma$ are used interchangeably for stress, and $\mu$ (or sometimes $\eta$) is viscosity.

Therefore, we have obtained the relationship between stress and strain rate:

(192)#\[\tau = \mu \dot{\epsilon}\]

Figure 7.2: the gradient of the figure on the left hand side represents the inversity of viscosity. Dash line on the figure at the right represents the effective viscosity that changes with stress.

(193)#\[ \dot{\epsilon}=A\sigma^n\]

where $A = \frac{1}{\mu}$.

If the strain rate is linearly proportional to the change of stress ($n = 1$), then the fluid is Newtonian (diffusion creep).

If $n>1$, the fluid is non-Newtonian (dislocation creep).

If $n\rightarrow\infty$, the fluid is pseudo-plastic.

Effective Viscosity#

The viscosity of drilling mud will change with a change in the shear rate because in some degree drilling fluids are shear-thinning. Because of its nature of drilling mud, people create the new term of viscosity called “Effective Viscosity” to compensate the change in shear rate of viscosity. By definition, the effective viscosity means the viscosity of Newtonian fluid that gives the same shear stress at the same shear rate.

As we can see from figure 7.2, the viscosity changes with stress in a non-Newtonian fluid. However, based on lab experiments, the variable viscosity (not constant anymore), or effective viscosity, is dependent on not only stress $\sigma$ but also many factors (e.g. temperature $T$):

(194)#\[ \mu = C \left(\frac{\sigma}{G}\right)^{1-n}\exp\left[\frac{E^{*}+pV^{*}}{RT}\right]\]

where the details of the equation can be found online and

(195)#\[C=\frac{1}{2}\left(\frac{G}{A}\frac{d}{b}\right)^{m}\]

In this chapter, we mainly discuss the channel flow of stress-dependent viscosity and temperature-dependent viscosity.

Problem 7.10 (Stress-dependent viscosity)#

Suppose we assume the viscosity of the fluid is only stress-dependent, then the equation (194) can be simplified to a power-law (equation (193) ) with respect to stress $\sigma$. Then, we can get the analytical solution of effective viscosity and other quantities for the fluid in channel flow. Suppose we have no-slip boundary at $y=\pm h/2$ and the strain rate obeys the power law.

Figure 7.3: The velocity profile is symmetric due to the boundary conditions and the geometry of the channel.

Balancing viscous shear and pressure forces, we have

(196)#\[\frac{d\tau}{dy}=\frac{dp}{dx}\]

Therefore, we can obtain the equation for stress by integrating wrt. y

(197)#\[\tau(y)=\frac{dp}{dx}y+C\]

Since we know that the stress can be expressed as

(198)#\[\tau = \mu \frac{du}{dy}\]

and that the $\frac{du}{dy}=0$ at $y=0$ by inspecting the velocity profile in figure 7.3, we can obtain

(199)#\[\tau (y=0) = 0\]

So we know that the integral constant $C$ in (197) is zero.

Also, we have an non-linear relationship between the stress and strain rate for the non-Newtonian fluid as follow

(200)#\[\frac{du}{dy}=C_1\tau^{n}\]

We can plug in the $\tau$ from (200) into (196) and arrive at

(201)#\[\frac{du}{dy}=C_1\left(\frac{dp}{dx}\right)^{n}y^n\]

Further integrating (201), and using the boundary condition $u(y=\pm h/2)=0$, we have the velocity profile

(202)#\[u=-\frac{C_1}{n+1}\left(\frac{dp}{dx}\right)^n \left[\left(\frac{h}{2}\right)^{n+1}-y^{n+1}\right]\]

We can also obtain the mean velocity across the channel

(203)#\[\bar{u}=\frac{C_1}{n+2}\left(\frac{dp}{dx}\right)^n\left(\frac{h}{2}\right)^{n+1}\]

Now, using the idea that the effective viscosity is the ratio between stress and strain rate, we have

(204)#\[\mu_{eff}=\frac{\tau}{\frac{du}{dy}}=\frac{\frac{dp}{dx}y}{C_1\left(\frac{dp}{dx}\right)^{n}y^n} = \frac{dp}{dx}\frac{h^2}{4(n+2)\bar{u}}\left(\frac{2y}{h}\right)^{1-n}\]

where $\frac{dp}{dx} = \frac{p_1-p_0}{L}$

And using the boundary condition, we can derive the scaled effective viscosity (as a function of $y/h$ or scaled $y$ distance where $-h/2 < y < h/2$).

(205)#\[\frac{\mu_{eff}}{\mu_{eff,wall}}=\left(\frac{2y}{h}\right)^{1-n}\]

We can plot (205) for $n=1,3,5$ and for $-0.5<y/h<0.5$.

import numpy as np
import matplotlib.pyplot as plt
N  = 10000    # number of points to plot (ignore)

n_list = [1,3,5]

def main():
    y_h = np.linspace(-0.5, 0.5, N)
    
   

    for n in n_list:
        scal_eff_mu = (2*y_h)**(1-n)
        plt.figure()
        plt.plot(y_h, scal_eff_mu)
        plt.legend(str(n))
        plt.xlabel('$y/h$')
        plt.ylabel('$scaled effective viscosity$')
        plt.title('Scaled Effective Viscosity Profile')
        plt.ylim(0, 5)
        plt.show()

if __name__ == '__main__':
    main()

../../_images/Lecture6_Non_Linear_Viscosity_2D_Stokes_50_0.png

../../_images/Lecture6_Non_Linear_Viscosity_2D_Stokes_50_1.png

../../_images/Lecture6_Non_Linear_Viscosity_2D_Stokes_50_2.png

Problem 7.12 (T-dependent viscosity)#

Now, let’s assume that the viscosity of the fluid is only temperature-dependent. We can determine the velocity profile $u(y)$ for a classic model as shown below:

Figure 7.4: Icesheet with thickness $h$ on sloping bedrock, creeping downhill under its own weight. The temperature profile is linear and $(T_1-T_0)/T_0$ is small. It is no-slip at base.

First, we can write out the temperature profile

(206)#\[T(y) = T_0+(T_1-T_0)\frac{y}{h}\]

We assume the temperature dependence of the pre-exponential factor as negligible compared with the temperature dependence of the exponential factor, so the equation (194) simplifies to

(207)#\[\mu=Ce^{E^{*}/RT}\]

Plugging the equation (206) in (207), and using Taylor expansion for the exponential term, we have

$$\mu=C \exp{\left[\frac{E^{*}}{RT_0}\left(1-\frac{T_1-T_0}{T_0}\frac{y}{h}\right)\right]}$$ (equation6_19)}

Also, from the forces analysis shown by figure 6.3 in lecture 5, we have

(208)#\[\frac{d\tau}{dy}= \rho g \sin\alpha\]

Integrating it with the boundary condition $\tau = 0$ at $y=0$, we have

(209)#\[\tau(y)= \rho gy \sin\alpha\]

Since effective viscosity is the ratio between stress $\tau$ and strain rate $du/dy$, we can combine equations equation6_19 and (209) and have

(210)#\[\frac{du}{dy}=\frac{\rho gy \sin\alpha}{Ce^{\frac{E^{*}}{RT_0}}e^{-\frac{E^{*}(T_1-T_0)y}{RT_0^2h}}}=G\frac{y}{e^{-By}}=Gye^{By}\]

where we defined constant $G$ and $B$ for simplicity.

After integration, we have

(211)#\[u(y)=-\frac{Ge^{By}}{B^2}+\frac{Gye^{By}}{B}+C_1\]

Since it is no-slip base, we have the boundary condition $u=0$ at $y=h$. Apply this boundary condition we can solve the constant

(212)#\[C_1=\frac{Ge^{Bh}}{B^2}-\frac{Ghe^{Bh}}{B}\]

Hence, the equation of the velocity profile can be obtained

(213)#\[u(y)=-\frac{Ge^{By}}{B^2}+\frac{Gye^{By}}{B}+\frac{Ge^{Bh}}{B^2}-\frac{Ghe^{Bh}}{B}\]

import numpy as np
import matplotlib.pyplot as plt
N  = 10000    # number of points to plot (ignore)

h = 1    # ice sheet thickness
alpha = 20    # slope 
rho = 1000    # density of ice
g = 10    # gravitational acceleration
C = 1000000  # term before expential in equation 6, tbh i didnt check the value for ice, just a random value here
E = 75    # creep activation energies of ice
R = 8     # universal gas constant
T0 = 9    # surface temperature 
T1 = 10   # base temperature hotter because heat absorption and release

G = np.sin(alpha)*rho*g/(C*np.exp(E/(R*T0)))
B = E*(T1-T0)/(R*T0^2*h)
C1 = G*np.exp(B*h)/(B*B) - G*h*np.exp(B*h)/B


def main():
    y = np.linspace(0, h, N)
    u = G*np.exp(B*y)/(B*B) - G*y*np.exp(B*y)/B + C1
    plt.figure()
    plt.plot(y, u)
    
    plt.xlabel('$y(m)$')
    plt.ylabel('$u(y)$')
    plt.title('Velocity Profile for Ice-sheet on a Slope')
    #plt.ylim(-10000, 10)
    plt.show()
    
if __name__ == '__main__':
    main()

2-D Stokes#

Problem 6.10#

Determine the stream function for the general one-dimensional channel flow in lecture 5.

(214)#\[u(y)=\frac{1}{2\mu}\frac{dp}{dx}(y^2-hy)-\frac{u_0y}{h}+u_0=-y+5\]

where we let $\frac{1}{2\mu}\frac{dp}{dx}=0$, $h=5$ and $u_0=5$ for simplicity.

$y$	$u$
0	5
1	4
2	3
3	2
4	1
5	0

Figure 7.5: One-dimensional Channel Flow.

The stream function is a function in connnection with a fluid flow denoted by $\psi$ who partial differentiation with $y$ gives the velocity in $x$ direction (i.e. $u$), and negative partial differentiation with $x$ gives the velocity of fluid in $y$ direction (i.e. $v$).

This is very confusing as in how the differentiation wiht $y$ gives the velocity in $x$ direction? It would certainly make more sense if it gives the velocity in $y$ direction…Do not worry, let’s try to understand the stream function by looking at it in a one-dimensional channel flow before we move in to the two-dimensional case.

But let’s take a look at the definition of stream function $\psi(x,y)$ in a two-demensional case

(215)#\[u=\frac{d\psi}{dy}\]

(216)#\[v=-\frac{d\psi}{dx}\]

For the one-dimensional channel flow problem I have shown above, there is no velocity in the vertical direction, and so we only consider the velocity in the x-direction.

(217)#\[u=\frac{d\psi}{dy}=-y+5\]

If we integrate it wrt. $y$, then we can obtain the equation for the stream function

(218)#\[\psi=-\frac{1}{2}y^2+5y+C(x)\]

Discuss the physical meaning of $\psi$ for this case.

So what does that $\psi$ function mean?

Figure 7.6: Volumetric Flow Rate of One-dimensional Channel Flow (with the width of z-direction equals 1).

Let’s look at the tiny volumetric flow rate $dQ$ through the tiny slice of thickness $dy$ as shown above.

(219)#\[dQ=u\,dy \cdot 1 = u\,dy\]

Then the overall volumetric flow rate through the whole section (from $y=1$ to $y=2$) as shown in between the blue lines is

(220)#\[Q=\int^2_1 u\,dy=\int^2_1 (-y+5)\,dy=-\frac{1}{2}y^2+5y+C(x)|^2_1\]

Compare the final part of equation (220) with (218), we realise that the stream function represents volumetric flow rate in a sense

(221)#\[Q = \psi_2 - \psi_1\]

Now, we should finally understand the very confusing equation (215) we mentioned before, with the help of figure 7.6. The horizontal velocity $u$ equals the volumetric flow rate through the thin section $dQ$ divided by the vertical thickness $dy$, in this case

(222)#\[u=\frac{d\psi}{dy}=\frac{dQ}{dy} = \frac{d}{dy}\left(-\frac{1}{2}y^2+5y+C(x)\right)= -y+5\]

Sketch stream function and streamlines for the two end-member cases (Couette and P-driven). Pay attention to stream line orientation and spacing (in terms of colorbars).

So the streamlines represent the path of the particle in the channel flow, as we can imagine the path in this one-dimensional flow problem should be linear and horizontal. In another word, the particle should only flow in the x-direction as in there is no force acting on the y-direction.

from mpl_toolkits import mplot3d
import numpy as np
import matplotlib.pyplot as plt
N  = 100    # number of points to plot (ignore)
mu = 0.001    # viscosity of water
dpdx = 0    # when pressure difference = 0, known as Couette ﬂow
h = 5    # thickness of channel (m)
u0 = 5   # boundary (max) velocity (m/s)
x_range = 10    # length of the channel

def psi(x, y):
    return 1/(6*mu)*dpdx*y**3 - h/(4*mu)*dpdx*y**2 - u0/(2*h)*y**2 + u0*y


x = np.linspace(0, h, N)
y = np.linspace(0, h, N)

X, Y = np.meshgrid(x, y)
Z = psi(X, Y)

fig = plt.figure(figsize=(7, 5))
ax = plt.axes(projection='3d')
ax.contour3D(X, Y, Z, 50, cmap='binary')
ax.set_xlabel('x')
ax.set_ylabel('y')
ax.set_zlabel('$\psi$')
ax.set_title('Stream Function $\psi$ in 3D Contours (Couette)')
plt.show()

../../_images/Lecture6_Non_Linear_Viscosity_2D_Stokes_103_0.png

import numpy as np
import matplotlib.pyplot as plt
import matplotlib.gridspec as gridspec

N  = 100    # number of points to plot (ignore)
mu = 0.001    # viscosity of water
dpdx = 0    # when pressure difference = 0, known as Couette ﬂow
h = 5    # thickness of channel (m)
u0 = 5   # boundary (max) velocity (m/s)


Y, X = np.mgrid[0:h:100j, 0:h:100j]
U = (1/(2*mu)) * dpdx * (Y*Y - h*Y) - (u0*Y)/(h) + u0
V = Y*0
speed = np.sqrt(U**2 + V**2)

fig = plt.figure(figsize=(14, 18))
gs = gridspec.GridSpec(nrows=3, ncols=2, height_ratios=[1, 1, 2])


# Varying color along a streamline
ax1 = fig.add_subplot(gs[0, 1])
strm = ax1.streamplot(X, Y, U, V, color=U, linewidth=2, cmap='autumn')
fig.colorbar(strm.lines)
ax1.set_title('Streamlines (Couette)')
ax1.set_xlabel('x')
ax1.set_ylabel('y')


plt.tight_layout()
plt.show()

../../_images/Lecture6_Non_Linear_Viscosity_2D_Stokes_104_0.png

from mpl_toolkits import mplot3d
import numpy as np
import matplotlib.pyplot as plt
N  = 100    # number of points to plot (ignore)
mu = 0.001    # viscosity of water
dpdx = -10    # when pressure difference = 0, known as Couette ﬂow
h = 5    # thickness of channel (m)
u0 = 5   # boundary (max) velocity (m/s)
x_range = 10    # length of the channel

def psi(x, y):
    return 1/(6*mu)*dpdx*y**3 - h/(4*mu)*dpdx*y**2 - u0/(2*h)*y**2 + u0*y


x = np.linspace(0, h, N)
y = np.linspace(0, h, N)

X, Y = np.meshgrid(x, y)
Z = psi(X, Y)

fig = plt.figure(figsize=(7, 5))
ax = plt.axes(projection='3d')
ax.contour3D(X, Y, Z, 50, cmap='binary')
ax.set_xlabel('x')
ax.set_ylabel('y')
ax.set_zlabel('$\psi$')
ax.set_title('Stream Function $\psi$ in 3D Contours (P-driven)')
plt.show()

../../_images/Lecture6_Non_Linear_Viscosity_2D_Stokes_105_0.png

import numpy as np
import matplotlib.pyplot as plt
import matplotlib.gridspec as gridspec

N  = 100    # number of points to plot (ignore)
mu = 0.001    # viscosity of water
dpdx = -10    # pressure driven flow
h = 5    # thickness of channel (m)
u0 = 5   # boundary (max) velocity (m/s)


Y, X = np.mgrid[0:h:100j, 0:h:100j]
U = (1/(2*mu)) * dpdx * (Y*Y - h*Y) - (u0*Y)/(h) + u0
V = Y*0
speed = np.sqrt(U**2 + V**2)

fig = plt.figure(figsize=(14, 18))
gs = gridspec.GridSpec(nrows=3, ncols=2, height_ratios=[1, 1, 2])


# Varying color along a streamline
ax1 = fig.add_subplot(gs[0, 1])
strm = ax1.streamplot(X, Y, U, V, color=U, linewidth=2, cmap='autumn')
fig.colorbar(strm.lines)
ax1.set_title('Streamlines (P-driven)')
ax1.set_xlabel('x')
ax1.set_ylabel('y')


plt.tight_layout()
plt.show()

../../_images/Lecture6_Non_Linear_Viscosity_2D_Stokes_106_0.png

ESE Jupyter Material

Non Linear Viscosity 2D Stokes

Contents