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.

```
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 (*x*_{source},*y*_{source}) 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.

```
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()`

.

```
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:

```
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.

```
# 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**.

```
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.

```
# 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 C _{p}*. 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, *C _{p}*=1 at a stagnation point.

```
# 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');
```