Starting from what is written here, we will try to modularize our code using function definitions. This will make things easier to manage.
We start by importing the libraries that we will be using in our code: the NumPy array library, the Matplotlib plotting library and the mathematical functions in the math
module.
1 2 3 4 | import numpy as np import math as mt from matplotlib import pyplot from matplotlib import cm |
To visualize the streamlines, we need to create a grid of points where we’ll compute the velocity.
We have our set of points now, and the two arrays X
and Y
contain their x– and y– coordinates (respectively) of every point on the rectangular grid.
Source in a uniform flow
We will first superimpose a source on a uniform flow and see what happens.
The streamlines of a freestream with speed U_\infty and angle of attack α are given by:
\psi_{\text {freestream }}(x, y)=U_{\infty}(y \cos \alpha-x \sin \alpha)
Note: the streamlines are all straight, parallel lines that make an angle α with the x-axis. If the flow is completely horizontal, ψ=U∞y. Integrate, and you get that u=U∞ and v=0.
Let’s write some code that will fill the arrays containing the u-velocity, the v-velocity and the stream function of a uniform horizontal flow (U∞,α=0), on every point of our grid. Note the handy NumPy functions ones()
, which creates a new array and fills it up with the value 1 everywhere, and zeros()
, which creates an array filled with 0.
The stream function of a source flow located at (xsource,ysource) is:
\psi_{\text {source }}(x, y)=\frac{\sigma}{2 \pi} \arctan \left(\frac{y-y_{\text {source }}}{x-x_{\text {source }}}\right)
and the velocity components are:

And remember that the stream function and velocity field of a source and a sink are exactly the same except one has positive strength while the other has negative strength.
We can write functions that serve a double purpose: with σ positive, they give the velocity and stream function of a source; with σ negative, they give them for a sink.
1 2 3 4 5 6 | def get_velocity(strength, xs, ys, X, Y): u = strength / ( 2 * np.pi) * (X - xs) / ((X - xs) * * 2 + (Y - ys) * * 2 ) v = strength / ( 2 * np.pi) * (Y - ys) / ((X - xs) * * 2 + (Y - ys) * * 2 ) return u, v |
Note that the output of the function consists of two arrays: u
and v
. They are calculated inside the function, which is indicated by the indentation of the lines after the colon. The final line indicates with the return
keyword that the arrays u, v
are sent back to the statement that called the function.
Similarly, we define another function to compute the stream-function of the singularity (source or sink) on the mesh grid, and call it get_stream_function()
.
1 2 3 4 5 | def get_stream_function(strength, xs, ys, X, Y): psi = strength / ( 2 * np.pi) * np.arctan2((Y - ys), (X - xs)) return psi |
Let’s use our brand new functions for the first time:
1 2 3 4 5 6 7 8 | strength_source = 5.0 # strength of the source x_source, y_source = - 1.0 , 0.0 # location of the source # compute the velocity field u_source, v_source = get_velocity(strength_source, x_source, y_source, X, Y) # compute the stream-function psi_source = get_stream_function(strength_source, x_source, y_source, X, Y) |
Let’s again use our superposition powers. The streamlines of the combination of a freestream and a source flow are:
\psi=\psi_{\text {freestream }}+\psi_{\text {source }}=U_{\infty} y+\frac{\sigma}{2 \pi} \arctan \left(\frac{y-y_{\text {source }}}{x-x_{\text {source }}}\right)
And since differentiation is linear, the velocity field induced by the new flow pattern is simply the sum of the freestream velocity field and the source velocity field:

The stagnation points in the flow are points where the velocity is zero. To find their location, we solve the following equations:
u=0, v=0
which leads to:

The streamline containing the stagnation point is called the dividing streamline. It separates the fluid coming from the freestream and the fluid radiating from the source flow. On the streamline plot, we’ll add a red curve to show the dividing streamline, and we’ll use the contour()
function for that.
We will also draw a red circle to show the location of the stagnation point, using the scatter()
function.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 | # plot the streamlines width = 10 height = (y_end - y_start) / (x_end - x_start) * width pyplot.figure(figsize = (width, height)) 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' ) # calculate the stagnation point x_stagnation = x_source - strength_source / ( 2 * np.pi * u_inf) y_stagnation = y_source # display the stagnation point pyplot.scatter(x_stagnation, y_stagnation, color = 'g' , s = 80 , marker = 'o' ) # display the dividing streamline pyplot.contour(X, Y, psi, levels = [ - strength_source / 2 , strength_source / 2 ], colors = '#CD2305' , linewidths = 2 , linestyles = 'solid' ); pyplot.show() |

If we ignore the flow inside the dividing streamline, we can consider that a solid body. In fact, this body has a name: it is called a Rankine half body.
Source-sink pair in a uniform flow
Now we can add a sink to our flow pattern without too much extra coding.
1 2 3 4 5 6 7 8 | strength_sink = - 5.0 # strength of the sink x_sink, y_sink = 1.0 , 0.0 # location of the sink # compute the velocity field on the mesh grid u_sink, v_sink = get_velocity(strength_sink, x_sink, y_sink, X, Y) # compute the stream-function on the grid mesh psi_sink = get_stream_function(strength_sink, x_sink, y_sink, X, Y) |
The superposition of the freestream, the source and the sink is just a simple addition.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 | # superposition of a source and a sink on the freestream u = u_freestream + u_source + u_sink v = v_freestream + v_source + v_sink psi = psi_freestream + psi_source + psi_sink # plot the streamlines width = 10 height = (y_end - y_start) / (x_end - x_start) * width pyplot.figure(figsize = (width, height)) 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, x_sink], [y_source, y_sink], color = '#CD2305' , s = 80 , marker = 'o' ) pyplot.contour(X, Y, psi, levels = [ 0. ], colors = '#CD2305' , linewidths = 2 , linestyles = 'solid' ); |

We can can look at the elliptical closed streamline as a solid surface and imagine that this is the flow around an egg-shaped object. It is called a Rankine oval.
Bernoulli’s equation and the pressure coefficient
A very useful measurement of a flow around a body is the coefficient of pressure Cp. To evaluate the pressure coefficient, we apply Bernoulli’s equation for an incompressible flow:

We define the pressure coefficient in the following way:

i.e.,

Note that in an incompressible flow, Cp=1 at a stagnation point.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 | # compute the pressure coefficient field cp = 1.0 - (u * * 2 + v * * 2 ) / u_inf * * 2 # plot the pressure coefficient field width = 10 height = (y_end - y_start) / (x_end - x_start) * width pyplot.figure(figsize = ( 1.1 * width, height)) pyplot.xlabel( 'x' , fontsize = 16 ) pyplot.ylabel( 'y' , fontsize = 16 ) pyplot.xlim(x_start, x_end) pyplot.ylim(y_start, y_end) contf = pyplot.contourf(X, Y, cp, levels = numpy.linspace( - 2.0 , 1.0 , 100 ), extend = 'both' , cmap = cm.gray) cbar = pyplot.colorbar(contf) cbar.set_label( '$C_p$' , fontsize = 16 ) cbar.set_ticks([ - 2.0 , - 1.0 , 0.0 , 1.0 ]) pyplot.scatter([x_source, x_sink], [y_source, y_sink], color = '#CD2305' , s = 80 , marker = 'o' ) pyplot.contour(X, Y, psi, levels = [ 0. ], colors = '#CD2305' , linewidths = 2 , linestyles = 'solid' ); |
