Source & Sink

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