A source sheet is an infinite row of sources along one line. You have already created a vortex sheet by placing a row of vortices on one line, so you are ready for this.
Make sure that you have studied AeroPython Lesson 7, on the Method of Images and the use of classes in Python. From now on, classes will feature prominently!
We start by importing our libraries and creating the grid of points, as we've done so many times before.
In [1]:
import numpy
import math
from matplotlib import pyplot
In [2]:
N = 100 # Number of points in each direction
x_start, x_end = -1.0, 1.0 # x-direction boundaries
y_start, y_end = -1.5, 1.5 # y-direction boundaries
x = numpy.linspace(x_start, x_end, N) # computes a 1D-array for x
y = numpy.linspace(y_start, y_end, N) # computes a 1D-array for y
X, Y = numpy.meshgrid(x, y) # generates a mesh grid
We will include a uniform flow of magnutide $U_\infty = 1$ parallel to the horizontal axis. The arrays uFreestream and vFreestream contain the velocity components of the free stream. Let's fill these up:
In [3]:
u_inf = 1.0 # free-stream speed
# calculates the free-stream velocity components
u_freestream = u_inf * numpy.ones((N, N), dtype=float)
v_freestream = numpy.zeros((N, N), dtype=float)
We consider first a finite number of sources arranged along a vertical line, normal to the freestream. The streamlines will come out of each source and be deflected by the streamlines coming from the freestream.
From a computational point of view, the finite number of sources can be represented by a 1D-array containing objects, each one having been created by the class Source. This class will need as attributes the strength of a source, strength, and its location at (x,y). Its methods will compute the velocity components and stream function on a given mesh grid (X,Y), and of course, its constructor:
In [4]:
class Source:
"""Contains information related to a source/sink."""
def __init__(self, strength, x, y):
"""Initializes the source/sink.
Arguments
---------
strength -- strength of the source/sink.
x, y -- coordinates of the source/sink.
"""
self.strength = strength
self.x, self.y = x, y
def velocity(self, X, Y):
"""Computes the velocity field generated by the source/sink.
Arguments
---------
X, Y -- mesh grid.
"""
self.u = self.strength/(2*math.pi)*(X-self.x)/((X-self.x)**2+(Y-self.y)**2)
self.v = self.strength/(2*math.pi)*(Y-self.y)/((X-self.x)**2+(Y-self.y)**2)
def stream_function(self, X, Y):
"""Computes the stream-function generated by the source/sink.
Arguments
---------
X, Y -- mesh grid.
"""
self.psi = self.strength/(2*math.pi)*numpy.arctan2((Y-self.y), (X-self.x))
Now you'll really see how useful classes are! We will use our class Source to create as many sources as we want to place on a source sheet. Pick a number, N_sources, and call the class constructor that many times, moving along the $y$ axis for the location of the sources. Then call the function to calculate the velocity of each source.
In [5]:
N_sources = 11 # number of sources
strength = 5.0 # sum of all source strengths
strength_source = strength/N_sources # strength of one source
x_source = numpy.zeros(N_sources, dtype=float) # horizontal position of all sources (1D array)
y_source = numpy.linspace(-1.0, 1.0, N_sources) # vertical position of all sources (1D array)
# creates a source line (Numpy array of Source's instantiations)
sources = numpy.empty(N_sources, dtype=object)
for i in xrange(N_sources):
sources[i] = Source(strength_source, x_source[i], y_source[i])
sources[i].velocity(X, Y)
# superposition of all sources to the free-stream flow
u = u_freestream.copy()
v = v_freestream.copy()
for source in sources:
u += source.u
v += source.v
Notice that the variable sources is a NumPy array, that is, a collection of items of the same type. What type? You're probably used to arrays of numbers, but we can also have arrays of objects of any type created by a class (in this case, of type Source).
We are creating an empty NumPy array called sources, and telling Python that its items will be some object that is not a built-in data type, like int or float for integers or floating-point real numbers. The number of items contained in the array will be N_sources.
In the first loop, we fill the sources array, calling the Source-class constructor for each item. We also go ahead and compute the velocity of each source (using the method velocity() of each source). In the second loop—after creating the velocity arrays u and v by copying the free-stream velocity components—, we're adding all the velocity contributions for every source source in the array sources.
The final flow pattern corresponds to the superposition of a uniform flow and N_sources identical sources of strength strength_source equally spaced along a vertical line normal to the uniform flow. In our plot, we'll mark the location of the sources with red points and also define a filled contour of velocity magnitude to visualize the location of the stagnation point, defined by $u=0$ and $v=0$.
In [6]:
# plots the streamlines
%matplotlib inline
size = 6
pyplot.figure(figsize=(size, (y_end-y_start)/(x_end-x_start)*size))
pyplot.grid(True)
pyplot.xlabel('x', fontsize=16)
pyplot.ylabel('y', fontsize=16)
pyplot.xlim(x_start, x_end)
pyplot.ylim(y_start, y_end)
pyplot.streamplot(X, Y, u, v, density=2, linewidth=1, arrowsize=1, arrowstyle='->')
pyplot.scatter(x_source, y_source, color='#CD2305', s=80, marker='o')
velocity = pyplot.contourf(X, Y, numpy.sqrt(u**2+v**2), levels=numpy.linspace(0.0, 0.1, 10))
cbar = pyplot.colorbar(velocity, ticks=[0, 0.05, 0.1], orientation='horizontal')
cbar.set_label('Velocity magnitude', fontsize=16);
Now try to change the total strength of the sources. What is the minimal total source strength so that the dividing streamline ($\psi = 0$) encloses all the singularities?
By definition, a source sheet is an infinite collection of side-by-side sources of infinitesimal equal strength distributed along a given path.
Consider $s$ to be the running coordinate along the sheet: we define $\sigma = \sigma(s)$ as the strength per unit length.
Recall that the strength of a single source represents the volume flow rate per unit depth (i.e., per unit length in the $z$-direction). Therefore, $\sigma$ represents the volume flow rate per unit depth and per unit length (in the $s$-direction). Thus, $\sigma ds$ is the strength of an infinitesimal portion $ds$ of the source sheet. This infinitesimal portion is so small that it can be treated as a distinct source of strength $\sigma ds$.
Following this, the stream-function at point $\left(r,\theta\right)$ of this infinitesimal portion is
$$d\psi\left(r,\theta\right) = \frac{\sigma ds}{2\pi} \theta$$Integrating along $s$, we find the stream-function of the entire source sheet:
$$\psi\left(r,\theta\right) = \frac{\sigma}{2\pi}\int_{\text{sheet}}\theta\left(s\right)ds$$In the previous section, we considered a vertical finite distribution of sources. Similarly, the stream function of a vertical source sheet at $x=0$ between $y_{\text{min}}$ and $y_{\text{max}}$ in Cartesian coordinates is
$$\psi\left(x,y\right) = \frac{\sigma}{2\pi}\int_{y_{\text{min}}}^{y_{\text{max}}} \tan^{-1}\left(\frac{y-\xi}{x}\right)d\xi$$And the velocity components are
$$u\left(x,y\right) = \frac{\sigma}{2\pi} \int_{y_{\text{min}}}^{y_{\text{max}}} \frac{x}{x^2+\left(y-\xi\right)^2}d\xi$$$$v\left(x,y\right) = \frac{\sigma}{2\pi} \int_{y_{\text{min}}}^{y_{\text{max}}} \frac{y-\xi}{x^2+\left(y-\xi\right)^2}d\xi$$We need to calculate the two integrals above to obtain the velocities. Bummer, you say. Calculating integrals is so passé. Don't worry! We have SciPy: a powerful collection of mathematical algorithms and functions. It includes the module integrate, a fundamental set of functions for scientific computing including basic integration, quadrature rules and numerical integrators for ordinary differential equations. How helpful is that?
Let's import the module we need from SciPy:
In [7]:
from scipy import integrate
We are going to use the function quad(func,a,b) of the module integrate to compute the definite integral of a function of one variable between two limits:
The first argument to the function quad needs to be a Python function ... you already know how to create a Python function using a def statement, but you'll learn a new way to create Python functions: the lambda statement. Why do I need two ways to define functions, you may ask? The answer is rather subtle, and the subject of full blog posts! (See below)
In short, a function created with a lambda statement is a single expression that returns a value (but without a return statement!). It is often called an anonymous function, because we don't need to give it a name.
Let's see how we use a lambda statement to integrate a mathematical function. Suppose we want to integrate the function $f:x\rightarrow x^2$ from $0$ to $1$. You can manually do this one, right? It gives the value $\frac{1}{3}$.
To use the function quad(), we give it as first parameter the expression lambda x : x**2:
In [8]:
print integrate.quad(lambda x : x**2, 0.0, 1.0)
You see, here we used the lambda statement to pass a function argument to quad(), without going through a function definition (in this sense, we talk of anonymous Python functions).
Note that the function quad() returns a list: the first element is the result of the integral, and the second element is an estimate of the error. If you just want the value of the definite integral, you need to specify the index [0], for the first element of the list!
We are ready to use this for our velocity calculations of the source sheet.
In [9]:
sigma = 2.5 # strength of the source-sheet
# boundaries of the source-sheet
y_min, y_max = -1.0, 1.0
# computes the velcity field generated by the source-sheet
u_sheet = numpy.empty((N, N), dtype=float)
v_sheet = numpy.empty((N, N), dtype=float)
for i in xrange(N):
for j in xrange(N):
integrand = lambda s : X[i,j]/(X[i,j]**2+(Y[i,j]-s)**2)
u_sheet[i,j] = sigma/(2*math.pi)*integrate.quad(integrand, y_min, y_max)[0]
integrand = lambda s: (Y[i,j]-s)/(X[i,j]**2+(Y[i,j]-s)**2)
v_sheet[i,j] = sigma/(2*math.pi)*integrate.quad(integrand, y_min, y_max)[0]
# superposition of the source-sheet to the uniform flow
u = u_freestream + u_sheet
v = v_freestream + v_sheet
Let's now visualize the streamlines (blue) on the mesh grid. We draw a red line to show the source sheet and a filled contour to visualize the location of the stagnation point.
In [10]:
# plots the streamlines
%matplotlib inline
size = 8
pyplot.figure(figsize=(size, (y_end-y_start)/(x_end-x_start)*size))
pyplot.grid(True)
pyplot.xlabel('x', fontsize=16)
pyplot.ylabel('y', fontsize=16)
pyplot.xlim(x_start, x_end)
pyplot.ylim(y_start, y_end)
pyplot.streamplot(X, Y, u, v, density=2, linewidth=1, arrowsize=1, arrowstyle='->')
pyplot.axvline(0.0, (y_min-y_start)/(y_end-y_start), (y_max-y_start)/(y_end-y_start),
color='#CD2305', linewidth=4)
velocity = pyplot.contourf(X, Y, numpy.sqrt(u**2+v**2), levels=numpy.linspace(0.0, 0.1, 10))
cbar = pyplot.colorbar(velocity, ticks=[0, 0.05, 0.1], orientation='horizontal')
cbar.set_label('Velocity magnitude', fontsize=16);
We said that the strength $\sigma$ represents a volume flow rate emanating from the source sheet. If you play with this parameter, you will see that the stagnation point comes closer and closer to the source sheet with decreasing strength.
If we wanted to use multiple source sheets to represents the streamlines around a given body shape, we need to make each source sheet part of the dividing streamline. The question will be What is the source strength to make this happen?
The volume flow rate on just the left side of the sheet is $\frac{\sigma}{2}$, flowing in the opposite direction of the freestream velocity $U_\infty$. Therefore, the flow-tangency boundary condition required is $\frac{\sigma}{2} = U_\infty$.
Now go back to the code above an replace the sigma by the correct value. Where is the stagnation point? Where is the dividing streamline?
The lambda statement can be especially confusing when you are starting out with Python. Here are a couple of places to dig deeper:
In [11]:
from IPython.core.display import HTML
def css_styling():
styles = open('../styles/custom.css', 'r').read()
return HTML(styles)
css_styling()
Out[11]: