Code
# import necessary libraries
import plotly.graph_objects as go
import numpy as np
import plotly.io as pio
pio.renderers.default = 'iframe'Flow around a bend
Introduction
In this notebook we will study how a fluid acts when it flows around a bend opposed to moving through a straight section. As the fluid is moving around a bend it will be subject to centripetal force, a force that acts within a rotating reference frame and keeps the particles moving in circular motion. Therefore, we will need to use Newton’s second law of motion
\[ \large F={\color[rgb]{0.501961,0.250953,0.010028}m}a \]
to see how the pressure gradient affects the velocity throughout the bend. Newton published this law along with his others in his 1687 book the Principia and is pictured here:
To understand the relationship between pressure and velocity we need Bernoulli’s equation which was qualitatively postulated by Daniel Bernoulli in 1738 and formally written by Leonhard Euler in 1752 as
\[ \large \require{color}{\color[rgb]{0.315209,0.728565,0.037706}p}_{0}={\color[rgb]{0.315209,0.728565,0.037706}p}+\frac{1}{2}{\color[rgb]{0.814433,0.253157,0.091125}\rho} {\color[rgb]{0.059472,0.501943,0.998465}v}^{2} \]
where \({\color[rgb]{0.315209,0.728565,0.037706}p}_0\) is the total or stagnation pressure, \({\color[rgb]{0.315209,0.728565,0.037706}p}\) is the static pressure, and \(\frac{1}{2}{\color[rgb]{0.814433,0.253157,0.091125}\rho} {\color[rgb]{0.059472,0.501943,0.998465}v}^{2}\) is the dynamic pressure. Bernoulli (left) and Euler (right) are pictured here:
 |
 |
Sources: Isaac Newton, Daniel Bernoulli, Leonhard Euler
The velocity of a flow moving through a bend can be represented in polar coordinates: \[ \large \require{color}\vec{{\color[rgb]{0.059472,0.501943,0.998465}v}} = {\color[rgb]{0.059472,0.501943,0.998465}v}_{r}\hat{e}_{r} + {\color[rgb]{0.059472,0.501943,0.998465}v}_{\theta}\hat{e}_{\theta}. \] As the flow is steady the streamlines will follow a circular path, therefore the radial component of velocity is zero, and the velocity will point in the angular direction
Defining the Velocity and Pressure
It will be useful to define the velocity and pressure and see how they are related. As we will be working in polar coordinates, we initialize \(\theta\) values for the plots below.
Code
# initialize theta values for graph
t1=np.linspace(-np.pi/8,np.pi/8,50)
t2=np.linspace(-np.pi/2, np.pi/2,50)Now, we will plot a free body diagram for a small elemental area of the bend.
Code
# plot duct and force diagram for small area of flow
fig = go.Figure()
fig.add_scatter(x=2*np.cos(t2), y=2*np.sin(t2), fill='tozeroy', fillcolor='rgba(133, 163, 201, 0.5)', line_color='black')
fig.add_scatter(x=np.cos(t2), y=np.sin(t2), fill='tozeroy', fillcolor='white', line_color='black')
fig.add_scatter(x=1.25*np.cos(t1), y=1.25*np.sin(t1), line_color='black')
fig.add_scatter(x=1.75*np.cos(t1), y=1.75*np.sin(t1), fill='tonexty', line_color='black', fillcolor = 'steelblue')
fig.add_scatter(x=np.array([1.25,1.75])*np.cos(np.pi/8), y=np.array([1.25,1.75])*np.sin(np.pi/8), line_color='black',
mode='lines')
fig.add_scatter(x=np.array([1.25,1.75])*np.cos(-np.pi/8), y=np.array([1.25,1.75])*np.sin(-np.pi/8), line_color='black',
mode='lines')
fig.add_annotation(x=1.2, y=0, text=r'$\require{color}{\color[rgb]{0.315209,0.728565,0.037706}p}$', ax=0.5, ay=0, xref='x',
yref='y', axref='x', ayref='y', arrowhead=5)
fig.add_annotation(x=1.8, y=0, text=r"$\require{color} {\color[rgb]{0.315209,0.728565,0.037706}p} + \
\frac{d {\color[rgb]{0.315209,0.728565,0.037706}p}}{dr} \delta r$", ax=3, ay=0, xref='x', yref='y',
axref='x', ayref='y', arrowhead=5)
fig.add_annotation(x=0.6, y=0.75, ax=-0.05, ay=-0.09, xref='x', yref='y', axref='x', ayref='y',arrowhead=5)
fig.add_scatter(x=[0.25], y=[0.5], mode='text', text=r'$r_1$')
fig.add_annotation(x=0.6, y=-1.9, ax=-0.05, ay=0.09, xref='x', yref='y', axref='x', ayref='y', arrowhead=5)
fig.add_scatter(x=[0.25], y=[-0.5], mode='text', text=r'$r_2$')
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 how the velocity changes with radius we need to use Newton’s second law of motion. The only force acting on the flow is pressure since the flow is inviscid and therefore does not experience shear force. The flow changes direction as it moves around the bend meaning there must be a radial force causing the flow to move in a circular path around the origin. This force is centripetal and is given by the acceleration
\[ \large a_{c} = -\frac{{\color[rgb]{0.059472,0.501943,0.998465}v}^{2}(r)}{r} \hat{r}. \]
To find the mass we can use the fact that mass is equivalent to density times volume, when done per unit width this is
\[ \large \require{color}{\color[rgb]{0.501961,0.250953,0.010028}m}={\color[rgb]{0.814433,0.253157,0.091125}\rho}(\delta l \delta r) \]
Then we can use the force diagram above and write out Newton’s second law
\[ \large F=\require{color}{\color[rgb]{0.501961,0.250953,0.010028}m}a=-{\color[rgb]{0.814433,0.253157,0.091125}\rho}(\delta l \delta r) \frac{{\color[rgb]{0.059472,0.501943,0.998465}v}^{2}(r)}{r} = {\color[rgb]{0.315209,0.728565,0.037706}p} \delta l -\left({\color[rgb]{0.315209,0.728565,0.037706}p} + \frac{d{\color[rgb]{0.315209,0.728565,0.037706}p}}{dr} \delta r\right)\delta l. \]
This can then be solved for
\[ \large \require{color}\frac{d{\color[rgb]{0.315209,0.728565,0.037706}p}}{dr} \]
yielding
\[ \large \require{color}\frac{d{\color[rgb]{0.315209,0.728565,0.037706}p}}{dr} = {\color[rgb]{0.814433,0.253157,0.091125}\rho} \frac{{\color[rgb]{0.059472,0.501943,0.998465}v}^{2}(r)}{r}. \]
From this we can see that the normal pressure gradient is what causes the flow to move in a circular path around the bend. To solve for velocity we will need Bernoulli’s equation:
\[ \large \require{color}{\color[rgb]{0.315209,0.728565,0.037706}p}_{0}={\color[rgb]{0.315209,0.728565,0.037706}p}+\frac{1}{2}{\color[rgb]{0.814433,0.253157,0.091125}\rho} {\color[rgb]{0.059472,0.501943,0.998465}v}^{2} \]
where \(\require{color}{\color[rgb]{0.315209,0.728565,0.037706}p}_{0}\) is a constant. To relate this equation to the pressure gradient equation we just solved for we must take the derivative of Bernoulli’s equation so that it involves the pressure gradient. The derivative should be taken with respect to the radius as that is the only variable that changes velocity and pressure, remembering that this is incompressible flow we get
\[ \large \require{color}0=\frac{d{\color[rgb]{0.315209,0.728565,0.037706}p}}{dr}+ {\color[rgb]{0.814433,0.253157,0.091125}\rho} {\color[rgb]{0.059472,0.501943,0.998465}v}\frac{d{\color[rgb]{0.059472,0.501943,0.998465}v}}{dr}. \].
Plugging in our previous result for \[ \large \require{color}\frac{d{\color[rgb]{0.315209,0.728565,0.037706}p}}{dr} \]
we have
\[ \large \require{color}0={\color[rgb]{0.814433,0.253157,0.091125}\rho} \frac{{\color[rgb]{0.059472,0.501943,0.998465}v}^{2}(r)}{r} + {\color[rgb]{0.814433,0.253157,0.091125}\rho} {\color[rgb]{0.059472,0.501943,0.998465}v}(r) \frac{d{\color[rgb]{0.059472,0.501943,0.998465}v}(r)}{dr}. \]
Next, we can separate the varaibles and integrate both sides:
\[ \large \require{color}-\frac{dr}{r} = \frac{d{\color[rgb]{0.059472,0.501943,0.998465}v}(r)}{{\color[rgb]{0.059472,0.501943,0.998465}v}(r)} \longrightarrow {\color[rgb]{0.059472,0.501943,0.998465}v}(r)=\frac{C}{r} \]
where \(C\) is a constant to be determined by initial conditions.
From this we can see that the velocity is only dependent on the radius so velocity will be constant along at each radius, as we expected. Using the boundary conditions that
\[ \large \require{color}{\color[rgb]{0.059472,0.501943,0.998465}v}(r_1)={\color[rgb]{0.059472,0.501943,0.998465}v}_0 \]
we can get that
\[ \large \require{color}C={\color[rgb]{0.059472,0.501943,0.998465}v}_{0}r_{1}. \]
Therefore
\[ \large \require{color}{\color[rgb]{0.059472,0.501943,0.998465}v}(r)={\color[rgb]{0.059472,0.501943,0.998465}v}_{0}\frac{r_{1}}{r} \]
where \(r_{1}\) is the inner radius of the bend.
Let us now define the necessary variables to calculate the speed at each radii.
Code
# define necessary variables to calculate the speed at each radii
v_0 = 100 # initial speed at r_1
r_1 = 1 # inner radius
r_2 = 2 # outer radius
# calculate the speed between r_1 and r_2
r = np.linspace(r_1, r_2, 9)
v = v_0*r_1*(1/r)To solve for pressure we must use Bernoulli’s equation again. If we assume that at \(r_1\) the static pressure is atmospheric \(\require{color}{\color[rgb]{0.315209,0.728565,0.037706}p}_a = 101.325 \ kPa\) then we can find the stagnation pressure, this stagnation pressure remains constant throughout the bend. To do so we use Bernoulli’s equation at \(r_1\) as we know both the pressure and velocity at that point, using those values we get
\[ \large \require{color}{\color[rgb]{0.315209,0.728565,0.037706}p}_0 = {\color[rgb]{0.315209,0.728565,0.037706}p}_{\color[rgb]{0.989013,0.435749,0.811750}a} + \frac{1}{2}{\color[rgb]{0.814433,0.253157,0.091125}\rho} {\color[rgb]{0.059472,0.501943,0.998465}v}_0^2. \]
Now we can solve for the pressure at an radius by using the stagnation pressure we just solved for and the equation for velocity that we previously solved for:
\[ \large \require{color}{\color[rgb]{0.315209,0.728565,0.037706}p}(r) = {\color[rgb]{0.315209,0.728565,0.037706}p}_0 - \frac{1}{2}{\color[rgb]{0.814433,0.253157,0.091125}\rho} {\color[rgb]{0.059472,0.501943,0.998465}v}^2(r) = {\color[rgb]{0.315209,0.728565,0.037706}p}_0 - \frac{1}{2}{\color[rgb]{0.814433,0.253157,0.091125}\rho} \frac{{\color[rgb]{0.059472,0.501943,0.998465}v}_0^2r_1^2}{r^2} \]
\[ \large \Rightarrow {\color[rgb]{0.315209,0.728565,0.037706}p} = {\color[rgb]{0.315209,0.728565,0.037706}p}_a + \frac{1}{2} {\color[rgb]{0.814433,0.253157,0.091125}\rho} {\color[rgb]{0.059472,0.501943,0.998465}v}_0^2\left(1-\frac{r_1^2}{r^2}\right) \]
Next, we define some variables required to work out the pressure
Code
# define variables needed to calculate pressure
rho = 1.293 # denisty of air in kg/m^3
p_a = 101325 # pressure of air (atmospheric pressure) in Pa
p_0 = p_a + (1/2)*rho*v_0**2
# calculate pressure between r_1 and r_2
p = p_0 - (1/2)*rho*(v_0**2)*(r_1**2)/(r**2)
# define variables needed to make a contour plot
theta = np.linspace(0, np.pi/2, 50)
x=r_2*np.cos(theta)
y=np.linspace(-2,2,50)Visualizing the Flow
We are now ready to visualize the velocity and pressure contours!
Code
# define pressure and velocity for contour plot
v_z = []
p_z = []
for y_val in y:
v_y_list = []
p_y_list = []
for x_val in x:
if (x_val**2 + y_val**2) < r_1:
p_xy = p_a
v_xy = v_0
else:
v_xy = v_0/np.sqrt(x_val**2 + y_val**2)
p_xy = p_0 - (1/2)*rho*(v_0**2)*(r_1**2)/(x_val**2 + y_val**2)
v_y_list.append(v_xy)
p_y_list.append(p_xy)
v_z.append(v_y_list)
p_z.append(p_y_list)Finally, plotting this yields:
Code
# make contour plots for velocity and pressure
fig = go.Figure()
fig.add_contour(z=v_z, x=x, y=y, contours_coloring='heatmap', line_width=0, colorbar=dict(title='Velocity (m/s)',
titleside='right'), colorscale='PuBuGn', reversescale=True)
fig.add_contour(z=p_z, x=x, y=y, contours_coloring='heatmap', line_width=0, visible=False, colorbar=dict(title=
'Pressure (Pa)', titleside='right'), colorscale='dense', reversescale=True)
fig.add_shape(type='circle', xref='x', yref='y', x0=-1, y0=-1, x1=1, y1=1, fillcolor='white', line_color='white')
fig.add_scatter(x=2*np.cos(t2), y=2*np.sin(t2), line_width=0)
fig.add_scatter(x=3*np.cos(t2), y=3*np.sin(t2), line_color='white', fill='tonexty', fillcolor='white')
fig.update_layout(yaxis=dict(scaleanchor="x", scaleratio=1, showticklabels=False, range=[-2,2]),
xaxis=dict(showticklabels=False), plot_bgcolor='rgba(0, 0, 0, 0)', autosize=False, width=800, height=800,
showlegend=False, updatemenus=[dict(active=0, buttons=[dict(label='Velocity', method='update',
args=[{'visible':[True, False, True, True, True]}]), dict(label='Pressure', method='update',
args=[{'visible':[False, True, True, True, True]}])], y=1.2, x=0.41)])
fig.show()As you can see in the contour plots above the pressure and velocity are inversely proportional to eachother which is what we would expect based on Bernoulli’s equation. Additionally, the velocity decreases with increasing radius which we would also expect because centripetal acceleration decreases as the radius of curvative increases.
Code
# define lines where we want to visualize the flow and time and distance information for those lines
total_distance = np.pi*r
total_time = total_distance / v
dist = total_time[7] * v
r_lines = np.array([r[1], r[4], r[7]])
r_dist = np.array([dist[1], dist[4], dist[7]])
rad = -np.pi/2 + (r_dist / r_lines)
# create animated plot to show ideallized flow movement
fig = go.Figure()
frames = [go.Frame(data=[]) for k in range(50)]
for i, r_line in enumerate(r_lines):
fig.add_scatter(x=r_line*np.cos(t2), y=r_line*np.sin(t2), line_dash='dash', line_color='lightsteelblue')
x_pos = r_line*np.cos(np.linspace(-np.pi/2, rad[i], 50))
y_pos = r_line*np.sin(np.linspace(-np.pi/2, rad[i], 50))
try:
end = np.where(x_pos < 0)[0][0]
x_pos[end:len(x_pos)] = 'NaN'
except IndexError:
None
for j, f in enumerate(frames):
f.data += (go.Scatter(x=[x_pos[j]], y=[y_pos[j]], 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':100}}])])], xaxis=dict(range=[0,2], showticklabels=False),
yaxis=dict(scaleanchor="x", scaleratio=1, showticklabels=False), showlegend=False,
plot_bgcolor='rgba(0, 0, 0, 0)', autosize=False, width=900, height=700)
fig.add_scatter(x=2*np.cos(t2), y=2*np.sin(t2), fill='tozeroy', fillcolor='rgba(133, 163, 201, 0.5)', line_color='black')
fig.add_scatter(x=np.cos(t2), y=np.sin(t2), fill='tozeroy', fillcolor='white', line_color='black')
fig.show()In the above animation we can visually see how an idealized fluid would move around a bend, and we can see how the particles at inner radii move faster than those at outer radii