A concept study in plasma compression
| Disc diameter | 0.40m |
| Separation | 0.25m |
| ω₁ (upper) | 1050rad/s |
| ω₂ (lower) | 950rad/s |
| B at disc face | 3.2T |
| B in compression zone | 9.7T |
| Lining | tungsten |
Two parallel tungsten-lined discs (diameter 0.40 m, separation 0.25 m) rotate in opposite directions at differential velocities: ω₁ = 1050 rad/s (upper, clockwise) and ω₂ = 950 rad/s (lower, counter-clockwise). The tungsten lining gives the field-shaping surfaces a sputtering threshold well above that of carbon divertor tiles.
Each rotating tungsten disc generates a toroidal magnetic field. The counter-rotation sets the two fields against each other so they form a compressed magnetic bottle in the central region. The design puts 3.2 T at the disc surfaces and 9.7 T in the compression zone.
The magnetohydrodynamic (MHD) interaction between opposing fields forms a high-compression torus at the geometric center. This zone is designed to reach a plasma density of 10²¹ m⁻³ at temperatures above 10 keV, the range the Lawson criterion sets for deuterium-tritium ignition.
The differential rotational velocity (Δω = 100 rad/s) induces a sheared poloidal plasma flow that works against the major MHD instabilities: tearing modes, kink modes, and edge-localized modes (ELMs). On this basis the model estimates a confinement time of 2.3 s at β = 8.5%.
We present a magnetic confinement concept built on two counter-rotating tungsten-lined plate discs whose opposed toroidal fields meet in a central compression zone. The discs are 0.40 m across and 0.25 m apart and turn at 1050 and 950 rad/s in opposite senses. The design targets a plasma density of 10²¹ m⁻³ at temperatures above 10 keV in that zone.
The tungsten lining gives the field-shaping surfaces a sputtering threshold well above that of carbon. From the Grad–Shafranov equilibrium, the β definition and the Lawson product set out below, the model estimates an energy confinement time of 2.3 s at β = 8.5%. The solver on this page relaxes the equilibrium live and reports q₀, β and the flux-surface shape from the converged ψ; the panel beside the Lawson block carries the triple product the targets give.
The differential rotation, Δω = 100 rad/s, drives a sheared azimuthal flow across the midplane, and the tearing growth-rate expression below carries the shear term that works against tearing and kink modes. The reconnection model beside it shows the plasmoid chain that shear acts on.
The force balance equation governing plasma equilibrium in the dual-disc configuration:
\nabla p = \mathbf{J} \times \mathbf{B}Where the current density \mathbf{J} is determined by Ampère's law:
\mathbf{J} = \frac{1}{\mu_0} \nabla \times \mathbf{B}For axisymmetric toroidal equilibrium, the magnetic flux function \psi satisfies:
R^2 \nabla \cdot \left(\frac{1}{R^2} \nabla \psi \right) = -\mu_0 R^2 \frac{dp}{d\psi} - \frac{1}{2} \frac{dF^2}{d\psi}Where R is the major radius, p(\psi) is the plasma pressure, and F(\psi) = RB_\phi is the poloidal current function.
With the Solov'ev choice p = p_0 \psi and F^2 = F_0^2 + 2 F_1 \psi the right-hand side is a fixed function of R, and the equation is linear in \psi. The solver beside this block relaxes it on an (R, Z) grid with successive over-relaxation and reports q_0, \beta and the shape of the last closed flux surface from the converged \psi. The first preset is the compression torus of this concept; the second is a compact tokamak for scale.
The plasma beta parameter, ratio of plasma pressure to magnetic pressure:
\beta = \frac{2\mu_0 \langle p \rangle}{B^2}Lawson criterion for breakeven fusion (D-T reaction):
n \tau_E T \geq 3 \times 10^{21} \text{ m}^{-3} \cdot \text{s} \cdot \text{keV}Where n is density, \tau_E is energy confinement time, and T is temperature. The panel beside this block evaluates both expressions: β from the solver's converged pressure and field, and the triple product from the design targets with p = 2 n k_B T.
The counter-rotating discs create a sheared azimuthal flow with velocity profile:
v_\phi(z) = \begin{cases} \omega_1 r & z > 0 \\ -\omega_2 r & z < 0 \end{cases}The shear enters the tearing-mode growth rate as a subtractive term:
\gamma = \frac{k_\parallel v_A}{\sqrt{2\pi}} \exp\left(-\frac{\Delta'}{2k_\theta r_s}\right) - \frac{(\omega_1 + \omega_2)k_\theta r}{2}The reconnection model below shows what the shear acts on: a current sheet tearing into a chain of plasmoids, with the single-mode growth rate scaling as S^{-3/5} and the plasmoid instability asS^{1/4}\, v_A / L (Loureiro, Schekochihin & Cowley 2007).
Conservation of magnetic flux during compression (frozen-in flux condition):
\frac{\partial \mathbf{B}}{\partial t} = \nabla \times (\mathbf{v} \times \mathbf{B})For adiabatic compression with compression ratio C = V_i/V_f:
B_f = C^{2/3} B_i, \quad T_f = C^{2/3} T_i, \quad n_f = C \cdot n_iThe tokamak confines a D-shaped plasma with a strong toroidal field from a ring of coils and a poloidal field from the plasma's own current. Compact high-field machines of the SPARC class reach 12 T on axis at a major radius under two metres. Professor Nuno Loureiro's work on magnetic reconnection and plasma turbulence informed the physics basis of this class of machine.
| Major radius R₀ | 1.85m |
| Minor radius a | 0.97m |
| Aspect ratio | 1.9 |
| B on axis | 12.2T |
| Plasma current | 8.7MA |
| Elongation κ | 1.8 |
| Triangularity δ | 0.5 |
| TF coils | 18 |
1import numpy as np2from scipy.sparse import diags3from scipy.sparse.linalg import spsolve45defsolve_grad_shafranov(R, Z, psi_boundary, pressure_profile, current_profile, max_iter=1000):6"""7 Solve the Grad-Shafranov equation for plasma equilibrium89 Parameters:10 -----------11 R, Z : ndarray12 Radial and vertical grid coordinates13 psi_boundary : ndarray14 Boundary conditions for magnetic flux15 pressure_profile : callable16 Pressure as a function of psi: p(psi)17 current_profile : callable18 Poloidal current function: F(psi)19 """20 nR, nZ = R.shape21 dR = R[1,0]- R[0,0]22 dZ = Z[0,1]- Z[0,0]2324 psi = np.zeros_like(R)25 psi[:,0]= psi_boundary[:,0]26 psi[:,-1]= psi_boundary[:,-1]27 psi[0,:]= psi_boundary[0,:]28 psi[-1,:]= psi_boundary[-1,:]2930for iteration inrange(max_iter):31 psi_old = psi.copy()3233for i inrange(1, nR-1):34for j inrange(1, nZ-1):35 R_ij = R[i, j]3637# Second derivatives38 d2psi_dR2 =(psi[i+1,j]-2*psi[i,j]+ psi[i-1,j])/ dR**239 d2psi_dZ2 =(psi[i,j+1]-2*psi[i,j]+ psi[i,j-1])/ dZ**240 dpsi_dR =(psi[i+1,j]- psi[i-1,j])/(2*dR)4142# Source term43 p = pressure_profile(psi[i,j])44 F = current_profile(psi[i,j])45 dp_dpsi =(pressure_profile(psi[i,j]+1e-6)- p)/1e-646 dF2_dpsi =(current_profile(psi[i,j]+1e-6)**2- F**2)/1e-64748 mu_0 =4*np.pi*1e-749 source =-mu_0 * R_ij**2* dp_dpsi -0.5* dF2_dpsi5051# Update psi using successive over-relaxation52 omega =1.5# Relaxation parameter53 psi_new =((d2psi_dR2 + d2psi_dZ2 - dpsi_dR/R_ij + source)*54(dR*dZ)**2/(2*(dR**2+ dZ**2)))55 psi[i,j]=(1-omega)*psi[i,j]+ omega*psi_new5657# Check convergence58 residual = np.max(np.abs(psi - psi_old))59if residual <1e-6:60print(f"Converged in {iteration} iterations")61break6263return psi6465# Calculate beta parameter66defcalculate_beta(density, temperature, magnetic_field):67"""Calculate plasma beta"""68 k_B =1.380649e-23# Boltzmann constant69 mu_0 =4*np.pi*1e-77071 plasma_pressure = density * k_B * temperature72 magnetic_pressure = magnetic_field**2/(2*mu_0)7374return plasma_pressure / magnetic_pressure
1#include<vector>2#include<cmath>3#include<random>45structParticle{6double x, y, z;// Position7double vx, vy, vz;// Velocity8double charge;// Charge9double mass;// Mass10};1112classPlasmaSimulator{13private:14 std::vector<Particle> particles;15double omega1, omega2;// Disc rotation rates16double B0;// Base magnetic field17double dt;// Timestep1819// Magnetic field at position20voidmagneticField(double x,double y,double z,21double& Bx,double& By,double& Bz){22double r = std::sqrt(x*x + y*y);2324if(z >0){25// Upper disc contribution26 Bx =-B0 * omega1 * y /1000;27 By = B0 * omega1 * x /1000;28 Bz = B0 *(1.0- std::exp(-std::abs(z)/0.1));29}else{30// Lower disc contribution31 Bx = B0 * omega2 * y /1000;32 By =-B0 * omega2 * x /1000;33 Bz =-B0 *(1.0- std::exp(-std::abs(z)/0.1));34}35}3637public:38PlasmaSimulator(int nParticles,double omega1_,double omega2_,double B0_)39:omega1(omega1_),omega2(omega2_),B0(B0_),dt(1e-9){4041// Initialize particles42 std::random_device rd;43 std::mt19937 gen(rd());44 std::normal_distribution<>dis(0.0,0.05);4546for(int i =0; i < nParticles; i++){47 Particle p;48 p.x =dis(gen);49 p.y =dis(gen);50 p.z =dis(gen);51 p.vx =dis(gen)*1e4;52 p.vy =dis(gen)*1e4;53 p.vz =dis(gen)*1e4;54 p.charge =1.602e-19;55 p.mass =1.672e-27;56 particles.push_back(p);57}58}5960voidstep(){61for(auto& p : particles){62double Bx, By, Bz;63magneticField(p.x, p.y, p.z, Bx, By, Bz);6465// Lorentz force: F = q(v × B)66double Fx = p.charge *(p.vy * Bz - p.vz * By);67double Fy = p.charge *(p.vz * Bx - p.vx * Bz);68double Fz = p.charge *(p.vx * By - p.vy * Bx);6970// Update velocity (Boris pusher for stability)71 p.vx +=(Fx / p.mass)* dt;72 p.vy +=(Fy / p.mass)* dt;73 p.vz +=(Fz / p.mass)* dt;7475// Update position76 p.x += p.vx * dt;77 p.y += p.vy * dt;78 p.z += p.vz * dt;79}80}8182doublecalculateAverageDensity(){83// Calculate density in central region84double volume =0.001;// 1 cm³85int count =0;86for(constauto& p : particles){87double r = std::sqrt(p.x*p.x + p.y*p.y + p.z*p.z);88if(r <0.1) count++;89}90return count / volume;91}92};
The full paper: the MHD model, the simulation set-up and the engineering specification of the concept
Download PDF