Solution (source code)

= Solution

For axisymmetric incompressible flow without swirl, the <Stokes streamfunction> satisfies
$$
\boxed{u_R=\frac{\Psi_\theta}{R^2\sin\theta},\qquad
u_\theta=-\frac{\Psi_R}{R\sin\theta}.}
$$
Substitution into the divergence in spherical coordinates makes <incompressibility> an identity. Its only <vorticity> component is
$$
\omega_\phi=\frac1R[\partial_R(Ru_\theta)-\partial_\theta u_R]
=-\frac{\mathcal D^2\Psi}{R\sin\theta}.
$$
Taking the curl of the <Stokes equation> eliminates pressure and gives $\nabla^2\boldsymbol\omega=0$. For an axisymmetric azimuthal vector, its scalar component obeys $(\Delta-1/(R^2\sin^2\theta))\omega_\phi=0$. Direct differentiation verifies
$$
\left(\Delta-\frac1{R^2\sin^2\theta}\right)
\frac{h}{R\sin\theta}
=\frac{\mathcal D^2h}{R\sin\theta}.
$$
Apply this to $h=\mathcal D^2\Psi$ to obtain
$$
\boxed{\mathcal D^2(\mathcal D^2\Psi)=0.}
$$