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*_{2},0) and a sink of opposite strength located at (*l*_{2},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()
```