# Magnetohydrodynamic simulations in COMSOL Multiphysics

If you don’t know what magnetohydrodynamics (henceforce “MHD”, because I’m boring) is, this is basically all you need to know.

If you don’t know what COMSOL is… this.

In all seriousness, though, MHD is basically just Euler’s equations for fluid flow generalised to a highly conductive fluid where Maxwell’s equations also apply. It has a handful of applications, mostly relating to astrophysical or nuclear reactor plasmas. A simple example of a use case of MHD is in predicting stability and optimising design of magnetically confining fusion reactors like tokamaks (a type of reactor that has a fat donut-shaped form and spirally magnetic fields) without having to build them.

COMSOL Multiphysics is a finite-element analysis modelling program, and one of the most powerful in the business. It can simulate anything from chemical reactions in a supersonic flow to material stress in a changing magnetic field. Assuming, that is, you have the right modules. And you can find the right buttons in the endless submenus. The documentation for COMSOL is not the best. There’s a huge gallery of remarkably well-documented example simulations, but making anything outside of those requires a lot of troubleshooting with very little companionship.

So naturally, I decided to calculate the MHD equilibrium configuration of the tokamak known as JET using COMSOL’s AC/DC and fluid flow modules.

Depending on your familiarity with MHD and COMSOL (I can’t imagine you’re reading this if you know nothing of either), making a tokamak in COMSOL may seem like a strange thing to want to do. Well, it is. COMSOL is typically used for industrial-type things—systems with simple physics but complicated geometries where analytic solutions are infeasible and writing new code isn’t worth it.

Nuclear fusion reactors are neither of those things. MHD hardly qualifies as “simple physics”, but even it is only an approximation of plasma behaviour, based on assumptions that aren’t actually true in modern reactors. The actual dynamics are far more complex, making them ill-suited for COMSOL. The geometry of fusion reactors, on the other hand, is pretty simple—just an axisymmetric vacuum chamber with some carefully positioned coils. Nothing like the heat sinks and brackets for which COMSOL is usually used. Furthermore, analytic solutions actually kind of exist. Nothing exact or very general, and nothing that wasn’t the result of several PhD theses’ worth of work, but enough that partial solutions can usually be obtained far faster and more easily than via a full simulation. Despite this, full simulation codes also exist, specifically tailored to fusion reactors. TRANSP and NUBEAM are currently some of the more notable ones. These are extremely slow programs, but far more accurate than a simple MHD solver in COMSOL could ever be.

I did this anyway for several reasons. For one thing, I, unlike anyone else I know, like COMSOL. I thought it would be kind of cool to try to do MHD in it. For another, and more importantly, this was the final project of an independent study I did in MHD, so naturally the more complex Particle-in-Cell and Monte Carlo simulations that I really should be using were overscoped for my purposes. I didn’t have time to write code of my own, so making some MHD in COMSOL seemed like a fun, well-scoped, educational challenge.

## Geometry

The reactor after which I chose to model my simulation is the Joint European Torus, in England. The boundary of JET’s plasma can be represented by a rotational sweep of the curve

$R = R_0(1 + \epsilon \cos( t + \delta_0 \sin{t})$
$z = \kappa \epsilon \sin(t)$

where

$R_0 = 3 \mathrm{m};\ \epsilon = 0.33;\ \kappa = 1.7;\ \delta_0 = 0.25;\ 0 \le t \le 2 \pi$

Thus, the geometry was pretty straightforward to define. I went with a 2D axisymmetric model characterised by my parametric curve. A 2D model can’t fully capture all of the interesting plasma dynamics that go on in a tokamak, but then, neither can ideal MHD. As a first-pass model, two dimensions are enough. I surrounded this curve with a large cylinder of air, to give the electromagnetic fields space to resolve. I keep hearing that there’s a way to do an infinite boundary cell in COMSOL, but I don’t know how, so a large cylinder it is. I also added an arbitrarily thin wire in the centre to carry current for the toroidal magnetic field. In real tokamaks, there are coils that go all the way around the plasma for this, but you know, I don’t need your judgement.

## Physics

The basic equations of ideal MHD are as follows:

$\nabla \cdot \vec{B} = 0$
$\nabla \times \vec{B} = \mu_0 \vec{J}$
$\nabla \times \vec{E} = -\frac{\partial \vec{B}}{\partial t}$
$\vec{E} + \vec{v} \times \vec{B} = \eta \vec{J} \approx 0$
$\frac{\partial \rho}{\partial t} + \nabla \cdot (\rho \vec{v}) = 0$
$\rho \left(\frac{\partial \vec{v}}{\partial t} + \vec{v} \cdot \nabla \vec{v} \right) = - \nabla p + \vec{B} \times \vec{J}$
$\frac{3}{2} \left( \frac{\partial p}{\partial t} + \nabla \cdot (p \vec{v}) \right) = -p \nabla \cdot \vec{v} + \eta J^2 - \nabla \cdot \vec{q} \approx 0$

The first three are Maxwell’s equations. Gauss’s Law is omitted because ideal MHD doesn’t concern itself with charge density, and the plasma is assumed to be quasineutral. The fourth equation is Ohm’s Law for a perfect conductor. The last three are Navier-Stokes equations. Together, they fully define the dynamics of a perfectly conductive fluid. Specifically in my case, I cared about equilibrium configurations where $\frac{\partial}{\partial t}$ and $\vec{v} \rightarrow 0$.

At first, I thought since I’m simulating a plasma, that the plasma module would be the most appropriate way to go for this simulation. I even got COMSOL to give my school a free trial of the plasma module so that I could try it.

As it turns out, though, the plasma module is way overkill for simple MHD. It simulates different species of particles and the collisions and reactions between them, making it take upwards of 20 minutes to compute the simplest of problems. If you don’t get all the species and reactions just right, it won’t converge at all. After trying to make a box of plasma do absolutely noting for a few days, I decided to move to AC/DC and fluid flow.

### AC/DC

I started with a Coil domain in the plasma to represent the toroidal current, which produced my poloidal field. IRL, solenoids and other tricks are used to produce this current; putting one of those in COMSOL would probably have been more realistic than defining a constant current as I did, but would have also taken much more time and effort. I initially faced issues where the Coil would simply induce a contrary current in the space around the plasma, cancelling out any fields. I had to turn the conductivity of the outer region way down to prevent this.

I couldn’t figure out how to get the Coil domain to point axially instead of azimuthally, so for the toroidal field, I filled the central wire with an External Current Density domain. It was a little annoying to have to calculate the current density based on the area and total current, which I had to compute from the desired field ($\vec{J} = \frac{\frac{1}{\mu_0} 2 \pi B_0 R_0}{\pi r_\mathrm{wire}^2} \hat{z}$), but it worked easily enough.

The vertical field that many tokamaks employ to keep their plasma centred turned out to be surprisingly difficult. Magnetisation (a submenu of the Ampere’s Law domain) didn’t work out for various reasons, so I ended up defining a magnetic potential on all the outer boundaries as $\vec{A} = B_V\! R \hat{\theta}$. That did the trick, fully defining all of my magnetic fields. I got some pretty cool streamline plots out of that.

### Fluid Flow

The fluid flow turned out to be both surprisingly difficult and deceptively simple. Since I didn’t anticipate any plasma outside of that parametric domain, I limited the fluid flow physics to that region and set all the boundaries to Walls. The only other boundary condition I needed was a Pressure Point Constraint to fully define it.

Plasmas are characterised by extremely small viscosities; so small that treating them as fluids often isn’t even accurate, since velocities at a point are so poorly correlated. Nevertheless, I turned the viscosity way up, while also setting the fluid model to a Creeping Flow, to help it find equilibrium more quickly without wasting time on inertia and turbulence. After all, once it reaches equilibrium and there’s no velocity, the viscosity won’t matter, right?

In any case, my next step was to add the Lorentz force to my Navier-Stokes equations to couple the two modules. There’s a multiphysics called “Lorentz Force” for this very situation, but for whatever reason, I couldn’t get it to apply to Coil domains, so I had to do it manually. I simply added a Volume Force and plugged in the pre-calculated “mf.FLtzr” and “mf.FLtzz” values.

## Results and shortcomings

The resulting model was pretty bad.

If you think for a bit, you’ll quickly notice some important physics that are completely absent. I didn’t forget them; I just couldn’t figure out how to do them. For one thing, there’s no poloidal current, so the plasma has no toroidal magnetic response. As far as I can tell, these currents usually arise in real reactors due to the way it starts up, with the toroidal field coils (which I didn’t properly model in the first place) starting off and ramping up. The plasma induces poloidal currents to impede this, which continue indefinitely because the plasma is so conductive. While my model was time-dependent, I did not give it any time-dependent inputs, so such phenomena could not occur.

The second thing is related to how, as you may notice, my magnetic field lines, and thus plasma boundary, don’t line up with the parametric curve I drew. I’m not entirely sure how to rectify this, but I think I would need to start by implementing Faraday’s law with the fluid velocity. Currently, the magnetic field affects the fluid, but the fluid does not affect the magnetic field. Instead, currents should be induced whenever there’s a $\vec{v} \times \vec{B}$ to make the magnetic field lines move with the plasma. That would help in shaping the magnetic field and stop the plasma from moving so far.

The result is a plasma that always accelerates in the same direction, magnifying the same vortices indefinitely instead of finding equilibrium. I see the path forward, but don’t know enough of the technical details for fixing these problems to want to proceed. Given how long it took me to do such basic things as adding a uniform magnetic field within COMSOL, I don’t think completing the model here is feasible. If I really wanted to do a simulation like this, it would be more straightforward to write my own FEA or PIC.

## Conclusion

COMSOL is not the place to model a tokamak. Tokamaks are so well studied and so niche that working within the framework of existing analytical solutions and either using existing codes or making one from scratch makes way more sense than setting up all the physics in COMSOL just right.

If you have some more complex geometry and not too much going on in terms of weird magnetic fields, MHD in COMSOL might be worth a try. It won’t be easy—MHD is evidently not something COMSOL expects its users to need—but it will certainly be easier for you than it was for me if you don’t spend so much time on weird currents and magnetic fields.

As a project, I call this a success. I learned some fun COMSOL tricks, fastened my understanding of magnetohydrodynamics (What? I finished boring. You would not believe the sound the lathe was making the whole time!), and got better acquainted with the state of the art and with unique challenges in modelling for nuclear fusion. As far as actual applications, though, I see none. Hopefully if anyone ever finds themself in a situation as I did this semestre, this post can offer them some guidance.

This site uses Akismet to reduce spam. Learn how your comment data is processed.