TORAX equation summary
TORAX simulates the time evolution of core plasma profiles using a coupled set of 1D transport PDEs, discretized in space and time, and solved numerically. The PDEs arise from flux-surfaced averaged moments of the underlying kinetic equations, and from Faraday’s law.
Governing equations
TORAX solves coupled 1D PDEs in normalized toroidal flux coordinates, \(\hat{\rho}\). The set of 1D PDEs being solved are as follows:
Ion heat transport, governing the evolution of the ion temperature \(T_i\).
\[\begin{split}\begin{split} \frac{3}{2} V'^{-5/3} \left(\frac{\partial }{\partial t}- \frac{\dot{\Phi}_b}{2\Phi_b}\frac{\partial}{\partial\hat{\rho}} \hat{\rho}\right)\left[V'^{5/3} n_i T_i\right] &= \\ \frac{1}{V'} \frac{\partial}{\partial \hat{\rho}} \left[ \chi_i n_i \frac{g_1}{V'} \frac{\partial T_i}{\partial \hat{\rho}} - g_0q_i^{\mathrm{conv}}T_i\right] &+ Q_i \end{split}\end{split}\]If multiple main ion species are present (e.g., a D-T mix), then \(n_i\) represents the sum of all main ions, and ion attributes like charge and mass are averaged values for the mixture, weighted by fractional abundance.
Electron heat transport, governing the evolution of the electron temperature \(T_e\).
\[\begin{split}\begin{split} \frac{3}{2} V'^{-5/3} \left(\frac{\partial }{\partial t}- \frac{\dot{\Phi}_b}{2\Phi_b}\frac{\partial}{\partial\hat{\rho}} \hat{\rho}\right)\left[V'^{5/3} n_e T_e\right] &= \\ \frac{1}{V'} \frac{\partial}{\partial \hat{\rho}} \left[ \chi_e n_e \frac{g_1}{V'} \frac{\partial T_e}{\partial \hat{\rho}} - g_0q_e^{\mathrm{conv}}T_e \right] &+ Q_e \end{split}\end{split}\]Electron particle transport, governing the evolution of the electron density \(n_e\).
\[\begin{split}\begin{split} \left(\frac{\partial}{\partial t}- \frac{\dot{\Phi}_b}{2\Phi_b}\frac{\partial} {\partial\hat{\rho}}\hat{\rho}\right)\left[ n_e V' \right] &= \\ \frac{\partial}{\partial \hat{\rho}} \left[D_e \frac{g_1}{V'} \frac{\partial n_e}{\partial \hat{\rho}} - g_0V_e n_e \right] &+ V'S_n \end{split}\end{split}\]Current diffusion, governing the evolution of the poloidal flux \(\psi\).
\[\begin{split}\begin{split} \frac{16 \pi^2 \sigma_{||}\mu_0 \hat{\rho} \Phi_b^2}{F^2} \left(\frac{\partial \psi}{\partial t}-\frac{\hat{\rho}\dot{\Phi}_b} {2\Phi_b}\frac{\partial \psi}{\partial \hat{\rho}}\right) &= \\ \frac{\partial}{\partial \hat{\rho}} \left( \frac{g_2 g_3}{\hat{\rho}} \frac{\partial \psi}{\partial \hat{\rho}} \right) &- \frac{8\pi^2 V' \mu_0 \Phi_b}{F^2} \langle \mathbf{B} \cdot \mathbf{j}_{n_i} \rangle \end{split}\end{split}\]
where \(T_{i,e}\) are ion and electron temperatures, \(n_{i,e}\) are ion and electron densities, and \(\psi\) is poloidal flux. \(\chi_{i,e}\) are ion and electron heat conductivities, \(q_{i,e}^{\mathrm{conv}}\) are ion and electron heat convection terms, \(D_e\) is electron particle diffusivity, and \(V_e\) is electron particle convection. \(Q_{i,e}\) are the total ion and electron heat source densities, and \(S_n\) is the total electron particle source density. \(V' \equiv dV/d\hat{\rho}\), i.e. the flux surface volume derivative. \(\sigma_{||}\) is the plasma neoclassical conductivity, and \(\langle \mathbf{B} \cdot \mathbf{j}_{n_i} \rangle\) is the flux-surface-averaged non-inductive current density (comprised of the bootstrap current and external current drive) projected onto and multiplied by the magnetic field. \(F \equiv RB_\varphi\) is the toroidal field flux function, where \(R\) is major radius, and \(B_\varphi\) is the toroidal magnetic field. \(\Phi_b\) is the toroidal flux enclosed by the last closed flux surface, and \(\dot{\Phi}_b\) is its time derivative, non-zero with time-varying toroidal magnetic field and/or last closed flux surface shape. \(\mu_0\) is the permeability of free space. The geometric quantities \(g_0\), \(g_1\), \(g_2\) and \(g_3\) are defined as follows.
where \(\nabla V\) is the radial derivative of the plasma volume, and \(\langle \rangle\) denotes flux surface averaging.
where \(R\) is the major radius along the flux surface being averaged.
The geometry terms \(V'\), \(g_0\), \(g_1\), \(g_2\) and \(g_3\) are calculated from pre-calculated flux-surface-averaged outputs of a Grad-Shafranov equilibrium code (see Physics models). Future work will allow for the option of self-consistent magnetic equilibrium calculated during the solver step.
TORAX can incorporate a dynamic magnetic equilibrium through inputting a sequence of precalculated magnetic geometries. Time derivative terms are calculated by interpolating the geometry data at the time points across a time interval at each solver step, and carrying out numeric differentiation.
Boundary conditions
All equations have a zero-derivative boundary condition at \(\hat{\rho}=0\). The \(T_i\), \(T_e\), \(n_e\) equations have Dirichlet (fixed) boundary conditions at \(\hat{\rho}=1\), which are user-defined. The \(\psi\) equation has a Neumann (derivative) boundary condition at \(\hat{\rho}=1\), which sets the total plasma current through the relation:
An alternative option for the \(\psi\) equation is to set a boundary condition based on a user-provided loop voltage at the last closed flux surface. This option sets the \(\psi\) Dirichlet boundary at time \(t+dt\). Total plasma current is then emergent from the evolution of the \(\psi\) equation.
Future work can extend the governing equations to include momentum transport, and multiple ion species including impurities. Details of the physics models underlying the PDE coefficients is provided in Physics models.