Diary

Lei crede in Dio? «Dio per me non è neanche un’ ipotesi». In un lontano incontro, Cabibbo disse che la fede è comunque un bel vantaggio. «Non dubito, se uno crede. Ma uno non è che può credere per vivere meglio. Cabibbo è stato Presidente dell’ Accademia Pontificia, ma non so quanto credesse».

• Da Giorgio Parisi: “Ecco perché la fisica moderna è assurda come un testo di Beckett“, la Repubblica, 31 dicembre 2010.

A doublet alone does not give so much information about how it can be used to represent a practical flow pattern in aerodynamics. But let’s use our superposition powers: our doublet in a uniform flow turns out to be a very interesting flow pattern. Let’s first define a uniform horizontal flow.

u_inf = 1.0        # freestream speed

Now, we can calculate velocities and stream function values for all points in our grid. And as we now know, we can calculate them all together with one line of code per array.

u_freestream = u_inf * numpy.ones((N, N), dtype=float)
v_freestream = numpy.zeros((N, N), dtype=float)

psi_freestream = u_inf * Y

Below, the stream function of the flow created by superposition of a doublet in a free stream is obtained by simple addition.

The plot shows that this pattern can represent the flow around a cylinder with center at the location of the doublet. All the streamlines remaining outside the cylinder originated from the uniform flow. All the streamlines inside the cylinder can be ignored and this area assumed to be a solid object. This will turn out to be more useful than you may think.

# superposition of the doublet on the freestream flow
u = u_freestream + u_doublet
v = v_freestream + v_doublet
psi = psi_freestream + psi_doublet

# 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.contour(X, Y, psi,
               levels=[0.], colors='#CD2305', linewidths=2, linestyles='solid')
pyplot.scatter(x_doublet, y_doublet, color='#CD2305', s=80, marker='o')

# calculate the stagnation points
x_stagn1, y_stagn1 = +math.sqrt(kappa / (2 * math.pi * u_inf)), 0.0
x_stagn2, y_stagn2 = -math.sqrt(kappa / (2 * math.pi * u_inf)), 0.0

# display the stagnation points
pyplot.scatter([x_stagn1, x_stagn2], [y_stagn1, y_stagn2],
               color='g', s=80, marker='o');

Bernoulli’s equation and the pressure coefficient

A very useful measurement of a flow around a body is the coefficient of pressureCp. To evaluate the pressure coefficient, we apply Bernoulli’s equation for ideal flow and with simple mathematical steps:

\(C_{p}=1-\left(\frac{U}{U_{\infty}}\right)^{2}\)

In an incompressible flow, C_p=1 at a stagnation point. Let’s plot the pressure coefficient in the whole domain.

# 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')
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_doublet, y_doublet,
               color='#CD2305', s=80, marker='o')
pyplot.contour(X,Y,psi,
               levels=[0.], colors='#CD2305', linewidths=2, linestyles='solid')
pyplot.scatter([x_stagn1, x_stagn2], [y_stagn1, y_stagn2],
               color='g', s=80, marker='o');
pyplot.streamplot(X, Y, u, v,
                  density=2, linewidth=1, arrowsize=1, arrowstyle='->')

Below, we report the complete Python code:

import math as mt
import numpy as np
from matplotlib import pyplot as pyp

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)    # creates a 1D-array for x
y = np.linspace(y_start, y_end, N)    # creates 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

def get_velocity_doublet(strength, xd, yd, X, Y):
    
    u = (- strength / (2 * mt.pi) *
         ((X - xd)**2 - (Y - yd)**2) /
         ((X - xd)**2 + (Y - yd)**2)**2)
    v = (- strength / (2 * mt.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):
    
    psi = - strength / (2 * mt.pi) * (Y - yd) / ((X - xd)**2 + (Y - yd)**2)
    
    return psi

# compute the velocity field on the mesh grid
u_doublet, v_doublet = get_velocity_doublet(kappa, x_doublet, y_doublet, X, Y)

# compute the stream-function on the mesh grid
psi_doublet = get_stream_function_doublet(kappa, x_doublet, y_doublet, X, Y)

# plot the streamlines
width = 10
height = (y_end - y_start) / (x_end - x_start) * width
pyp.figure(figsize=(width, height))
pyp.xlabel('x', fontsize=16)
pyp.ylabel('y', fontsize=16)
pyp.xlim(x_start, x_end)
pyp.ylim(y_start, y_end)
pyp.streamplot(X, Y, u_doublet, v_doublet,
                  density=2, linewidth=1, arrowsize=1, arrowstyle='->')
pyp.scatter(x_doublet, y_doublet, color='#CD2305', s=80, marker='o');



u_inf = 1.0        # freestream speed

u_freestream = u_inf * np.ones((N, N), dtype=float)
v_freestream = np.zeros((N, N), dtype=float)

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

# plot the streamlines
width = 10
height = (y_end - y_start) / (x_end - x_start) * width
pyp.figure(figsize=(width, height))
pyp.xlabel('x', fontsize=16)
pyp.ylabel('y', fontsize=16)
pyp.xlim(x_start, x_end)
pyp.ylim(y_start, y_end)
pyp.streamplot(X, Y, u, v,
                  density=2, linewidth=1, arrowsize=1, arrowstyle='->')
pyp.contour(X, Y, psi,
               levels=[0.], colors='#CD2305', linewidths=2, linestyles='solid')
pyp.scatter(x_doublet, y_doublet, color='#CD2305', s=80, marker='o')

# calculate the stagnation points
x_stagn1, y_stagn1 = +mt.sqrt(kappa / (2 * mt.pi * u_inf)), 0.0
x_stagn2, y_stagn2 = -mt.sqrt(kappa / (2 * mt.pi * u_inf)), 0.0

# display the stagnation points
pyp.scatter([x_stagn1, x_stagn2], [y_stagn1, y_stagn2],
               color='g', s=80, marker='o');

# 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
pyp.figure(figsize=(1.1 * width, height))
pyp.xlabel('x', fontsize=16)
pyp.ylabel('y', fontsize=16)
pyp.xlim(x_start, x_end)
pyp.ylim(y_start, y_end)
contf = pyp.contourf(X, Y, cp,
                        levels=np.linspace(-2.0, 1.0, 100), extend='both')
cbar = pyp.colorbar(contf)
cbar.set_label('$C_p$', fontsize=16)
cbar.set_ticks([-2.0, -1.0, 0.0, 1.0])
pyp.scatter(x_doublet, y_doublet,
               color='#CD2305', s=80, marker='o')
pyp.contour(X,Y,psi,
               levels=[0.], colors='#CD2305', linewidths=2, linestyles='solid')
pyp.scatter([x_stagn1, x_stagn2], [y_stagn1, y_stagn2],
               color='g', s=80, marker='o');
pyp.streamplot(X, Y, u, v,
                  density=2, linewidth=1, arrowsize=1, arrowstyle='->')


pyp.show()

Before we do the Python code, let’s talk about the heat equation and finite-difference method. Heat equation is basically a partial differential equation, it is

If we want to solve it in 2D (Cartesian), we can write the heat equation above like this:

where u is the quantity that we want to know, t is for temporal variable, x and y are for spatial variables, and α is diffusivity constant. So basically we want to find the solution u everywhere in x and y, and over time t.

We can write the heat equation above using finite-difference method like this:

If we arrange the equation above by taking Δx = Δy, we get this final equation:

where

\(\gamma = \alpha \frac{{\Delta t}}{{\Delta {x^2}}}\)

We use explicit method to get the solution for the heat equation, so it will be numerically stable whenever

\(\Delta t \leq \frac{{\Delta {x^2}}}{{4\alpha }}\)

Everything is ready. Now we can solve the original heat equation approximated by algebraic equation above, which is computer-friendly.

Let’s suppose a thin square plate with the side of 50 unit length. The temperature everywhere inside the plate is originally 0 degree (at t = 0), let’s see the diagram below:

For our model, let’s take Δx = 1 and α = 2.0.
Now we can use Python code to solve this problem numerically to see the temperature everywhere (denoted by i and j) and over time (denoted by k). Let’s first import all of the necessary libraries, and then set up the boundary and initial conditions.

import numpy as np
import matplotlib.pyplot as plt
import matplotlib.animation as animation
from matplotlib.animation import FuncAnimation

print("2D heat equation solver")

plate_length = 50
max_iter_time = 750

alpha = 2
delta_x = 1

delta_t = (delta_x ** 2)/(4 * alpha)
gamma = (alpha * delta_t) / (delta_x ** 2)

# Initialize solution: the grid of u(k, i, j)
u = np.empty((max_iter_time, plate_length, plate_length))

# Initial condition everywhere inside the grid
u_initial = 0

# Boundary conditions
u_top = 100.0
u_left = 0.0
u_bottom = 0.0
u_right = 0.0

# Set the initial condition
u.fill(u_initial)

# Set the boundary conditions
u[:, (plate_length-1):, :] = u_top
u[:, :, :1] = u_left
u[:, :1, 1:] = u_bottom
u[:, :, (plate_length-1):] = u_right

We’ve set up the initial and boundary conditions, let’s write the calculation function based on finite-difference method that we’ve derived above.

def calculate(u):
    for k in range(0, max_iter_time-1, 1):
        for i in range(1, plate_length-1, delta_x):
            for j in range(1, plate_length-1, delta_x):
                u[k + 1, i, j] = gamma * (u[k][i+1][j] + u[k][i-1][j] + u[k][i][j+1] + u[k][i][j-1] - 4*u[k][i][j]) + u[k][i][j]

    return u

Let’s prepare the plot function so we can visualize the solution (for each k) as a heat map. We use Matplotlib library, it’s easy to use.

def plotheatmap(u_k, k):
    # Clear the current plot figure
    plt.clf()

    plt.title(f"Temperature at t = {k*delta_t:.3f} unit time")
    plt.xlabel("x")
    plt.ylabel("y")

    # This is to plot u_k (u at time-step k)
    plt.pcolormesh(u_k, cmap=plt.cm.jet, vmin=0, vmax=100)
    plt.colorbar()

    return plt

One more thing that we need is to animate the result because we want to see the temperature points inside the plate change over time. So let’s create the function to animate the solution.

def animate(k):
  plotheatmap(u[k], k)

anim = animation.FuncAnimation(plt.figure(), animate, interval=1, frames=max_iter_time, repeat=False)
anim.save("heat_equation_solution.gif")

That’s it! And here’s the result:

Think about the Source & Sink again, and now imagine that you are looking at this flow pattern from very far away. The streamlines that are between the source and the sink will be very short, from this vantage point. And the other streamlines will start looking like two groups of circles, tangent at the origin. If you look from far enough away, the distance between source and sink approaches zero, and the pattern you see is called a doublet.

Let’s see what this looks like. First, load our favorite libraries.

First, consider a source of strength \(\Lambda\) at (−l₂,0) and a sink of opposite strength located at (l₂,0). At any point P in the flow, the stream funcion is

\(\psi(x, y)=\frac{\Lambda}{2 \pi}\left(\theta_{1}-\theta_{2}\right)=-\frac{\Lambda}{2 \pi} \Delta \theta \)

where \(\Delta \theta = \theta_2 – \theta_1\). Let the distance l between the two singularities approach zero while the strength magnitude is increasing so that the product \(\Lambda l\) remains constant. In the limit, this flow pattern is a doublet and we define its strength by \(k=l \Lambda\). The stream function for a doublet is obtained from previous equation as follows:

\(\psi \left( {x,y} \right) = \mathop {\lim }\limits_{l \to 0} \left( { – \frac{\Lambda }{{2\pi }}d\theta } \right){\text{ and }}k = l\Lambda = {\text{cost}}\)

where in the limit \(\Delta \theta \to d\theta \to 0\).

Considering the case where dθ is infinitesimal, we deduce from the figure above that:

Hence the stream function becomes:

i.e.

\(\psi(r, \theta)=-\frac{\kappa}{2 \pi} \frac{\sin \theta}{r}\)

In Cartesian coordinates, a doublet located at the origin has the stream function

\(\psi(x, y)=-\frac{\kappa}{2 \pi} \frac{y}{x^{2}+y^{2}}\)

from which we can derive the velocity components

Now we have done the math, it is time to code and visualize what the streamlines look like. We start by creating a mesh grid.

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 = numpy.linspace(x_start, x_end, N)    # creates a 1D-array for x
y = numpy.linspace(y_start, y_end, N)    # creates a 1D-array for y
X, Y = numpy.meshgrid(x, y)              # generates a mesh grid

We consider a doublet of strength κ=1.0 located at the origin.

kappa = 1.0                        # strength of the doublet
x_doublet, y_doublet = 0.0, 0.0    # location of the doublet

We play smart by defining functions to calculate the stream function and the velocity components that could be re-used if we decide to insert more than one doublet in our domain.

Once the functions have been defined, we call them using the parameters of the doublet: its strength kappa and its location x_doublet, y_doublet.

# compute the velocity field on the mesh grid
u_doublet, v_doublet = get_velocity_doublet(kappa, x_doublet, y_doublet, X, Y)

# compute the stream-function on the mesh grid
psi_doublet = get_stream_function_doublet(kappa, x_doublet, y_doublet, X, Y)

We are ready to do a nice visualization.

# 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_doublet, v_doublet,
                  density=2, linewidth=1, arrowsize=1, arrowstyle='->')
pyplot.scatter(x_doublet, y_doublet, color='#CD2305', s=80, marker='o');

pyplot.show()

The combinatorial calculus studies the groupings that can be obtained with a given number n of objects arranged on a given number k of places.

The groupings can be formed without repetitions or with repetitions of the n objects.

According to Wikipedia: Factorial of a positive integer n, denoted by n! is the product of all positive integers less than or equal to n.

You can calculate factorials with the following formula:

\(n! = n \times \left( {n – 1} \right) \times \ldots \times 2 \times 1 = \prod\limits_{i = 1}^n i \)

Now you might be wondering how you would go about calculating factorials in Python. While I’m certain in the existence of out-of-the-box functions in various libraries, it’s really easy to define your own function, and that’s exactly what we’ll do.

Here’s a simple recursive function which will get the job done:

def factorial(n):
    if n == 1: return 1
    else: return n * factorial(n-1)

In mathematics, a permutation of a set is, loosely speaking, an arrangement of its members into a sequence or linear order, or if the set is already ordered, a rearrangement of its elements. The word “permutation” also refers to the act or process of changing the linear order of an ordered set. In other words, they are the groupings made when the number of objects equals the number of places and counts the order in which they are arranged. Permutations can be without repetitions of objects or with repetition of objects. In particular, we will talk about permutations if the number of objects corresponds to the number of places; if the number of objects is different from the number of places, we will speak of dispositions instead.

A combination is a selection of items from a collection, such that (unlike permutations) the order of selection does not matter.

Let’s work this out with an example.

You have a website on which users can register. They need to provide a password that needs to be exactly 8 characters long, and characters cannot repeat. We first need to determine how many characters and digits there are in the English alphabet:

• the number of letters: 26
• the number of digits: 10

Which is 36 in total. So n = 36. r would then be 8, because the password needs to be 8 characters long. Once we know that, it’s easy to calculate the number of unique passwords, given the following formula:

\(D_{n,k}=\frac{{n!}}{{\left( {n – r} \right)!}}\)

If you went ahead and calculated by hand:

\(\frac{{36!}}{{\left( {36 – 8} \right)!}} = 1220096908800\)

Or even in Python, it really is a trivial task:

def factorial(n):
    if n == 1: return 1
    else: return n * factorial(n-1)

def disposition_without_repetition(n,r):
    return (factorial(n) / factorial(n-r))


m= disposition_without_repetition(36, 8)
print(m)

Okay, cool, but I want to allow my users to repeat characters. No problem, in that case, we’re talking about dispositions with repetition, and the formula is even simpler:

\(D_{n,k}^r = {n^k}\)

You already know what n is 36, and what r is 8, so here’s the solution:

def factorial(n):
    if n == 1: return 1
    else: return n * factorial(n-1)

def disposition_without_repetition(n,r):
    return (factorial(n) / factorial(n-r))

def disposition_with_repetition(n,r):
    return n ** r


m= disposition_with_repetition(36, 8)
print(m)

or:

\(D_{36,8}^r=36^8=2821109907456\)

That’s a whole lot of password options.

How many anagrams, even without meaning, can be formed with the word “MAMMA”?

In this case, the order counts, but the five letters are not all distinct: M is repeated 3 times; A is repeated 2 times.

\(P_n^r = \frac{{5!}}{{3! \cdot {2_2}!}} = \frac{{120}}{{6 \cdot 2}} = 10\)

or, in different words, there are 10 words that can be formed with the letters of the word MAMMA.

Now, consider the following sentence: A group of people selected for a team is the same group, the order doesn’t matter. That’s the whole idea behind combinations. If you select 5 members for the team, you could order them by name, height, or something else, but essentially you would still have the same team — ordering is irrelevant.

So, let’s formalize this idea with a formula. The number of combinations C of a set of n objects taken r at a time is calculated as follows:

\(\boldsymbol{C}(\boldsymbol{n}, \boldsymbol{r})=\frac{\boldsymbol{n} !}{\boldsymbol{r} !(\boldsymbol{n}-\boldsymbol{r}) !}\)

Now you can take that equation and solve the following task: On how many ways can you choose 5 people from a group of 10 for a football team?

The group would be the same, no matter the ordering. So let’s see, n would be equal to 10, and r would be 5:

\(C\left( {10,5} \right) = \frac{{10!}}{{5!\left( {10 – 5} \right)!}} = 252\)

This can once again be easily done with Python:

def factorial(n):
    if n == 1: return 1
    else: return n * factorial(n-1)

def disposition_without_repetition(n,r):
    return (factorial(n) / factorial(n-r))

def disposition_with_repetition(n,r):
    return n ** r

def combinations_without_repetition(n,r):
    return (factorial(n)/(factorial(r) * (factorial(n-r))))



m= combinations_without_repetition(36, 8)
print(m)

Great! But now you might be wondering if there exists a version of combinations which allows repetition. The answer is yes. I’ll explain now.

Assigned two droppers, the first contains 5 drops of white color and the second 5 drops of black color. By mixing 5 drops chosen between the two, how many different colors can be formed?

Answer:

\(C_{2,5}^r = \frac{{\left( {2 + 5 – 1} \right)!}}{{5!\left( {2 – 1} \right)!}} = \frac{{6!}}{{5! \cdot 1!}} = 6\)

def factorial(n):
    if n == 1: return 1
    else: return n * factorial(n-1)

def disposition_without_repetition(n,r):
    return (factorial(n) / factorial(n-r))

def disposition_with_repetition(n,r):
    return n ** r

def combinations_without_repetition(n,r):
    return (factorial(n)/(factorial(r) * (factorial(n-r))))

def combinations_with_repetition(n,r):
    return ((factorial(n+r-1))/(factorial(r)*(factorial(n-1))))


m= combinations_with_repetition(2, 5)
print(m)

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

The definition of circulation:

In words, the circulation is the line integral of velocity around a closed contour.

In accord with Stockes’ teorem, this line integral is equal to the flux through the contour of the curl of velocity, which is the vorticity, ω=∇×v:

If the vorticity is zero (irrotational flow), so is the circulation around any closed contour equal to zero. This means that the line integral of velocity for any curve going from A to B must be equal and opposite to that of any curve going back from B to A. Expand the dot product in the integral, where the velocity is v=(u,v,w) :

In irrotational flow, it doesn’t matter what path you take, this line integral from A to B is always the same value.

Remeber that u dx+v dy+w dz  is an exact differential of a potential ϕ, where:

Or, for short: v=∇ϕ. Applying the continuity equation for incompressible flow, ∇⋅v=0, we get the beautifully simple governing equation of potential flow:

Laplace’s equation! So any solution to Laplace can be a potential flow.

We want to numerically express the flow field of a source and a sink, two potential flow solutions, so we can plot these flows and admire them.

Let’s start importing some useful Python libraries for our program:

• NumPy is a scientific library to create and manage multi-dimensional arrays and matrices.
• Matplotlib is a 2D plotting library that we will use to visualize our results.
• the math module provides the mathematical functions defined by the C standard.

import math
import numpy as np
from matplotlib import pyplot

The objective is to visualize the streamlines corresponding to a source and a sink. To do that, we need to first define a set of points where the velocity components will be computed.

Let’s define an evenly spaced Cartesian grid of points within a spatial domain that is 4 units of length wide in the x-direction and 2 units of length wide in the y-direction, i.e. x,y∈[−2,2],[−1,1].

The variable N will be the number of points we want in each direction, and we define the computational boundaries by the variables x_start, x_end, y_start and y_end.

N = 50                                	 # number of points in each direction
x_start, x_end = -2.0, 2.0            	 # boundaries in the x-direction
y_start, y_end = -1.0, 1.0            	 # boundaries in the y-direction

x = np.linspace(x_start, x_end, N)    # creates a 1D-array with the x-coordinates
y = np.linspace(y_start, y_end, N)    # creates a 1D-array with the y-coordinates

X, Y = np.meshgrid(x, y)              # generates a mesh grid

The last line of the code block calls the meshgrid() function, which generates arrays containing the coordinates of points where the numerical solution will be calculated.

Now that the mesh grid has been generated, it is time to visualize it with the module pyplot from the library matplotlib using the function scatter().

# plot the grid of points
width = 10.0
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.scatter(X, Y, s=5, color='#CD2305', marker='x')

On all of those nicely ordered points, we now will calculate the velocity vector corresponding to a source flow. Then we’ll plot the streamlines.

Source flow

A very important quality of the potential flow is that the governing equation is linear and the solutions can be constructed by superposition. For this reason it is very useful to have a toolbox of elementary solutions that we can use as building blocks. Sources and sinks are such elementary solutions.

A source is a point from which we imagine that fluid is flowing out, uniformly. Thus, all the streamlines radiate from a single point as straight lines and the radial velocity decreases with the distance from the source point.

Let’s consider first the purely two-dimensional case. Because of the radial symmetry, it is convenient to use a cylindrical coordinate system, (r,θ). The angle θ is tan⁻¹(y/x). The velocity components (radial and tangential) are:

\(u_{r}(r, \theta)=\frac{\Lambda}{2 \pi r}, \quad u_{\theta}(r, \theta)=0\)

where \(\Lambda\) represents the source strength.

For the stream function we have:

\(\psi = \frac{\Lambda}{2 \pi} \theta+\text{cost}\)

In practical problems, we are more interested in the velocity components that are obtained by differentiation of the stream function, so that the constant can be dropped.

In Cartesian coordinates, the velocity field (u,v) at position (x,y) corresponding to a source of strength \(\Lambda\) located at (x_source,y_source) is given by:

\(u=\frac{\partial \psi}{\partial y}=\frac{\Lambda}{2 \pi} \frac{x-x_{\text {source }}}{\left(x-x_{\text {source }}\right)^{2}+\left(y-y_{\text {source }}\right)^{2}}\)

and

\(v=-\frac{\partial \psi}{\partial x}=\frac{\Lambda}{2 \pi} \frac{y-y_{\text {source }}}{\left(x-x_{\text {source }}\right)^{2}+\left(y-y_{\text {source }}\right)^{2}}\)

Let’s calculate the velocity field for our grid of points. We’ll place the source at the location (−1,0) and give it a strength \(\Lambda=5\).

Instead of picking one point on the grid and calculate its velocity (which means that we would have to iterate over all positions [i,j]), we directly compute velocity arrays (u_source, v_source) using arithmetic operators on arrays. Yes, with Numpy, arithmetic operators on array apply elementwise and a new array is created and filled with the result.

Lambda= 5.0                                # source strength
x_source, y_source = -1.0, 0.0             # location of the source

# compute the velocity field on the mesh grid
u_source = (Lambda / (2 * math.pi) *
            (X - x_source) / ((X - x_source)**2 + (Y - y_source)**2))
v_source = (Lambda / (2 * math.pi) *
            (Y - y_source) / ((X - x_source)**2 + (Y - y_source)**2))

Let’s plot the stream lines already:

# plot the streamlines
width = 10.0
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_source, v_source,
                  density=2, linewidth=1, arrowsize=2, arrowstyle='->')
pyplot.scatter(x_source, y_source,
               color='#CD2305', s=80, marker='o');

pyplot.show()

Sink flow

In the source flow, the strength \(\Lambda\) was chosen to be positive. A source with a negative strength is called a sink. Instead of radiating from a single point, the straight streamlines are now converging to a single point.

The velocity field corresponding to a sink looks similar to that of a source, except for the direction of the flow. Thus, the Python code requires very few modifications.

We will place the sink at the location (1,0) and give it an equal strength to our source, but negative of course.

Lambda_s = -5.0                          # strength of the sink
x_sink, y_sink = 1.0, 0.0                # location of the sink

# compute the velocity on the mesh grid
u_sink = (Lambda_s / (2 * math.pi) *
          (X - x_sink) / ((X - x_sink)**2 + (Y - y_sink)**2))
v_sink = (Lambda_s / (2 * math.pi) *
          (Y - y_sink) / ((X - x_sink)**2 + (Y - y_sink)**2))

and then,

# plot the streamlines
width = 10.0
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_sink, v_sink,
                  density=2, linewidth=1, arrowsize=2, arrowstyle='->')
pyplot.scatter(x_sink, y_sink,
               color='#CD2305', s=80, marker='o');

pyplot.show()

Source-sink pair

Now, let’s to see the superposition’s powers. We already have the velocity field of the source and the velocity field of the sink. We can just add these velocity fields, point wise, to get a new solution of potential flow: the source-sink pair.

# compute the velocity of the pair source/sink by superposition
u_pair = u_source + u_sink
v_pair = v_source + v_sink

# plot the streamlines of the pair source/sink
width = 10.0
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_pair, v_pair,
                  density=2.0, linewidth=1, arrowsize=2, arrowstyle='->')
pyplot.scatter([x_source, x_sink], [y_source, y_sink], 
               color='#CD2305', s=80, marker='o');

pyplot.show()