Vortex Loops for Jupiters’ Stripes

Project Members: Aleksa Milovanović, Reid Tang

Mentor: Sadashige Ishida

Introduction to Fluid Dynamics

Fluid motion in computer graphics is typically modeled by the incompressible Navier-Stokes equations:

\frac{\partial \vec{u}}{\partial t} + \vec{u} \cdot \nabla \vec{u} + \frac{1}{\rho} \nabla p = \vec{g} + \nu \nabla \cdot \nabla \vec{u},

\nabla\cdot\vec{u} = 0,

where is the velocity of the fluid, density, pressure, body force and kinematic viscosity. The first equation (momentum) is Newton’s second law applied to fluid elements, and the second one constrains the flow to be volume-preserving.

An important concept is the material derivative, which measures the rate of change of a quantity following a fluid parcel (Lagrangian) rather than at a fixed point in space (Eulerian).

This allows us to write the momentum equation compactly as

\frac{D\vec{u}}{Dt} + \frac{1}{\rho}\nabla p = \vec{g} + \nu\nabla\cdot\nabla\vec{u},

and defines advection s.t. means the quantity is carried by the flow without changing.

Vorticity

Vorticity is defined as the curl of velocity , and measures the local rate of rotation of the fluid. For inviscid flow with no body force, taking the curl of the momentum equation gives the vorticity transport equation

\frac{D\vec{\omega}}{Dt} = (\vec{\omega}\cdot\nabla)\vec{u}

where the right-hand side is the vortex-stretching term (as shown in Section 2.1, this term vanishes identically in the 2D setting used throughout this report, so vorticity there is materially conserved exactly, ); it is simply carried by the flow. This makes vorticity important for simulating turbulent and swirl-dominated phenomena, since instead of solving for pressure and velocity everywhere, one only tracks regions where vorticity is concentrated.

The key problem is recovering the velocity field given a vorticity field. In 3D, this is posed as a least-squares problem, i.e. find minimizing

\iiint \|\nabla\times\vec{u} – \vec{\omega}\|^2\,dV

solved via calculus of variations, yielding the Poisson problem , whose solution is the Biot-Savart law:

\vec{u}(\vec{x}) = \iiint -\nabla\Phi_3(\vec{x} – \vec{p})\times\vec{\omega}(\vec{p})\,d\vec{p}, \quad \Phi_3(\vec{x}) = -\frac{1}{4\pi\|\vec{x}\|}

Deriving the 2D analogue

Jupiter’s bands and jet streams are a 2D phenomenon (a thin shell on the planet’s surface), so we re-derive the above for a 2D velocity field , where vorticity reduces to a scalar:

\omega = \frac{\partial v}{\partial x} – \frac{\partial u}{\partial y} =: \mathrm{curl}_2(\vec{u})

In this 2D setting points purely out of the plane while is purely in-plane, so the vortex-stretching term from Section 2 vanishes identically — not as an extra assumption, but automatically, because there is no out-of-plane direction left for vorticity to be stretched or tilted into. The vorticity transport equation therefore simplifies exactly to in 2D, which is what we use throughout.

As in 3D, we seek minimizing

\iint \big(\mathrm{curl}_2(\vec{u}) – \omega\big)^2\,dA

Perturbing and using linearity of , the stationarity condition gives, with ,

\iint f\, \mathrm{curl}_2(\vec{p})\,dA = 0 \ \text{, for all } \vec{p}

Unlike 3D, maps a vector field to a scalar, so its adjoint maps a scalar back to a vector. Integrating by parts term by term (boundary terms vanish as fields decay at infinity),

\iint f\,\mathrm{curl}_2(\vec{p})\,dA = \iint \vec{p}\cdot\nabla^\perp f\,dA, \quad \nabla^\perp f \equiv \left(\frac{\partial f}{\partial y},-\frac{\partial f}{\partial x}\right)

so that, since this must vanish for all ,

\nabla^\perp\big(\mathrm{curl}_2(\vec{u})\big) = \nabla^\perp\omega

The same vector identity used in 3D,

\nabla\times\nabla\times\vec{u} = -\nabla\cdot\nabla\vec{u} + \nabla(\nabla\cdot\vec{u})

restricted to 2D becomes . With the constraint , we obtain a Poisson problem:

-\nabla\cdot\nabla\vec{u} = \nabla^\perp\omega

The 2D fundamental solution of the Laplacian is

\Phi_2(\vec{x}) = \frac{1}{2\pi}\log\|\vec{x}\|

satisfying , verified by the divergence theorem on the unit circle:

\oint \nabla\Phi_2\cdot\hat{n}\,d\ell = \oint \frac{1}{2\pi}\,d\ell = 1

The solution to the Poisson problem, via convolution with :

\vec{u}(\vec{x}) = \iint \Phi_2(\vec{x} – \vec{p})\,\nabla^\perp\omega(\vec{p})\,d\vec{p}

Integrating by parts once again (moving off onto ) gives the 2D Biot-Savart law in the form:

\vec{u}(\vec{x}) = \iint \omega(\vec{p})\,\nabla^\perp\Phi_2(\vec{x}-\vec{p})\,d\vec{p}

Discretizing as a sum of Dirac deltas collapses the integral to a sum:

\vec{u}(\vec{x}) = \sum_i \frac{\Gamma_i}{2\pi} \frac{(\vec{x} – \vec{x}_i)^\perp}{\|\vec{x} – \vec{x}_i\|^2}

Since inviscid vorticity is materially conserved, each point vortex simply moves with the velocity induced by all other vortices at its own location (the self-term is excluded), giving the point vortex dynamics equation:

\frac{d\vec{x}_i}{dt} = \sum_{j\neq i}\frac{\Gamma_j}{2\pi\|\vec{x}_i-\vec{x}_j\|^2}\big(-(y_i-y_j),\,x_i-x_j\big)

This is the equation we implemented and verified on various examples (e.g. leapfrogging etc.) before extending it to the sphere.

Implementation and Validation in the Plane

The point vortex dynamics equation above was implemented in Houdini using a custom SOP Solver: each timestep, a VEX wrangle reads the previous frame’s point cloud (position and circulation per point), evaluates the Biot–Savart sum, and advances each point’s position. Before trusting this on any nontrivial configuration, it was validated against two exactly-solvable two-vortex cases: equal circulations , where both vortices orbit a common, stationary centroid on a circular path, and opposite circulations , where the pair instead translates together in a straight line at constant speed. Both matched their analytic period/speed formulas.

As a further, qualitative check, a classical four-vortex leapfrogging configuration (two counter-rotating pairs, arranged so the trailing pair periodically overtakes the leading pair) was simulated and reproduced the expected periodic behavior.

Trajectories of the four-vortex leapfrogging configuration at three points in the cycle (frames 20, 80, 140). Each mark traces the recent path of one vortex; the pairs periodically exchange places as they orbit and overtake one another.

A Vortex-Ring Model of Kelvin–Helmholtz Instability

To move from the plane to Jupiter’s spherical atmosphere, the point vortex equation was extended to the unit sphere by replacing the planar kernel with the spherical Biot–Savart law and integrating positions via the exponential map, , so that points remain exactly on the sphere regardless of step size. Because has no boundary, the total vorticity must also integrate to zero, — a constraint with no analogue in the plane, enforced here by construction (equal numbers of positive and negative vortices, or an explicit normalization step subtracting the mean circulation from every point).

On this spherical setup, a shear layer was modeled as a vortex sheet: a ring of point vortices along a line of latitude, one sign on each side of the interface, with a small sinusoidal perturbation added to seed the instability. With no additional forcing, mutual induction alone causes the shear layer to roll up into a periodic array of vortices, as expected from classical vortex-sheet theory.

A single vortex-sheet interface at three stages: the initial perturbed line (frame 50), partial roll-up into discrete billows (frame 250), and a more developed, mixed state (frame 450).

This was then generalized to several alternating-sign rings restricted to a narrower latitude band, approximating Jupiter’s alternating belt/zone structure. Passive marker particles (zero circulation, contributing nothing to the velocity field but advected by it) were added and colored with a Jupiter-like palette for visualization.

Multi-band configuration (M=3) at frames 20, 250, and 500. Marker particles, seeded uniformly, are progressively advected and mixed along the band-induced shear flow.

Modeling the Coriolis Effect

To incorporate planetary rotation, the momentum equation was written in a frame rotating with angular velocity . Taking the curl and combining terms gives the vorticity equation with a rotational source,

\frac{D\omega}{Dt} = -2\vec{u}\cdot\vec{\Omega},

in place of the inviscid used above. This was implemented by updating each vortex’s circulation every step, , so that circulation is no longer materially conserved. This is a first, direct discretization of the effect; a more careful treatment (in particular, its interaction with the constraint over long integration times) is left for future work.

Discussion and Limitations

The point-vortex approach reproduces known analytic solutions in the plane, the classical leapfrogging trajectory, and qualitatively correct Kelvin–Helmholtz roll-up on the sphere, both for a single interface and for a multi-band configuration resembling Jupiter’s belts and zones. The Coriolis implementation here is a direct, unrefined discretization rather than a fully validated one, and the multi-band results are qualitative: parameters such as band spacing, perturbation amplitude, and desingularization radius were chosen for a visually clear roll-up rather than calibrated against a specific analytic growth rate or against Jupiter’s actual physical parameters. Extending the integrator beyond forward Euler, and validating the Coriolis term against a known conserved quantity, are the natural next steps.