Poiseuille flow

Introduction

In this notebook we will be studying a incompressible fluid experiencing steady Poiseuille flow. Poiseuille flow is named after Jean Poiseuille (left) a French physicist who derived the equation for how fluids in a pipe move, this equation was also derived by Gotthilf Hagen (middle) at the same time and later justified by George Stokes (right) in 1845.

Poiseuille flow occurs when a fluid driven by a constant pressure gradient flows between two non-moving plates and is both axisymmetric and swirl free. Blood and pipe flow are examples of Poiseuille flow.

Sources: Jean Poiseuille, Gotthilf Hagen, George Stokes

Code

# import necessary libraries 
import plotly.graph_objects as go
import plotly.figure_factory as ff
import plotly.express as px
import numpy as np
import pandas as pd
import plotly.io as pio
pio.renderers.default = 'iframe'

Defining the Velocity

It will be very useful to consider a free body diagram of a fluid blob.

Code

# plot force diagram
fig = go.Figure()
fig.add_shape(type='rect', x0=0, y0=1, x1=4, y1=3, fillcolor='rgba(133, 163, 201, 0.5)')
fig.add_scatter(x=np.linspace(-1,5, 20), y=np.ones(20)*3.5, line_color='black')
fig.add_scatter(x=np.linspace(-1,5, 20), y=np.ones(20)*0.5, line_color='black')
fig.add_annotation(x=-0.1, y=2, ax=-1, ay=2, xref='x', yref='y', axref='x', ayref='y',arrowhead=5, text=
                   r'$\require{color}{\color[rgb]{0.315209,0.728565,0.037706}p}$')
fig.add_annotation(x=4.1, y=2, ax=6, ay=2, xref='x', yref='y', axref='x', ayref='y', arrowhead=5, text=
                   r'$\require{color} {\color[rgb]{0.315209,0.728565,0.037706}p} + \frac{\partial \
                   {\color[rgb]{0.315209,0.728565,0.037706}p}}{\partial x} \delta x$')
fig.add_annotation(x=1, y=0.9, ax=3, ay=0.9, xref='x', yref='y', axref='x', ayref='y', arrowhead=5, 
                   text=r'$\require{color}{\color[rgb]{0.990308,0.800015,0.121270}\tau}$')
fig.add_annotation(x=3.5, y=3.25, ax=0.5, ay=3.25, xref='x', yref='y', axref='x', ayref='y', arrowhead=5, text=
                  r'$\require{color}{\color[rgb]{0.990308,0.800015,0.121270}\tau} + \frac{\partial \
                  {\color[rgb]{0.990308,0.800015,0.121270}\tau}}{\partial y} \delta y$')
fig.add_scatter(x=[2], y=[3.65], mode='text', text=r'$\require{color}{\color[rgb]{0.059472,0.501943,0.998465}v}=0$')
fig.add_scatter(x=[2], y=[0.4], mode='text', text=r'$\require{color}{\color[rgb]{0.059472,0.501943,0.998465}v}=0$')
fig.add_annotation(x=-1.2, y=3.5, ax=-1.2, ay=0.4, xref='x', yref='y', axref='x', ayref='y', arrowhead=5)
fig.add_scatter(x=[-1.3], y=[2], mode='text', text=r'$H$')
fig.update_layout(yaxis=dict(scaleanchor="x", scaleratio=1, showticklabels=False), xaxis=dict(showticklabels=False),
                  showlegend=False, plot_bgcolor='rgba(0, 0, 0, 0)')
fig.show()

To find the acceleration of a fluid particle we need to use the material derivative:

\[ \large \require{color} \vec{a} = \frac{D{\color[rgb]{0.059472,0.501943,0.998465}v}}{Dt} = \frac{\partial {\color[rgb]{0.059472,0.501943,0.998465}v}} {\partial t} + \left({\color[rgb]{0.059472,0.501943,0.998465}v} \cdot \nabla \right){\color[rgb]{0.059472,0.501943,0.998465}v} \]

Due to the fact that the flow is steady

\[ \large \frac{\partial {\color[rgb]{0.059472,0.501943,0.998465}v}}{\partial t} = 0. \]

Additionally, the velocity in the x-direction is uniform. This implies

\[ \large \frac{\partial {\color[rgb]{0.059472,0.501943,0.998465}v}}{\partial x}=0 \]

and

\[ \large {\color[rgb]{0.059472,0.501943,0.998465}v}_y=0, \]

therefore

\[ \left({\color[rgb]{0.059472,0.501943,0.998465}v} \cdot \nabla \right){\color[rgb]{0.059472,0.501943,0.998465}v} = \left(v_x \frac{\partial}{\partial x} + v_y \frac{\partial}{\partial y}\right)v = 0. \]

Thus, acceleration in both the x and y directions is zero which tells us that \(F_{x}=0\) and \(F_{y}=0\). There are two forces acting on the flow, shear and pressure, therefore the vectoral sum of these two forces must sum to zero. The shear stress on the fluid is a function only of y and the pressure is a function only of x, therefore in the y-direction

\[ \large \require{color}\frac{\partial {\color[rgb]{0.990308,0.800015,0.121270}\tau}}{\partial x} = \frac{\partial {\color[rgb]{0.315209,0.728565,0.037706}p}}{\partial y} = 0 \]

For the x components of shear and pressure to sum to zero they must both be constant along their respective axes. By using the above force diagram we can get a relationship between \(\require{color}\frac{\partial {\color[rgb]{0.315209,0.728565,0.037706}p}}{\partial x}\) and $ $.

By summing the forces, we write \[ \large \require{color}F_{x} = 0 \] \[ \large = {\color[rgb]{0.315209,0.728565,0.037706}p}\delta y - \left({\color[rgb]{0.315209,0.728565,0.037706}p}+\frac{\partial {\color[rgb]{0.315209,0.728565,0.037706}p}}{\partial x} \delta x\right)\delta y - {\color[rgb]{0.990308,0.800015,0.121270}\tau}\delta x + \left({\color[rgb]{0.990308,0.800015,0.121270}\tau} + \frac{\partial {\color[rgb]{0.990308,0.800015,0.121270}\tau}}{\partial y} \delta y\right)\delta x \] \[ \large \Rightarrow \frac{\partial {\color[rgb]{0.315209,0.728565,0.037706}p}}{\partial x} = \frac{\partial {\color[rgb]{0.990308,0.800015,0.121270}\tau}}{\partial y} \]

As the pressure is only a function of x and shear stress is only a function of y we can say the partial derivatives from above can be converted to total derivatives:

\[ \large \require{color}\frac{d{\color[rgb]{0.315209,0.728565,0.037706}p}}{dx} = \frac{d{\color[rgb]{0.990308,0.800015,0.121270}\tau}}{dy}. \]

Furthermore, we know the equation for shear stress

\[ \large \require{color}{\color[rgb]{0.990308,0.800015,0.121270}\tau} = {\color[rgb]{0.990448,0.502245,0.032881}\mu} \frac{d{\color[rgb]{0.059472,0.501943,0.998465}v}_{x}}{dy} \]

so we can get a relationship between velocity and the pressure gradient

\[ \large \require{color}\frac{d}{dy}\left({\color[rgb]{0.990448,0.502245,0.032881}\mu} \frac{d{\color[rgb]{0.059472,0.501943,0.998465}v}_{x}}{dy}\right) = {\color[rgb]{0.990448,0.502245,0.032881}\mu} \frac{d^{2}{\color[rgb]{0.059472,0.501943,0.998465}v}_{x}}{dy^{2}} = \frac{d{\color[rgb]{0.315209,0.728565,0.037706}p}}{dx}. \]

In Poiseuille flow the pressure gradient is constant but not zero, therefore we can solve the above equation for \(\require{color}{\color[rgb]{0.059472,0.501943,0.998465}v}_{x}\) by assuming \(\require{color}\frac{\partial {\color[rgb]{0.315209,0.728565,0.037706}p}}{\partial x}\) is some constant value thus the solution is in the form \(\require{color}{\color[rgb]{0.059472,0.501943,0.998465}v}_{x} = \frac{1}{2{\color[rgb]{0.990448,0.502245,0.032881}\mu}}\frac{d{\color[rgb]{0.315209,0.728565,0.037706}p}}{dx}y^{2} + By + C\) where B and C are determined by the boundary conditions.

The boundary conditions are: at \(\require{color}y=0\ {\color[rgb]{0.059472,0.501943,0.998465}v}=0\) and at \(\require{color}y=H\ {\color[rgb]{0.059472,0.501943,0.998465}v}=0\) from these we get that

\[ \large \require{color}{\color[rgb]{0.059472,0.501943,0.998465}v}_{x} = \frac{1}{2{\color[rgb]{0.990448,0.502245,0.032881}\mu}}\frac{d{\color[rgb]{0.315209,0.728565,0.037706}p}}{dx}y(y - H). \]

Given the above, let us now define the necessary variables.

Code

# define necessary variables 
H=10
dpdx = -1e-2
mu = 1e-3
x=np.zeros(11)
y=np.linspace(0, H, 11)
v_x= (1/(2*mu))*(dpdx)*y*(y-H)
v_y=np.zeros(11)

Visualizing the Flow

Visualizing How the Flow Moves

The code below is used for visualizing how Poiseuille flow behaves.

Code

# create a vector plot showing the velocity vector of particles along y direction each with same x position 
fig = ff.create_quiver(x=x, y=y, u=v_x, v=v_y, arrow_scale=0.04)
fig.update_layout(xaxis=dict(range=[-1, 20]))
fig.add_scatter(x=np.linspace(-1,15, 20), y=np.ones(20)*10, line_color='black')
fig.add_scatter(x=np.linspace(-1,15, 20), y=np.zeros(20), line_color='black')
fig.update_layout(yaxis=dict(scaleanchor="x", scaleratio=1, showticklabels=False), xaxis=dict(showticklabels=False),
                  showlegend=False, plot_bgcolor='rgba(0, 0, 0, 0)')
fig.show()

The plot above shows the velocity vectors at 9 points between the two non-moving walls and illustrates how fast the particles in each “layer” of fluid are moving

Code

# define how long it takes for the fastest particle to travel the length of the pipe
max_v = v_x[list(v_x).index(max(v_x))]
time = 20/max_v
distance = time*v_x

Code

# using this time animate how the particles would idealy move in Poiseuille flow
fig = go.Figure()
frames = [go.Frame(data=[]) for k in range(50)]
for i, y_pos in enumerate(y[1:len(y)-1]):
    fig.add_scatter(x=np.linspace(0, 20, 20), y=np.ones(20)*y_pos, line_dash='dash', line_color='lightsteelblue')
    x_pos = np.linspace(0, distance[i+1], 50)
    for j,f in enumerate(frames):
        f.data += (go.Scatter(x=[x_pos[j]], y=[(np.ones(50)[j]*y_pos)], mode='markers', marker=dict(color='teal', size=10)),)

fig.update(frames=frames)

fig.update_layout(updatemenus = [dict(type='buttons', buttons=[dict(label='play', method='animate', 
                  args=[None, {'frame':{'duration':50}}])])], showlegend=False, xaxis=dict(range=[0,20]))

fig.add_scatter(x=np.linspace(0, 20, 20), y=np.ones(20)*H, line_color='black')
fig.add_scatter(x=np.linspace(0, 20, 20), y=np.zeros(20), line_color='black')
fig.update_layout(yaxis=dict(scaleanchor="x", scaleratio=1, showticklabels=False), xaxis=dict(showticklabels=False),
                  showlegend=False, plot_bgcolor='rgba(0, 0, 0, 0)')
fig.show()

The animation above shows how an ideallized fluid would move through a pipe or non-moving boundary and demonstrates how the flow follows a parabolic velocity profile dependent on y (vertical) coordinate.

Adverse and Favorable Pressure Gradients

All of the figures in the above section, 3a, are based on the assumption that there is a negative pressure gradient. A negative pressure gradient is considered to be favorable because it allows the flow to move in the direction of initial velocity. However, when the pressure gradient is greater than zero it is considered to be adverse. If the flow were Couette+Poiseuille and the pressure gradient were adverse the flow near the bottom unoving plate would be pushed oppoiste the direction of the moving top plate. In just Poiseuille flow the the velocity profile reverses and the flow would move in the opposite direction to the initial velocity, which is in this case rightwards

Code

# define an array of pressure gradients
dpdxs = np.linspace(-1e-2, 1e-2, 17)

# create a dictionary of all the y position and velocity data at each pressure gradient 
all_data = {}
for p in dpdxs:
    all_data[p] = {'y':[], 'v':[]}
    for y_val in y:
        v_xy = ((1/(2*mu))*(p)*y_val*(y_val-H))
        all_data[p]['y'].append(y_val)
        all_data[p]['v'].append(v_xy)
        
# create a plot showing how velocity profile changes as the pressure gradient changes 
fig = go.Figure()
fig.add_scatter(x=np.linspace(-13, 13, 20), y=np.ones(20)*H, line_color='black')
fig.add_scatter(x=np.linspace(-13, 13, 20), y=np.zeros(20), line_color='black')

frames = []
for k in range(0,len(dpdxs)):
    figaux = ff.create_quiver(x=x, y=all_data[dpdxs[k]]['y'], u=all_data[dpdxs[k]]['v'], v = v_y, arrow_scale=0.04, 
                              line_color='steelblue')
    frames.append(go.Frame(data=[figaux.data[0]]))

fig.update_layout(updatemenus=[dict(type='buttons', buttons=[dict(label='Play', method='animate',
                  args=[None, dict(frame=dict(duration=500, redraw=True))])])], xaxis=dict(range=[-13,13], showticklabels=
                  False), yaxis=dict(showticklabels=False), showlegend=False, plot_bgcolor='rgba(0, 0, 0, 0)')
fig.update(frames=frames)

for k in range(len(fig.frames)):
    fig.frames[k]['layout'].update(title_text=f'Pressure Gradient {round(dpdxs[k],3)}')

fig.add_scatter(x=np.linspace(-13, 13, 20), y=np.ones(20)*H, line_color='black')
fig.add_annotation(x=5, y=-1, ax=-5, ay=-1, xref='x', yref='y', axref='x', ayref='y', arrowhead=5, text = r'$+x$')

fig.show()

In the above animation you can see that when

\[ \large \require{color}\frac{d{\color[rgb]{0.315209,0.728565,0.037706}p}}{dx} > 0 \]

the flow is moving to the left opposing the initial rightwards velocity of the flow. When the flow instead has

\[ \large \require{color}\frac{d{\color[rgb]{0.315209,0.728565,0.037706}p}}{dx} < 0 \]

the flow moves rightwards in the same direction as the initial velocity.