In [1]:
import numpy
from matplotlib import pyplot
%matplotlib inline
we can do the same thing for the code in VortexPanel. All we need is:
In [2]:
from VortexPanel import Panel, solve_gamma, plot_flow
This makes the Panel class, and the solveGamma and plotFlow functions available. We can test out all three by solving for the flow around a circle again.
In [3]:
def make_circle(N):
# define the end-points of the panels
x_ends = numpy.cos(numpy.linspace(0, -2*numpy.pi, N+1))
y_ends = numpy.sin(numpy.linspace(0, -2*numpy.pi, N+1))
# define the panels
circle = numpy.empty(N, dtype=object)
for i in xrange(N):
circle[i] = Panel(x_ends[i], y_ends[i], x_ends[i+1], y_ends[i+1])
return circle
In [4]:
N = 16
circle = make_circle(N) # set-up geom
solve_gamma(circle) # find gamma
plot_flow(circle) # plot the flow
Looks good.
To find the lift on a foil, we first need to define the geometry. Let's use a symmetric Jukowski foil for now. This foil shape is obtained by applying the following mapping
$$ x' = \frac x2 (1+1/r^2), \quad y' = \frac y2 (1-1/r^2),\quad r^2= x^2+y^2 $$to the points on a circle which intersects $x,y=1,0$. Shifting the circle an amount $\Delta x$ sets the thickness of the foil.
In [5]:
def make_jukowski( N, dx = 0.2, dy = 0, dr = 0 ):
# define the circle
x_ends = numpy.cos(numpy.linspace(0, -2*numpy.pi, N+1))
y_ends = numpy.sin(numpy.linspace(0, -2*numpy.pi, N+1))
# shift circle
r = numpy.sqrt((1+dx)**2+dy**2)+dr
x2_ends = r*x_ends-dx
y2_ends = r*y_ends-dy
r2_ends = x2_ends**2+y2_ends**2
# apply jukowski mapping
x3_ends = x2_ends*(1+1./r2_ends)/2
y3_ends = y2_ends*(1-1./r2_ends)/2
# define the panels
foil = numpy.empty(N, dtype=object)
for i in xrange(N):
foil[i] = Panel(x3_ends[i], y3_ends[i], x3_ends[i+1], y3_ends[i+1])
return foil
In [6]:
foil = make_jukowski(N) # make foil
solve_gamma(foil) # solve for gamma
plot_flow(foil) # plot the flow
The lift on a body in potential flow is given by
$$ L = -\rho U_\infty \Gamma $$where $\rho$ is the density and $\Gamma$ is the total circulation. For a vortex sheet we remember this is given by $\Gamma=\int\gamma ds $, and so for a vortex panel method this is simply
$$ \Gamma = \sum_{i=0}^{N-1} \gamma_i 2S_i $$since the width of each panel is $2S_i$. The lift coefficient is then
$$ C_L =\frac L{\tfrac 12 \rho U_\infty^2 c} = -\frac 4{U_\infty c} \sum_{i=0}^{N-1} \gamma_i S_i $$where $c$ is the coord length, the distance from the leading to trailing edge.
This is very simple to implement, so let's find the lift on the foil:
In [7]:
def lift(panels):
c = panels[0].x[0]-panels[len(panels)/2].x[0] # length scale
return -4./c*numpy.sum([p.gamma*p.S for p in panels])
print 'C_L =',lift(foil)
Zero lift...
Oh, of course. We need an angle of attack to get lift. No problem:
In [8]:
alpha = numpy.pi/16 # set angle of attack
solve_gamma(foil,alpha) # solve for gamma
plot_flow(foil,alpha) # plot the flow
print 'C_L =',lift(foil) # print the lift
Sigh. There are so many things wrong with this picture...
These are all physcially incorrect!
What can we do to the circle flow to generate lift force?
Note that the addition of a vortex of strength $\Gamma^+$ to the center of the circle has exactly the same effect on the external flow as changing $\gamma_i$ by an amount $\Delta\gamma=\frac {\Gamma^+}{2\pi R}$.
Try it out with the code below:
In [9]:
alpha = 0 # angle of attack
dgamma = 0./(2*numpy.pi) # vortex circulation
N = 16 # number of panels
circle = make_circle(N) # set-up geom
solve_gamma(circle,alpha) # find gamma
for p in circle: p.gamma += dgamma # increase gamma by dgamma
plot_flow(circle,alpha) # plot the flow
print 'C_L =',lift(circle) # print the lift
We see that we can add any amount $\Delta\gamma$, and generate any lift we want. This should be unsettling in at least two different ways:
Let get the first issue out of the way first. Indeed, $\gamma_i$ is not unique. This is because the matrix $\mathbf A$ for the vortex panel system is singular.
We didn't notice this in the last notebook, because the discretization errors in construct_A_b nudge $\mathbf A$ out of singularity. This is kind of nice since it lets numpy.linalg.solve find a solution. But it is also the reason why the solution for the foil at an angle of attack above is completely bogus - the discretization error doesn't nudge you towards the correct solution in general.
The fundamental problem with the solution for the foil above is that the flow is wrapping around the sharp trailing edge. This means that particles traveling on the body streamline instantly change direction - which is impossible for any object with mass.
Avoiding this impossibility is called the Kutta condition.
This uniquely determines the correct solution $\gamma_i$ and therefore $\Gamma$ and the $C_L$.
Additionally, this links the angle of attack $\alpha$ to the lift. As $\alpha$ is increased, the amount of circulation increases to enforce the Kutta condition.
To enfore the Kutta condition in the vortex panel method, we need to convert it into a statement about the strength $\gamma_i$.
What condition on $\gamma_i$ will enforce the Kutta conditions?
(Hint: Look at the sketch above. What does it tell you about $\gamma$ near the trailing edge?)
Now we know the additional condition which makes the solution unique, but we still need to enforce it. This is done by adjusting the linear system.
The way we have assembled the geometries, the trailing edge is the first and last end point. This means the panels on either side of the edge are $i=0,N-1$. The Kutta condition is therefore
$$ \gamma_0 +\gamma_{N-1} = 0.$$The easiest (and most numerically stable) way to include this condition is to add this equation to the no-slip equation at every panel, ie
$$ \sum_{j=0}^{N-1} a_{ij} \gamma_j +\gamma_0 +\gamma_{N-1} = b_i \quad i=0\ldots N-1$$This gives a new linear system $\mathbf{A' \gamma = b}$, where $$ a_{ij}' = a_{ij}+1 \quad j=0,{N-1} $$ $$ a_{ij}' = a_{ij} \quad \mbox{else} $$
We can implement this in a new solveGammaKutta function, which is identical to solveGamma but with the change to the matrix given above.
In [10]:
# determine the vortex panel strength with Kutta Condition
def solve_gamma_kutta(panels,alpha=0):
from VortexPanel import construct_A_b
A,b = construct_A_b(panels,alpha) # construct linear system
A[:, 0] += 1 # gamma[0]+ ...
A[:,-1] += 1 # gamma[N-1]=0
gamma = numpy.linalg.solve(A, b) # solve for gamma!
for i,p_i in enumerate(panels):
p_i.gamma = gamma[i] # update panels
And that's it! Let's try it out:
In [11]:
N = 16
alpha = numpy.pi/32
foil = make_jukowski(N) # make foil
solve_gamma_kutta(foil,alpha) # solve for gamma
plot_flow(foil,alpha) # plot the flow
print 'C_L =',lift(foil) # print the lift
Compare to the figure above. This one has:
So, is the prediction any good?
The exact solution for the lift on a Jukowski foil is
$$C_L = 2\pi \left(1+\frac 4{3\sqrt 3} \frac tc \right)\sin\alpha$$where $t/c \approx \Delta x+0.019$ is the thickness to coord ratio. ie
In [16]:
def jukowski_CL(alpha,t_c=0.2+0.019): return 2*numpy.pi*(1+4/3/numpy.sqrt(3)*t_c)*numpy.sin(alpha)
In [13]:
# your code here
Ignore the line below - it just loads the style sheet.
In [14]:
from IPython.core.display import HTML
def css_styling():
styles = open('../styles/custom.css', 'r').read()
return HTML(styles)
css_styling()
Out[14]: