Version 0.2 -- February 2014¶
Lift on a cylinder¶
Remember when we computed uniform flow past a doublet? The stream-line pattern produced flow around a cylinder. When studying the pressure coefficient, we realized that the drag on the cylinder was exactly zero, leading to the D'Alembert paradox.
What about lift? Is it possible for a perfectly circular cylinder to experience lift? What if the cylinder is rotating? Have you heard about the Magnus effect?
You might be surprised to learn that all we need to do is add a vortex in the center of the cylinder. Let's see how that looks.
First, we recall the equations for the flow of a doublet. In Cartesian coordinates, a doublet located at the origin has a stream function and velocity components given by
$$\psi\left(x,y\right) = -\frac{\kappa}{2\pi}\frac{y}{x^2+y^2}$$$$u\left(x,y\right) = \frac{\partial\psi}{\partial y} = -\frac{\kappa}{2\pi}\frac{x^2-y^2}{\left(x^2+y^2\right)^2}$$$$v\left(x,y\right) = -\frac{\partial\psi}{\partial x} = -\frac{\kappa}{2\pi}\frac{2xy}{\left(x^2+y^2\right)^2}$$Let's start computing!¶
We'll place a doublet of strength $\kappa=1$ at the origin, and add a free stream $U_\infty=1$ (yes, we really like the number one). We can re-use the code we have written before; this is always a good thing.
from math import *
import numpy as np
import matplotlib.pyplot as plt
N = 50 # Number of points in each direction
x_start, x_end = -2.0, 2.0 # x-direction boundaries
y_start, y_end = -1.0, 1.0 # y-direction boundaries
x = np.linspace(x_start, x_end, N) # computes a 1D-array for x
y = np.linspace(y_start, y_end, N) # computes a 1D-array for y
X, Y = np.meshgrid(x, y) # generates a mesh grid
kappa = 1.0 # strength of the doublet
x_doublet, y_doublet = 0.0, 0.0 # location of the doublet
u_inf = 1.0 # freestream speed
Here are our function definitions for the doublet:
def get_velocity_doublet(strength, xd, yd, X, Y):
"""Returns the velocity field generated by a doublet.
Arguments
---------
strength -- strength of the doublet.
xd, yd -- coordinates of the doublet.
X, Y -- mesh grid.
"""
u = - strength/(2*pi)*((X-xd)**2-(Y-yd)**2)/((X-xd)**2+(Y-yd)**2)**2
v = - strength/(2*pi)*2*(X-xd)*(Y-yd)/((X-xd)**2+(Y-yd)**2)**2
return u, v
def get_stream_function_doublet(strength, xd, yd, X, Y):
"""Returns the stream-function generated by a doublet.
Arguments
---------
strength -- strength of the doublet.
xd, yd -- coordinates of the doublet.
X, Y -- mesh grid.
"""
psi = - strength/(2*pi)*(Y-yd)/((X-xd)**2+(Y-yd)**2)
return psi
And now we compute everything to get the flow around a cylinder, by adding a free stream to the doublet:
# computes the velocity field on the mesh grid
u_doublet, v_doublet = get_velocity_doublet(kappa, x_doublet, y_doublet, X, Y)
# computes the stream-function on the mesh grid
psi_doublet = get_stream_function_doublet(kappa, x_doublet, y_doublet, X, Y)
# freestream velocity components
u_freestream = u_inf * np.ones((N, N), dtype=float)
v_freestream = np.zeros((N, N), dtype=float)
# stream-function of the freestream flow
psi_freestream = u_inf * Y
# superposition of the doublet on the freestream flow
u = u_freestream + u_doublet
v = v_freestream + v_doublet
psi = psi_freestream + psi_doublet
We are ready to do a nice visualization. To get the plots inside the notebook, we use the command %matplotlib inline :
%matplotlib inline
# plots the streamlines
size = 10
plt.figure(figsize=(size, (y_end-y_start)/(x_end-x_start)*size))
plt.xlabel('x', fontsize=16)
plt.ylabel('y', fontsize=16)
plt.xlim(x_start, x_end)
plt.ylim(y_start, y_end)
plt.streamplot(X, Y, u, v, density=2, linewidth=1, arrowsize=1, arrowstyle='->')
plt.scatter(x_doublet, y_doublet, color='#CD2305', s=80, marker='o')
# calculates the cylinder radius and adds it to the figure
R = sqrt(kappa/(2*pi*u_inf))
circle = plt.Circle((0, 0), radius=R, color='#CD2305', alpha=0.5)
plt.gca().add_patch(circle)
# calculates the stagnation points and adds it to the figure
x_stagn1, y_stagn1 = +sqrt(kappa/(2*pi*u_inf)), 0
x_stagn2, y_stagn2 = -sqrt(kappa/(2*pi*u_inf)), 0
plt.scatter([x_stagn1, x_stagn2], [y_stagn1, y_stagn2], color='g', s=80, marker='o');
Nice! We have cylinder flow.
Now, let's add a vortex located at the origin with a positive strength $\Gamma$. In Cartesian coordinates, the stream function and velocity components are given by:
$$\psi\left(x,y\right) = \frac{\Gamma}{4\pi}\ln\left(x^2+y^2\right)$$$$u\left(x,y\right) = \frac{\Gamma}{2\pi}\frac{y}{x^2+y^2} \qquad v\left(x,y\right) = -\frac{\Gamma}{2\pi}\frac{x}{x^2+y^2}$$Based on these equations, we define the functions get_velocity_vortex() and get_stream_function_vortex() to do ... well, what's obvious by the function names (you should always try to come up with obvious function names). Play around with the value of $\ \Gamma$ and recalculate the flow. See what happens.
gamma = 4.0 # strength of the vortex
x_vortex, y_vortex = 0.0, 0.0 # location of the vortex
def get_velocity_vortex(strength, xv, yv, X, Y):
"""Returns the velocity field generated by a vortex.
Arguments
---------
strength -- strength of the vortex.
xv, yv -- coordinates of the vortex.
X, Y -- mesh grid.
"""
u = + strength/(2*pi)*(Y-yv)/((X-xv)**2+(Y-yv)**2)
v = - strength/(2*pi)*(X-xv)/((X-xv)**2+(Y-yv)**2)
return u, v
def get_stream_function_vortex(strength, xv, yv, X, Y):
"""Returns the stream-function generated by a vortex.
Arguments
---------
strength -- strength of the vortex.
xv, yv -- coordinates of the vortex.
X, Y -- mesh grid.
"""
psi = strength/(4*pi)*np.log((X-xv)**2+(Y-yv)**2)
return psi
# computes the velocity field on the mesh grid
u_vortex, v_vortex = get_velocity_vortex(gamma, x_vortex, y_vortex, X, Y)
# computes the stream-function on the mesh grid
psi_vortex = get_stream_function_vortex(gamma, x_vortex, y_vortex, X, Y)
Now that we have all the necessary ingredients (uniform flow, doublet and vortex), we apply the principle of superposition, and then we make a nice plot.
# superposition of the doublet and the vortex on the freestream flow
u = u_freestream + u_doublet + u_vortex
v = v_freestream + v_doublet + v_vortex
psi = psi_freestream + psi_doublet + psi_vortex
# calculates cylinder radius
R = sqrt(kappa/(2*pi*u_inf))
# calculates the stagnation points
x_stagn1, y_stagn1 = +sqrt(R**2-(gamma/(4*pi*u_inf))**2), -gamma/(4*pi*u_inf)
x_stagn2, y_stagn2 = -sqrt(R**2-(gamma/(4*pi*u_inf))**2), -gamma/(4*pi*u_inf)
# plots the streamlines
size = 10
plt.figure(figsize=(size, (y_end-y_start)/(x_end-x_start)*size))
plt.xlabel('x', fontsize=16)
plt.ylabel('y', fontsize=16)
plt.xlim(x_start, x_end)
plt.ylim(y_start, y_end)
plt.streamplot(X, Y, u, v, density=2, linewidth=1, arrowsize=1.5, arrowstyle='->')
circle = plt.Circle((0, 0), radius=R, color='#CD2305', alpha=0.5)
plt.gca().add_patch(circle)
plt.scatter(x_vortex, y_vortex, color='#CD2305', s=80, marker='o')
plt.scatter([x_stagn1, x_stagn2], [y_stagn1, y_stagn2], color='g', s=80, marker='o');
Challenge task¶
The challenge task in the doublet notebook was to calculate the radius of the cylinder created by the doublet in a uniform flow. You should have gotten
$$R = \sqrt{\frac{\kappa}{2\pi U_\infty}}$$The new challenge is to find where the stagnation points are located on the surface of the cylinder, when there's a vortex. (You just need an expression for the angles.)
What hapens if $\ \frac{\Gamma}{4\pi U_\infty R} >1$?
Go back and experiment with a value of $\Gamma$ that causes this.
Pressure coefficient¶
Let's get the pressure coefficient on the surface of the cylinder and compare with the case with no vortex. The velocity components in polar coordinates for the combined freestream + doublet + vortex are given by
$$u_r\left(r,\theta\right) = U_\infty \cos\theta \left(1-\frac{R^2}{r^2}\right)$$$$u_\theta\left(r,\theta\right) = -U_\infty \sin\theta \left(1-\frac{R^2}{r^2}\right) - \frac{\Gamma}{2\pi r}$$where $R$ is the cylinder radius.
We see that the radial component vanishes on the surface of the cylinder whereas the tangential velocity is given by
$$u_\theta\left(R,\theta\right) = -2U_\infty \sin\theta - \frac{\Gamma}{2\pi R} .$$As a note, when there is no vortex, the tangential velocity on the cylinder becomes
$$u_\theta\left(R,\theta\right) = -2U_\infty \sin\theta .$$From the doublet notebook, we know that the pressure coefficient is defined by
$$C_p = 1-\frac{U^2}{U_\infty^2}$$where $U^2 = u^2+v^2 = u_r^2+u_\theta^2$.
Let's plot it!
# calculates the surface tangential velocity on the cylinder
theta = np.linspace(0, 2*pi, 100)
u_theta = -2*u_inf*np.sin(theta) - gamma/(2*pi*R)
# computes the surface pressure coefficient
cp = 1.0 - (u_theta/u_inf)**2
# if there was no vortex
u_theta_no_vortex = -2*u_inf*np.sin(theta)
cp_no_vortex = 1.0 - (u_theta_no_vortex/u_inf)**2
# plots the surface pressure coefficient
size = 6
plt.figure(figsize=(size, size))
plt.grid(True)
plt.xlabel('$\theta$', fontsize=18)
plt.ylabel('$C_p$', fontsize=18)
plt.xlim(theta.min(), theta.max())
plt.plot(theta, cp, color='#CD2305', linewidth=2, linestyle='-')
plt.plot(theta, cp_no_vortex, color='g', linewidth=2, linestyle='-')
plt.legend(['with vortex', 'without vortex'], loc='best', prop={'size':16});
Lift and Drag¶
The lift is the component of force perpendicular to $U_\infty$, while the drag is the component parallel to $U_\infty$. How could we get them with the information we have above?
Well, the force on the cylinder is a product of the pressure acting on its surface (there is no viscosity here: it's ideal flow). If you draw a free body diagram, you should see that:
$$D = -\int_0^{2\pi} p \ \cos\theta \ R \ d\theta$$$$L = -\int_0^{2\pi} p \ \sin\theta \ R \ d\theta$$Challenge Task¶
Using Bernoulli's equation, replace $p$ in the equations above to obtain the lift and drag.
What does this mean?
The Magnus effect¶
The force experienced by a rotating cylinder (or sphere, or any object) is known as the Magnus effect.
Believe it or not, someone actually tried to build an airplane with this concept: spinning cylinders as "wings." According to an article on PilotFriend, a plane called the 921-V was built in 1930 and flew "at least once" before crashing.
from IPython.display import Image
Image(url='http://upload.wikimedia.org/wikipedia/commons/7/78/Flettner_Rotor_Aircraft.jpg')
And nowadays, a handful of hobbyists build RC "rotorwings" taking advantage of the Magnus effect to collect views on YouTube ...
from IPython.display import YouTubeVideo
YouTubeVideo('POHre1P_E1k')
from IPython.core.display import HTML
def css_styling():
styles = open('../styles/custom.css', 'r').read()
return HTML(styles)
css_styling()