1D Pipe Flow

Geodynamics

(Lecture 5)

%%javascript
MathJax.Hub.Config({
    TeX: { equationNumbers: { autoNumber: "AMS" } }
});
MathJax.Hub.Queue(
  ["resetEquationNumbers", MathJax.InputJax.TeX],
  ["PreProcess", MathJax.Hub],
  ["Reprocess", MathJax.Hub]
);

Channel Flow

We consider the one-dimensional flow of a Newtonian viscous fluid in a channel between parallel plates as a model for channel flow.

Figure 6.1: (a) The force balance on a layer of fluid in a channel with an applied pressure gradient. (b) A typical velocity profile.

For a Newtonian fluid with constant viscosity \(\mu\), the shear stress \(\tau\) at any location in the channel is given by

(139)\[\tau =\mu \frac{d u}{dy}\]

It is not hard to prove that

(140)\[\frac{d \tau}{d y} = \frac{dp}{dx} = -\frac{p_1-p_0}{l}\]

Combining the equations (139) and (310), and integrating the result we have the velocity \(u\) profile as a function of \(y\):

(141)\[u=\frac{1}{2\mu}\frac{dp}{dx}y^2+c_1y+c_2\]

To evaluate the constants, we must satisfy the boundary conditions. For example, for a no-slip boundary conditions as shown in figure 6.1, a viscous fluid in contact with a solid boundary must have the same velocity as the boundary, so that \(u=0\) at \(y=h\) and \(u=u_0\) at \(y =0\). Hence, the equation becomes

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

import numpy as np
import matplotlib.pyplot as plt
N  = 10000    # number of points to plot (ignore)
mu = 4*(10**19)    # viscosity
dpdx = 0    # when pressure difference = 0, known as Couette flow
h = 200000    # thickness of channel (m)
u0 = 1.5844   # boundary (max) velocity (m/s)

def main():
    y = np.linspace(0, h, N)
    
    u = (1/(2*mu)) * dpdx * (y*y - h*y) - (u0*y)/(h) + u0

    plt.figure()
    plt.plot(y, u)
    plt.xlabel('$y (m)$')
    plt.ylabel('$u (m/s)$')
    plt.title('Velocity Profile')
    
    plt.show()

if __name__ == '__main__':
    main()

Problem 6.1

Show that the mean velocity in the channel is given by

(143)\[\bar{u}=-\frac{h^2}{12\mu}\frac{dp}{dx}+\frac{u_0}{2}\]

This can be easily proved by integrating velocity equation (312) wrt. \(y\) from \(0\) to \(h\), and then divided by the thickness of channel \(h\)

import numpy as np
import matplotlib.pyplot as plt
from matplotlib.patches import Polygon


def velocity_u(y):
    return (1/(2*mu)) * dpdx * (y*y - h*y) - (u0*y)/(h) + u0


a, b = 0, h  # integral limits
y = np.linspace(a, b, N)
print(y)
u = velocity_u(y)

fig, ax = plt.subplots()
ax.plot(y, u, 'r', linewidth=2)
ax.set_ylim(bottom=0)

def mean_u(a,b,y):
    bin_size = (b-a)/N
    area = 0
    for iy in y:
        u = velocity_u(iy)
        area += u * bin_size
    return area/(b-a)

x_coordinates = [a, b]
y_coordinates = [mean_u(a,b,y), mean_u(a,b,y)]

ax.plot(x_coordinates, y_coordinates, linestyle='dashed',label='Mean Velocity u')
ax.legend()


# Make the shaded region
ix = np.linspace(a, b)
iy = velocity_u(ix)
verts = [(a, 0), *zip(ix, iy), (b, 0)]
poly = Polygon(verts, facecolor='0.9', edgecolor='0.5')
ax.add_patch(poly)

ax.text(0.5 * (a + b), 0.5, r"$\int_0^h u(y)\mathrm{d}y$",
        horizontalalignment='center', fontsize=20)

fig.text(0.9, 0.05, '$y$')
fig.text(0.1, 0.9, '$u$')

ax.spines['right'].set_visible(False)
ax.spines['top'].set_visible(False)
ax.xaxis.set_ticks_position('bottom')

ax.set_xticks((a, b))
ax.set_xticklabels(('$0$', '$h$'))
ax.set_yticks([])


plt.show()

[0.00000000e+00 2.00020002e+01 4.00040004e+01 ... 1.99959996e+05
 1.99979998e+05 2.00000000e+05]

Problem 6.2

Derive a general expression for the shear stress \(\tau\) at any location \(y\) in the channel.

This can be easily derived by using equation (139), in which the velocity \(u\) is expressed by equation (312). Hence the result is

(144)\[\tau(y)=\frac{1}{2}\frac{dp}{dx}(2y-h) - \mu\frac{u_0}{h}\]

What are the simplified forms of \(\tau\) for Couette flow and for the case \(u_0 = 0\)?

Since we know that the pressure difference \(\frac{dp}{dx}=0\) for Couette flow, then we can obtain

(145)\[\tau=-\mu\frac{u_0}{h}\]

Also, if the boundary velocity \(u_0=0\), then

(146)\[\tau=\frac{1}{2}\frac{dp}{dx}(2y-h)\]

Problem 6.3

Find the point in the channel \(y\) at which the velocity is a maximum.

At this point \(y\), the derivative of the velocity (equation (312)) must be zero and the second derivative of it must be negative. Using this relationship, we can rearrange the equation and obtain

(147)\[y=\frac{h}{2}+\frac{u_0\mu_0}{h(dp/dx)}\]

Problem 6.4

Figure 6.3: Unidirectional flow of a constant thickness layer of viscous fluid down an inclined plane.

Consider the steady, unidirectional flow of a viscous fluid down the upper face of an inclined plane. Assume that the flow occurs in a layer of constant thickness h, as shown in figure 6.3. Show that the velocity profile is given by

(148)\[u(y)=\frac{\rho g \sin \alpha}{2\mu}(h^2-y^2)\]

where \(y\) is the coordinate measured perpendicular to the inclined plane (\(y = h\) is the surface of the plane), \(\alpha\) is the inclination of the plane to the horizontal, and \(g\) is the acceleration of gravity.

Imagine we have a tiny rectangle in the liquid as shown in figure 6.3, there would be three forces acting on this body: the dragging force \(F_{\tau(y)}\) to the downhill direction caused by shear stress \(\tau(y)\), the gravitational force \(F_g \sin\alpha\) pulling it to the downhill direction and the resisting drag force acting uphill \(-F_{\tau(y+\delta y)}\) that is due to the shear stress \(\tau(y+\delta y)\):

(149)\[F_{\tau(y)}=l \cdot \tau(y)\]

(150)\[-F_{\tau(y+\delta y)}=-l \cdot \tau(y+\delta y)=-l \cdot (\tau(y)+\frac{d\tau}{dy} \cdot \delta y)\]

(151)\[F_g \sin \alpha = mg\sin\alpha =v\rho g \sin\alpha= (\delta y) l \rho g \sin\alpha\]

If we assume that the fluid flows in a constant velocity, that is, there is zero acceleration, then the net forces in the x-direction would be 0, which means the sum of the equations (149), (150) and (151) would be zero.

(152)\[\sum F_x = F_{\tau(y)} - F{\tau(y+\delta y)} + F_g\sin\alpha = 0\]

Rearrange equation (152), we can obtain

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

Applying the shear stress from equation (294) into equation (153), we have

(154)\[\mu \frac{d^2 u}{d y^2} = -\rho g \sin\alpha\]

Integrate equation (154) and we have

(155)\[\frac{du}{dy}=-\frac{\rho g \sin\alpha}{\mu}y+C_1\]

Recall equation (294), \(\frac{du}{dy}\) is related to shear stress \(\tau\). So we can apply the free surface condition (\(\tau=0\) at \(y=0\)), and obtain the constant \(C_1 = 0\).

Integrate equation (154) again, we have

(156)\[u(y) = -\frac{\rho g \sin\alpha}{2\mu}y^2+C_1 y+C_2\]

where \(C_1=0\) from the free surface condition.

Then, we apply no-slip condition (\(u=0\)) at \(y=h\), and obtain \(C_2\).

(157)\[0 = -\frac{\rho g \sin\alpha}{2\mu}h^2+C_2\]

Finally, we plug in \(C_1\) and \(C_2\) into equation (156), we can show that:

(158)\[u (y) = \frac{\rho g \sin\alpha}{2\mu}(h^2-y^2)\]

It is not hard to obtain the mean veolocity in the layer by integrating the equation (158) from \(0\) to \(h\) and then divide it by \(h\):

(159)\[ \bar{u} = \frac{\rho g h^2 \sin\alpha}{3\mu} \]

What is the thickness of the layer whose rate of flow down the incline is \(Q\) (per unit width in the direction perpendicular to the figure plane)?

(160)\[Q = \bar{u} h\]

Rearrange we get

(161)\[h = \left(\frac{3\mu Q}{\rho g \sin\alpha}\right)^{\frac{1}{3}}\]

Pipe Flow

With subsequent applications to flows in aquifers and volcanic conduits in mind, we next consider viscous flow through a circular pipe.

Figure 6.6: Poiseuille flow through a circular pipe.

If we assume the fluid flows in constant velocity, then the driving force and resisting force should be balanced like this

(162)\[(p_1-p_0)\pi r^2 = -\tau \cdot 2\pi r l\]

If we rearrange the equation, we have the shear stress

(163)\[-\tau = \frac{(p_1-p_0)\pi r^2}{2\pi r l} = \frac{r}{2}\frac{p_1-p_0}{l} = \frac{r}{2}\frac{dp}{dx}\]

where dp/dx is the pressure gradient along the pipe.

In the cylindrical geometry in Figure 6.6, the shear stress \(\tau\) is directly proportional to the radial gradient of the velocity \(u\)

(164)\[ \tau=\mu\frac{du}{dr} \]

Combining equations (163) and (164), we have

(165)\[\frac{du}{dr}=\frac{r}{2\mu}\frac{dp}{dx}\]

which can be integrated to give

(166)\[u(r)=-\frac{1}{4\mu}\frac{dp}{dx}(R^2-r^2)\]

We used the condition \(u = 0\) at \(r = R\). The velocity profile in the pipe is a parabaloid of revolution; it is known as Poiseuille flow.

It is not hard to get the following quantities:

(167)\[Q=\int^R_0 2\pi r u \,dr = - \frac{\pi R^4}{8\mu}\frac{dp}{dx} = \pi R^2 \bar{\mu}\]

where \(\mu\) is the mean velocity:

(168)\[\bar{\mu}=-\frac{R^2}{8\mu}\frac{dp}{dx}=\frac{1}{2}u_{max}\]

and \(u_{max}\) is the maximum velocity.

(169)\[u_{max}=-\frac{R^2}{4\mu}\frac{dp}{dx}\]

It is often convenient in fluid mechanics to work in terms of dimensionless variables. The relation between the mean velocity in the pipe and the pressure gradient (168) can be put into standard dimensionless form by introducing two quantities: a dimensionless pressure gradient or friction factor f and the Reynolds number Re. The friction factor is defined as

(170)\[f=\frac{-4R}{\rho\bar{\mu}^2}\frac{dp}{dx}\]

and the Reynolds number is given by

(171)\[Re = \frac{\rho \bar{u}D}{\mu}\]

where \(D=2R\) is the pipe diameter.

Critical \(Re\) is around \(2200\) for pipe flow, above which the flow is turbulent and below which laminar. Although this number depends on the geometry of the pipe. Theoretically, for laminar flow, the relationship between \(f\) and \(Re\) is

(172)\[f = \frac{64}{Re} \]

Empirically, for turbulent flow, the relationship between \(f\) and \(Re\) is

(173)\[f=0.3164Re^{-\frac{1}{4}}\]

Problem 6.7

Determine \(Re\) for asthenoshpheric type channel flow, assuming: \(\mu=4 \times 10^{19}\,Pa \cdot s, \rho=3200\,kg/m^3\), channel thickness is \(200\,km\) and the flow is driven solely by \(50\,mm/yr\) of motion of the overlying plate.

(174)\[Re= \frac{3200 \times \bar{u} \times 200000}{4 \times 10^{19}}\]

dpdx = 0
mu = 4* 10**19
h = 200000
u0 = 1.5844 * 10 ** (-9)
a,b=0,h
u_bar = mean_u(a,b,y)
print(u_bar)

7.922000000000002e-10

(175)\[Re= \frac{3200 \times 7.92 \times 10^{-10} \times 200000}{4 \times 10^{19}}= 1.27 \times 10^{-20}\]

This problem illustrates that the viscosity of mantle rock is so high that the Reynolds number is generally small.

Problem 6.8

Figure 6.9: A semicircular aquifer with a circular cross section (a toroid). A hydrostatic head b is available to drive the flow.

With the parameters provided below, what is the radius of the channel \(R\), the average velocity \(\bar{u}\) and is the flow laminar or turbulent?

b=50    # m, because the entrance of spring is 50 m above it
Q = 1.667*10**(-3)    # Spring outputs: 100 L/min i.e.1.667*10**(-3) m3/s 
R_prime = 1000    # m, because entrance of spring 2km from outlet

Now, we do not know if the flow is laminar or turbulent, but we can start with assuming the flow is laminar. Hence, the laminar flux is

(176)\[Q = - \frac{\pi R^4}{8\mu} \times \frac{dp}{ds} = \frac{\rho g b R^4}{8\mu R^{'}}\]

Rearrange the equation we have the radius of the channel \(R\) (under the laminar flow assumption)

(177)\[R = \sqrt[4]{\frac{8Q\mu R^{'}}{\rho g b}} = 0.0133\,m\]

The average velocity is therefore

(178)\[\bar{u} = \frac{Q}{A} = \frac{\rho g b R^4}{8\mu R^{'} \pi R^2} = 3.01\,m/s\]

Now, we can verify that if the flow is laminar or turbulent by calculating the \(Re\), which is \(1.49 \times 10^7\). This is obviously larger than the threshold 2200, so the flow is actually turbulent. Thus, we need to recalculate the channel radius and mean velocity under the assumption that the flow is turbulent using the equation (173).

Plugging in the \(f\) (equation (170)) and \(Re\) (equation (171)) we defined earlier into the relationship in (173), we have the mean velocity \(\bar{u}\) after rearranging

(179)\[\bar{u} = \frac{2^{9/7}R^{5/7}\frac{dp}{dy}^{4/7}}{0.3164^{4/7}\mu^{1/7}\rho^{3/7}} = 0.84\,m/s\]

Using the simple relationship in the equation (178), we can obtain the area, and thus the radius (\(R=0.025\,m\)) of the channel.

Problem 6.9

Figure 6.10

Determine the rate at which magma flows up a twodimensional channel of width \(d\) under the buoyant pressure gradient \(-(\rho_s- \rho_l)g\). Assume laminar flow.

Again, we can assume that the fluid is steady and solve the force balance equation

(180)\[\sum F_y = F_{\tau(x)} - F{\tau(x+\delta x)} + F_p = 0\]

(181)\[F_{\tau(x)}=\tau(x) \cdot l =\mu \frac{d u}{dx} \cdot l\]

(182)\[F_{\tau(x+\delta x)}=\tau(x+\delta x) \cdot l= \tau(x)+\frac{d \tau}{d x}\delta x \cdot l =\mu (\frac{d u}{dx}+\frac{d^2u}{dx^2}) \cdot l\]

(183)\[F_p = (\rho_s-\rho_l)g \cdot l \cdot \delta x \]

Put equations (181), (182) and (183) into (180), we can obtain

(184)\[\frac{d^2u}{dx^2} = \frac{(\rho_s-\rho_l)g}{\mu}\]

Integrating it twice we have

(185)\[u(x) = \frac{(\rho_s-\rho_l)g}{2\mu}x^2 + C_1 x + C_2\]

Because of the symmetry of the velocity profile, we can believe that \(C_1\) is zero, otherwise the equation (185) would not be symmetric. And plugging in the no-slip boundary condition where \(u=0\) at \(x=\pm\frac{d}{2}\), we have \(C_2 = - \frac{(\rho_s-\rho_l)g}{8\mu}d^2\).

As a result, the velocity profile for the channel is

(186)\[u(x) = \frac{(\rho_s-\rho_l)g}{2\mu}x^2 - \frac{(\rho_s-\rho_l)g}{8\mu}d^2\]

(187)\[\bar{u}=\frac{1}{d}\int^{d/2}_{-d/2} u(x) \,dx = \frac{(\rho_s-\rho_l)g}{12\mu}d^2\]

(188)\[Q = \bar{u} d = \frac{(\rho_s-\rho_l)g}{12\mu}d^3\]