🔬 physics

Development of a Quasi-Static MHD Solver for Turbulent Flows with Non-Uniform Applied Magnetic Fields in OpenFOAM

To address severe discrete cancellation errors and solution divergence in non-uniform magnetic fields that plague standard low-magnetic-Reynolds-number MHD solvers, this paper introduces **mhdturbFoamQS**, a new OpenFOAM solver that eliminates these residuals by directly computing the Lorentz force from the electric-potential formulation and implicitly discretizing magnetic damping, thereby achieving high accuracy against analytical benchmarks without modifying standard turbulence closures.

Original authors: Wellington Fonseca

Published 2026-10-06

📖 4 min read☕ Coffee break read

CC BY 4.0

Original authors: Wellington Fonseca

Original paper licensed under CC BY 4.0 (https://creativecommons.org/licenses/by/4.0/). ✨ This is an AI-generated explanation of the paper below. It is not written or endorsed by the authors. For technical accuracy, refer to the original paper. Read full disclaimer

Imagine a world where invisible magnetic fields can steer the flow of liquid metal, much like a conductor guiding a river. This is the realm of magnetohydrodynamics, a field of physics that studies how electrically conducting fluids, such as liquid metals or salty water, behave when subjected to magnetic forces. These interactions are vital for technologies ranging from cooling systems in nuclear fusion reactors to methods for stirring molten steel without touching it. In many of these applications, the fluid moves through a magnetic field that is not perfectly even; instead, the field strength and direction shift from one spot to another, creating a complex, uneven landscape for the fluid to navigate. To predict how the fluid will move in these conditions, scientists rely on powerful computer simulations. However, a hidden flaw in the standard way these simulations are built has made them unreliable whenever the magnetic field changes across the space, leading to results that are not just slightly off, but completely wrong.

A researcher at the Federal University of Pará in Brazil has identified and fixed this specific problem by developing a new computer program designed to handle these uneven magnetic fields correctly. The issue arose because the standard method used to calculate the magnetic force on the fluid involves subtracting two very large numbers from each other. When the magnetic field is perfectly uniform, these numbers cancel out almost perfectly, leaving the correct, small force. But when the field varies from cell to cell in the computer model, the cancellation is no longer exact. The tiny leftover amount, which should be zero, acts like a phantom force that grows stronger as the magnetic field gets stronger. In tests with a strong ten-tesla magnetic field, this phantom force was found to be nearly one hundred times larger than the actual physical force the fluid should feel, causing the simulation to crash almost immediately.

To solve this, the researcher, Wellington da Silva Fonseca, created a new solver called mhdturbFoamQS. Instead of trying to subtract two large numbers to find the force, this new approach calculates the electric current flowing through the fluid first, and then uses that current to determine the force directly. This method avoids the dangerous subtraction entirely, ensuring that no phantom forces are created, even when the magnetic field is highly irregular. The new program also handles the boundaries of the container more carefully, ensuring that the electric current does not leak out through insulating walls, a detail that previous attempts often got wrong. Furthermore, it treats the magnetic damping effect—the way the magnetic field slows the fluid down—inside the calculation itself, allowing the simulation to run stably without needing to take impossibly small time steps.

The new solver was put through a series of rigorous tests to prove it works. First, it was compared against known mathematical solutions for fluid flowing through a flat channel and a circular pipe. In these tests, the computer results matched the exact theoretical answers with an error of less than half a percent, confirming that the underlying physics was being solved correctly. A particularly revealing test involved a fluid at rest in a magnetic field that had no electric currents flowing through it. Physics dictates that such a fluid should not move at all. The standard simulation method caused the fluid to accelerate to a significant speed due to the calculation error, while the new solver correctly kept the fluid perfectly still, returning a result of zero velocity.

Finally, the researcher used the new tool to simulate a turbulent flow of liquid metal inside a duct surrounded by a ring of magnets that created a six-fold alternating magnetic field. This is the exact type of complex, non-uniform environment where the old methods fail. The simulation ran successfully, showing how the magnetic field flattened the flow near the walls and accelerated the core of the fluid. The study also explored how different mathematical models for turbulence responded to the magnetic field, finding that they predicted different levels of pressure rise and energy, highlighting that while the magnetic force is now calculated correctly, the way turbulence interacts with it remains a complex area for future study. By fixing the fundamental calculation error, this work provides a reliable foundation for engineers and scientists to design and analyze systems involving liquid metals and strong, uneven magnetic fields, ensuring that their computer models reflect reality rather than numerical artifacts.

Technical Summary: Development of a Quasi-Static MHD Solver for Turbulent Flows with Non-Uniform Applied Magnetic Fields in OpenFOAM

Problem Statement
Standard OpenFOAM magnetohydrodynamics (MHD) solvers typically utilize a magnetic-induction formulation where the Lorentz force is evaluated via the Maxwell-stress tensor. While robust for uniform applied magnetic fields, this approach suffers from a severe discrete cancellation error when the applied field ( $B_0$ ) is non-uniform. In the continuous limit, the divergence of the Maxwell stress tensor and the gradient of magnetic pressure cancel exactly. However, in a discrete finite-volume setting with spatially varying fields, this cancellation is imperfect. The residual acts as a spurious body force that scales with the square of the applied field ( $B_0^2$ ). For realistic engineering configurations (e.g., a 10 T field in a 0.15 m duct), this residual can reach accelerations of $\sim 10^3$ m s $^{-2}$ , two orders of magnitude larger than the physical Lorentz acceleration. This error causes solution divergence within a few time steps and is not mitigated by mesh refinement or time-step reduction, indicating a fundamental flaw in the formulation rather than a numerical resolution issue.

Methodology
To address this, the author developed mhdturbFoamQS, a finite-volume solver based on the quasi-static (low magnetic Reynolds number, $R_m \ll 1$ ) electric-potential formulation. The core methodology involves:

Formulation Shift: Instead of calculating the Lorentz force from the Maxwell stress tensor, the solver solves a Poisson equation for the electric potential ( $\phi$ ) derived from charge conservation ( $\nabla \cdot j = 0$ ). The current density ( $j$ ) is reconstructed algebraically from Ohm's law ( $j = \sigma(-\nabla\phi + u \times B_0)$ ), and the Lorentz force is evaluated directly as $j \times B_0$ . This avoids the subtraction of large $O(B_0^2)$ terms.
Generalized Field Handling: The applied field $B_0$ is treated as a general volVectorField rather than a constant, allowing for spatial variations in magnitude and direction.
Consistent Boundary Conditions: A specific implementation of the insulating-wall condition ( $j \cdot n = 0$ ) is applied to both sides of the discrete potential equation. This ensures the motional electromotive force flux is correctly removed from the right-hand side, preventing the inadvertent imposition of an electrode condition.
Implicit Magnetic Damping: The stiff linear damping term ( $-\sigma|B_0|^2 u$ ) arising from the Lorentz force decomposition is treated implicitly within the momentum matrix using fvm::Sp. This removes the restrictive time-step constraint associated with explicit treatment of magnetic damping, which would otherwise scale as $1/B_0^2$ .
Turbulence Coupling: The solver integrates with standard OpenFOAM turbulence closures (Spalart–Allmaras, $k$ – $\epsilon$ , $k$ – $\omega$ SST) without modifying their transport equations. This isolates the electromagnetic formulation from turbulence model changes.

Key Contributions

Identification and Quantification: The paper identifies and quantifies the discrete cancellation error in induction-based solvers under non-uniform fields, demonstrating that the error saturates rather than decays with mesh refinement.
Solver Development: The implementation of mhdturbFoamQS (344 lines of C++), which incorporates the quasi-static potential route, consistent insulating-wall treatment, and implicit magnetic sink handling.
Verification Framework: A comprehensive verification suite including:
- Comparison with the exact Hartmann solution for plane channels (friction factor error < 0.17%, velocity profile error < 0.5%).
- Comparison with the Shercliff–Gold solution for circular pipes with insulating walls.
- A "null test" using a current-free, non-uniform applied field where the fluid is initially at rest. The quasi-static solver correctly returns zero velocity, whereas the induction solver generates spurious motion.
- Manufactured solution tests confirming second-order accuracy for both the potential and reconstructed current.
Reproducibility: The solver, verification cases, and post-processing scripts are released under GPL v3 with a concept DOI on Zenodo, ensuring all reported numbers are independently reproducible.

Results

Verification: The solver reproduces analytical solutions for laminar flows with high precision. In the null test, the induction solver produced spurious velocities up to 0.36 m s $^{-1}$ (and 3.7 m s $^{-1}$ with exact wall values) in a field where the fluid should remain stationary, while the quasi-static solver maintained machine-precision zero.
Turbulent Demonstration: The solver was applied to a turbulent duct flow ( $Re \approx 10^4$ $R e \approx 1 0^{4}$ ) with a realistic, non-uniform alternating-magnet field (six-fold azimuthal structure). Simulations were run with three different turbulence closures at field strengths up to 10 T.
- The solver remained stable and did not exhibit the divergence seen in induction-based approaches.
- The results showed a monotonic increase in volume-averaged pressure and turbulent kinetic energy with field strength.
- Significant sensitivity to the choice of turbulence closure was observed (e.g., at 10 T, pressure rise varied from +17.3% for $k$ – $\epsilon$ to +30.1% for Spalart–Allmaras), highlighting that the electromagnetic formulation is decoupled from the turbulence modeling uncertainty.

Significance and Claims
The paper claims that the quasi-static electric-potential formulation provides a practical and robust method for simulating turbulent MHD flows with spatially varying magnetic fields in OpenFOAM, effectively eliminating the discrete cancellation pathology inherent in Maxwell-stress formulations. The work emphasizes that it evaluates the electromagnetic coupling and numerical formulation rather than proposing new turbulence physics; the standard closures are used unmodified to isolate numerical effects. The study establishes that the solver can handle engineering-scale, non-uniform magnetic fields without the spurious body forces that plague existing induction-based implementations, provided the low- $R_m$ approximation holds. The author notes limitations, including the lack of explicit magnetic terms in the turbulence closures, the restriction to insulating walls, and the first-order splitting error in the time-stepping scheme, defining these as directions for future work.

More like this