Divergence, surface integrals, and flux

Flux and surface integrals

Consider a surface \(S\). Let \(\bm{r}(u,v) = [x(u,v), y(u,v),z(u,v)]\) be a parametric position vector on a Euclidean vector space \(\mathbb{R}^3\). Furthermore, let \(\bm{F}: \mathbb{R}^3 \rightarrow \mathbb{R}^3\) be a vector-valued function of \(\bm{r}\) and let \(\bm{n}\) be a unit-normal vector on a surface \(S\). The \[\begin{align} \iint\limits_S \bm{F} \cdot \bm{n} \diff S \end{align}\] which integrates the normal of \(\bm{F}\) over the surface. We call this quantity the of \(\bm{F}\) out of the surface \(S\). This terminology comes from fluid flow, for which the flux is the mass flow rate out of \(S\). In general, the flux is a measure of a quantity (or field) passing through a surface. For more on computing surface integrals, see schey2005 and kreyszig2011.

Continuity

Consider the flux out of a surface \(S\) that encloses a volume \(\Delta V\), divided by that volume: \[\begin{align} \frac{1}{\Delta V} \iint\limits_S \bm{F} \cdot \bm{n} \diff S. \end{align}\] This gives a measure of flux per unit volume for a volume of space. Consider its physical meaning when we interpret this as fluid flow: all fluid that enters the volume is negative flux and all that leaves is positive. If physical conditions are such that we expect no fluid to enter or exit the volume via what is called a or a , and if we assume the density of the fluid is uniform (this is called an fluid), then all the fluid that enters the volume must exit and we get \[\begin{align} \frac{1}{\Delta V} \iint\limits_S \bm{F} \cdot \bm{n} \diff S = 0. \end{align}\]

This is called a , although typically this name is given to equations of the form in the next section. This equation is one of the governing equations in continuum mechanics.

Divergence

This is called the of \(\bm{F}\) and is defined at each point in \(\mathbb{R}^3\) by taking the volume to zero about it. It is given the shorthand \(\Div\bm{F}\).

What interpretation can we give this quantity? It is a measure of the vector field’s flux outward through a surface containing an infinitesimal volume. When we consider a fluid, a positive divergence is a local decrease in density and a negative divergence is a density increase. If the fluid is incompressible and has no sources or sinks, we can write the continuity equation \[\begin{align} \Div\bm{F} = 0. \end{align}\]

Exploring divergence

Divergence is perhaps best explored by considering it for a vector field in \(\mathbb{R}^2\). Such a field \(\bm{F} = F_x \bm{\hat{i}} + F_y \bm{\hat{j}}\) can be represented as a “quiver” plot. If we overlay the quiver plot over a “color density” plot representing \(\Div\bm{F}\), we can increase our intuition about the divergence.

First, load some Python packages.

import numpy as np
import sympy as sp
import matplotlib.pyplot as plt
from matplotlib.ticker import LogLocator
from matplotlib.colors import *
from sympy.utilities.lambdify import lambdify

Now we define some symbolic variables and functions.

x = sp.Symbol('x', real=True)
y = sp.Symbol('y', real=True)
F_x = sp.Function('F_x')(x, y)
F_y = sp.Function('F_y')(x, y)

Rather than repeat code, let’s write a single function quiver_plotter_2D() to make several of these plots.

def quiver_plotter_2D(
        field={},
        grid_width=3, grid_decimate_x=8, grid_decimate_y=8,
        norm=Normalize(), density_operation='div', 
        print_density=True):
    x, y = sp.symbols('x y', real=True)
    F_x, F_y = sp.Function('F_x')(x, y), sp.Function('F_y')(x, y)
    field_sub = field
    # Calculate density
    den = F_x.diff(x) + F_y.diff(y) if density_operation == 'div' else None
    if den is None:
        raise ValueError(f'Unknown density operation: {density_operation}')
    den_simp = den.subs(field_sub).doit().simplify()
    if den_simp.is_constant():
        print('Warning: density operator is constant (no density plot)')
    if print_density:
        print(f'The {density_operation} is:')
        print(den_simp)
    # Lambdify for numerics
    F_x_sub = F_x.subs(field_sub)
    F_y_sub = F_y.subs(field_sub)
    F_x_fun = sp.lambdify((x, y), F_x.subs(field_sub), 'numpy')
    F_y_fun = sp.lambdify((x, y), F_y.subs(field_sub), 'numpy')
    if F_x_sub.is_constant:
        F_x_fun1 = F_x_fun  # Dummy
        F_x_fun = lambda x, y: F_x_fun1(x, y) * np.ones(x.shape)
    if F_y_sub.is_constant:
        F_y_fun1 = F_y_fun  # Dummy
        F_y_fun = lambda x, y: F_y_fun1(x, y) * np.ones(x.shape)
    if not den_simp.is_constant():
        den_fun = sp.lambdify((x, y), den_simp, 'numpy')
    # Create grid
    w = grid_width
    Y, X = np.mgrid[-w:w:100j, -w:w:100j]
    # Evaluate numerically
    F_x_num = F_x_fun(X, Y)
    F_y_num = F_y_fun(X, Y)
    if not den_simp.is_constant():
        den_num = den_fun(X, Y)
    # Plot
    fig, ax = plt.subplots()
    if not den_simp.is_constant():
        cmap = plt.get_cmap('coolwarm')
        im = plt.pcolormesh(X, Y, den_num, cmap=cmap, norm=norm)
        plt.colorbar()
    dx, dy = grid_decimate_y, grid_decimate_x
    plt.quiver(X[::dx, ::dy], Y[::dx, ::dy], F_x_num[::dx, ::dy], 
               F_y_num[::dx, ::dy], units='xy', scale=10)
    plt.title(fr'$F(x, y) = \left[{sp.latex(F_x.subs(field_sub))},' +
              fr'{sp.latex(F_y.subs(field_sub))}\right]$')
    return fig, ax

Let’s inspect several cases.

fig, ax = quiver_plotter_2D(field={F_x: x**2, F_y: y**2})
plt.draw()

fig, ax = quiver_plotter_2D(field={F_x: x*y, F_y: x*y})
plt.draw()

fig, ax = quiver_plotter_2D(field={F_x: x**2 + y**2, F_y: x**2 + y**2})
plt.draw()

fig, ax = quiver_plotter_2D(
  field={F_x: x**2/sp.sqrt(x**2+y**2), F_y: y**2/sp.sqrt(x**2+y**2)},
  norm=SymLogNorm(linthresh=.3, linscale=.3)
)
plt.show()