3D Ferrofluid Rosensweig Instability: Spike Pattern on a Coil with Gauss Meter

Setting up the bench...

Presets

Live readouts

Graphs

Dispersion relation omega^2(k) for B = 0, B = Bc and the present field (negative = unstable)

Neutral stability curve: the field needed to destabilise each wavenumber (minimum = Bc at kc)

Measured power spectrum of the surface height from the running simulation

RMS surface height against simulated time (log scale); green: saturation height scale

Experiments

Ramp: watch the surface stay flat until B = Bc, then the spikes grow. Compare the measured wavelength with the critical one.

Numerical checks

running...

How to use

The preset Typical ferrofluid, ramp up raises the field slowly: the surface is flat, then at B = Bc the first peaks appear with spacing close to the critical wavelength (green bar). Change density, surface tension, permeability or gravity and see the threshold and wavelength move on the graphs. Drag to orbit and wheel to zoom; the gauss meter and the power-supply knob follow the field. Space pauses, R restarts.

The dish, coil, instruments and bench were built in Blender; all dynamics run in JavaScript.

The science

The phenomenon. A ferrofluid layer in a vertical magnetic field is flat until the field reaches a critical value, then the surface suddenly breaks into a pattern of peaks (the Rosensweig or normal-field instability). Gravity and surface tension flatten the surface; the magnetic normal traction pulls any bump upward, and above a threshold the pull wins.

Linear theory. For a thick layer with linear magnetisation (relative permeability mu_r) and a sharp interface, a ripple of wavenumber k obeys (Cowley and Rosensweig 1967) rho omega^2 = rho g k + sigma k^3 - (B^2 / mu0) F k^2, F = (mu_r - 1)^2 / (mu_r (mu_r + 1)) where B is the induction normal to the layer. The first two terms are gravity and capillarity; the last one is the magnetic destabilisation. Instability (omega^2 < 0) appears first at the wavenumber that minimises rho g / k + sigma k, which is independent of the field: k_c = sqrt(rho g / sigma), lambda_c = 2 pi sqrt(sigma / (rho g)), B_c^2 = 2 mu0 sqrt(g rho sigma) / F With rho = 1200 kg/m3, sigma = 28 mN/m, mu_r = 3 this gives lambda_c about 10 mm and B_c about 11 mT, the right order of magnitude for commercial ferrofluids. Above B_c a band of wavenumbers around k_c grows exponentially with rate sqrt(-omega^2); stronger fields shift the fastest mode to shorter wavelengths.

The simulation. The surface height eta(x, y, t) lives on a 64 x 64 grid of 1 mm cells. Each Fourier mode obeys the linear equation above plus viscous damping 2 nu k^2 d(eta)/dt, a quadratic term +alpha eta^2 that favours upward peaks (this is what makes hexagonal spot arrays) and a cubic term -beta eta^3 that limits the growth (beta is chosen so that the height saturates near the "saturation height scale"). The mean height is held constant (fluid volume). The display adds a mild sharpening of crests.

Honest simplifications. (1) The linear threshold and wavelength are the physics; the pattern after onset is a generic saturated one, not a solution of the full nonlinear magnetohydrodynamic problem. Real spikes are conical, hexagonally ordered, and the transition is subcritical with hysteresis (spikes persist when the field is lowered below B_c); this model has none of that. (2) Linear magnetisation: real fluids saturate at high field. (3) Infinite fluid depth and a uniform vertical field; a finite layer, the dish wall and field gradients change B_c and the pattern. (4) Damping is 2 nu k^2 plus a bulk friction 0.8 sqrt(g k_c) standing in for boundary-layer dissipation; it shapes the transient (how quickly spikes settle) but does not move the threshold. (5) Time is slowed by the factor shown (real growth takes a fraction of a second).

Checks run in the page. The minimum of the neutral curve equals the closed-form B_c; its position equals k_c; the growth rate of a single mode measured in the solver equals sqrt(-omega^2); the oscillation frequency of a stable mode at B = 0 equals sqrt(g k + sigma k^3 / rho); the 2D FFT round trip is exact to 1e-12.