How can one produce a proper streamline plot?

Question

I recently had trouble producing a proper streamline plot in Mathematica, and apparently the problem is a good bit harder than I appreciated. I would like to know if there exist general algorithms to produce streamline plots which deal appropriately with line spacings and terminations.

Let me make the problem specific. Consider a vector field $\vec F(\vec x)$, in two or in three dimensions. A streamline is an integral curve of the vector field, that is, a curve $\vec\alpha(s)$ whose derivative $\frac{d\vec\alpha}{ds}$ is proportional at every point to the vector field $\vec F(\vec\alpha(s))$ at that point. A streamline plot represents a collection of such streamlines, and it is an extremely useful tool for visualizing the field, and its direction in particular.

There are, however, additional ways to make the diagrams convey information about the field's magnitude and its divergence.

• If $\vec F$ is divergence-free (i.e. $\nabla\cdot\vec F=0$), then it is customary to draw streamlines that do not terminate unless they meet a zero or singularity of the field or go out of the plot, and to space the streamlines such that the vector field flow $$ \int_L\vec F\cdot \text d\vec x$$ is constant for $L$ connecting any two adjacent streamlines. (In three dimensions this is harder to define appropriately, but the flow $ \int_S\vec F\cdot \text d\vec a$ should be constant over any 'unit cell' whose corners are streamlines.)

  This allows for a direct reading of the intensity of the vector field from the line spacing, which is guaranteed to be smaller where the field is more intense.

• If $\vec F$ has a nonzero divergence, then it is still customary to keep the line spacing such that the flow between streamlines is approximately constant, though this requires that lines terminate or initiate in the middle of the diagram.

  This allows one to visualize regions with positive or negative divergence - volume sources and sinks of the field - as those in which streamlines initiate or terminate, respectively.

So: I would like to know if there exist algorithms that will take a vector field $\vec F$ and a region $R$, in two or three dimensions (mostly two), and which will produce a set of streamline termini such that the field flow in between streamlines is (approximately) constant? It should therefore have lines initiate and terminate at and only at regions with positive and negative field divergence.

Answer 1

You can calculate the stream function yourself, and from it you can draw contours or streamlines of constant $\psi$. Let us assume you are on a two-dimensional incompressible flow $\mathbf{u}=(u,v,0)^T$, then you can find an exact differential $d\psi$ that satisfies the mass conservation equation:

$ \frac{\partial u}{\partial x } + \frac{\partial v}{\partial y} = 0$

whence:

$d\psi = udy - vdx$

where:

$u = \frac{\partial \psi}{\partial y }, \ v = -\frac{\partial \psi}{\partial x }$

Over an arbitrary path immersed in the fluid, you can integrate the exact differential equation:

$\psi - \psi_0 = \int_{path} udy - vdx $

From the implementation perspective, if you are working on structured orthogonal grids you can implement a Simpson's rule to integrate the above equation in a given direction (have in mind that if you choose to integrate along one of the coordinates, one of the terms inside the integral vanishes identically).

For three dimensional flows, one can define two stream-functions instead of one. By assuming the velocity in the continuity equation ($\nabla \cdot \mathbf{u} = 0 $) as a vector vector potential $\mathbf{u}=\nabla \times \mathbf{B}$, where $\mathbf{B} = \phi \nabla \psi$. Thus, one can extend the above definition by noting:

$\mathbf{u} = \nabla \phi \times \nabla \psi $

The bulk flow over a volume enclosed by the surfaces of constant $\phi$ and $\psi$ can be calculated, using the Gauss Theorem, as :

$ M = \int \int_{\partial\Omega} \mathbf{u} \, d\mathbf{\nu} = \oint \phi \nabla \psi \, dl$

Finally, by taking the curl (vorticity) of the velocity potential $\mathbf{B}$ and after a bit of vector algebra you get:

$\mathcal{L}(\phi)\nabla \psi - \mathcal{L}(\psi)\nabla \phi = \mathbf{\omega}$

Where, $\mathcal{L}(\phi) = -\mathbf{I}\nabla\cdot (\nabla \phi) + \underline {\nabla}(\nabla \phi) $

Note: The underline indicates that the operator is a column vector.

The last equation can be reduced to simpler cases (irrotational flow, axisymmetric). You can extend the method further, by decomposing your vector field in a set of a curl-free and a divergence-free field by means of the Hemholtz's Theorem.

Answer 2

I have used matplotlib's plot_streamlines function in the past to create streamline plots. There is a nice introduction to the functionality here. If you want more control over the creation of the plot this thread at StackOverflow offers more flexible techniques (and source code) for creating streamline plots in Python.