ColdElectronTemperature

The DREAM.Settings.Equations.ColdElectronTemperature.ColdElectronTemperature holds settings for the (cold) electron temperature. The unknown quantity is called T_cold and can be evolved either according to a prescribed function, or according to an energy balance equation.

The electron temperature and electron thermal energy are related through

\[W_{\rm cold} = \frac{3}{2}n_{\rm cold}T_{\rm cold}.\]

Prescribed evolution

The temperature can be prescribed in both space and time. The full syntax for prescribing the electron temperature is

ds.eqsys.T_cold.setPrescribedData(temperature=T_cold, radius=T_cold_r, times=T_cold_t)

where T_cold is an array of shape (nt, nr) specifying the cold electron temperature in space and time, T_cold_r is a vector with nr elements representing the radial grid on which the temperature is prescribed, and T_cold_t is a vector with nt elements representing the time grid on which the temperature is prescribed.

For convenience, it is possible to prescribe a constant and radially uniform temperature by only specifying it as a scalar:

ds.eqsys.T_cold.setPrescribedData(1000)

This will cause the temperature to be \(T_{\rm cold} = 1\,\mathrm{keV}\) at every radius, in every time step, during the simulation.

Self-consistent evolution

The temperature can be evolved self-consistently, i.e. by calculating the energy loss of the plasma through various processes. In its most general form, the energy balance equation solved by DREAM is

(1)\[\begin{split}\frac{\partial W_{\rm cold}}{\partial t} &= \frac{j_\Omega}{B}\left\langle \boldsymbol{E}\cdot\boldsymbol{B} \right\rangle - \left\langle n_{\rm cold} \right\rangle\sum_i\sum_{j=0}^{Z_i-1} n_i^{(j)}L_i^{(j)} + \left\langle Q_c \right\rangle +\\ &+ \frac{1}{V'}\frac{\partial}{\partial r}\left[ V'\frac{3\left\langle n_{\rm cold} \right\rangle}{2}\left( -A_WT_{\rm cold} + D_W\frac{\partial T_{\rm cold}}{\partial r} \right) \right]\end{split}\]

Here, the terms included in the most complete form of the energy balance equation are

• \(\frac{j_\Omega}{B}\langle\boldsymbol{E}\cdot\boldsymbol{B}\rangle\): Ohmic heating (energy transferred to the plasma from the electric field).

• \(\langle n_{\rm cold} \rangle\sum_i\sum_{j=0}^{Z_i-1} n_i^{(j)} L_i^{(j)}\): Energy loss by inelastic atomic processes. The rate coefficient \(L_i^{(j)}\) is given by \(L_i^{(j)} = L_{\rm line} + L_{\rm free} + \Delta W_i^{(j)}(I_i^{(j)} - R_i^{(j)})\), where \(L_{\rm line}\) is the radiated power by line emission, \(L_{\rm free}\) is the radiated power by recombination emission and bremsstrahlung, \(\Delta W_i^{(j)}\) is the ionization threshold (taken from NIST), and \(I_i^{(j)}\) and \(R_i^{(j)}\) are the rates of ionization and recombination for the given ion charge state. Unless otherwise specified, coefficients are taken from OpenADAS.

• \(\langle Q_c \rangle\): Collisional heat transfer to the cold electrons, from hot and runaway electrons, as well as ions:

  \[\begin{split}\left\langle Q_c \right\rangle &= \int\mathrm{d}p \int_{-1}^1\mathrm{d}\xi_0 \frac{\mathcal{V}'}{V'}\Delta\dot{E}_{ee} f_{\rm hot/re} + \sum_i Q_{ei},\\ Q_{kl} &= \frac{\left\langle nZ^2 \right\rangle_k\left\langle nZ^2 \right\rangle_l e^4\ln\Lambda_{kl}}{\left(2\pi\right)^{3/2}\epsilon_0^2 m_k m_l} \frac{T_k - T_l}{\left(\frac{T_k}{m_k} + \frac{T_l}{m_l} \right)^{3/2}},\\ \Delta\dot{E}_ee &= 4\pi n_{\rm cold} r_0^2 \ln\Lambda_{ee}\frac{m_ec^4}{v},\end{split}\]

• Radial transport of heat. DREAM uses and advection-diffusion heat transport model with either arbitrary coefficients \(A_W\) and \(D_W\), or a Rechester-Rosenbluth diffusion model requiring the magnetic perurbation level \(\delta B/B\) to be specified.

To solve for the electron temperature self-consistently, the user should call the method setType() with the argument TYPE_SELFCONSISTENT. Additionally, the temperature profile at \(t=0\) must be specified:

import DREAM.Settings.Equations.ColdElectronTemperature as Tcold

ds = DREAMSettings()
...
ds.eqsys.T_cold.setType(Tcold.TYPE_SELFCONSISTENT)
ds.eqsys.T_cold.setInitialProfile(temperature=T0, radius=T0_r)

If a uniform initial temperature profile is desired, the parameter temperature in the call to setInitialProfile() can be a scalar, in which case the parameter radius can also be omitted.

Recombination radiation

By default, radiation due to recombination is neglected in the energy balance equation (1) (i.e. the coefficient \(L_{\rm free}\) contains only the contribution from bremsstrahlung, and the coefficient \(R_i^{(j)}\) is dropped altogether). To account for radiation due to recombination, the user must call the method setRecombinationRadiation with the option RECOMBINATION_RADIATION_INCLUDED or True.

import DREAM.Settings.Equations.ColdElectronTemperature as Tcold

ds = DREAMSettings()
...
ds.eqsys.T_cold.setRecombinationRadiation(Tcold.RECOMBINATION_RADIATION_INCLUDED)
# or...
ds.eqsys.T_cold.setRecombinationRadiation(True)

Neutral Beam Injection (NBI)

The cold electron temperature equation can include heating from Neutral Beam Injection (NBI). The NBI settings are handled through the NBISettings class. This class is used to define the beam geometry, physical tokamak properties, and optional radial or time-dependant profiles for current density and beam power. To include NBI heating, the user must call the method nbi.setEnabled() with the argument True.

import DREAM.Settings.Equations.ColdElectronTemperature as Tcold

ds = DREAMSettings()
...
ds.eqsys.T_cold.nbi.setEnabled(True)

The beam geometry can be set using the methods nbi.setOrigin(), nbi.setDirection(), and nbi.setBeamParameters(), which requires the following parameters:

...
# Set beam origin point in X, Y, Z (in meters)
ds.eqsys.T_cold.nbi.setOrigin(P0 = [X0, Y0, Z0])

# Set beam direction
ds.eqsys.T_cold.nbi.setDirection(n = [nX, nY, nZ])

# Set beam physical parameters, including:
# r_beam: beam radius (in meters)
# Ti_beam: beam energy (in Joules)
# m_i_beam: mass of beam ions (in kg)
# s_max: maximum path length of beam in plasma (in meters)
ds.eqsys.T_cold.nbi.setBeamParameters(r_beam = beamRadius, Ti_beam = beamEnergy, m_i_beam = beamMass, s_max = beamPathLength)

The beam and tokamak geometry can be visualized using the methods nbi.visualize_3d_tokamak() and nbi.visualize_flux_surfaces_top_view(),

...
# Visualize beam and tokamak in 3D
ds.eqsys.T_cold.nbi.visualize_3d_tokamak()

# Visualize flux surfaces in top view. Needs the radial grid as input
ds.eqsys.T_cold.nbi.visualize_flux_surfaces_top_view(radialGrid = radialGrid)

For visualization, it is necessary to set up the minor and major radius of the tokamak, which is done using the method nbi.setRadius(),

...
# Set major and minor radius of tokamak (in meters)
ds.eqsys.T_cold.nbi.setR0_NBI( R0 = majorRadius, a = minorRadius )

The NBI power and current density profiles can be defined using the methods nbi.setPowerProfile() and nbi.setCurrentProfile(), respectively. The power profile specifies the total injected beam power as a function of time, allowing the beam power to vary during the simulation. If a scalar value is provided, the specified power will be applied throughout the entire simulation. The current profile definesthe radial distribution of beam current density.

Predefined current profile types are available for typical TCV or ITER configurations via the method nbi.setCurrentProfile(), which accepts profile_type values: NBI_PROFILE_TCV (1), NBI_PROFILE_ITER (2), or NBI_PROFILE_CUSTOM (3) for user-defined profiles. These automatically apply beam shapes consistent with the experimental setups of the respective devices.

The energy partition between beam components (full, half, and third energy) can be configured using the method nbi.setEnergyFractions(). NBI beams contain a mix of atomic (full energy) and molecular species that dissociate into half and third energy components.

from DREAM.Settings.Equations.NBISettings import NBI_PROFILE_TCV, NBI_PROFILE_ITER, NBI_PROFILE_CUSTOM

...
# Set time-dependent power profile (in Watts)
ds.eqsys.T_cold.nbi.setPowerProfile(P_NBI_t = timeArray, P_NBI_x = powerArray)

# Set current profile type (TCV, ITER, or custom)
ds.eqsys.T_cold.nbi.setCurrentProfile(profile_type = NBI_PROFILE_TCV)

# For custom current profile, provide j_B_t and j_B_x arrays
ds.eqsys.T_cold.nbi.setCurrentProfile(profile_type = NBI_PROFILE_CUSTOM,
                                       j_B_t = timeArray, j_B_x = currentArray)

# Set energy fractions (must sum to 1.0)
ds.eqsys.T_cold.nbi.setEnergyFractions(f_full = 0.5, f_half = 0.4, f_third = 0.1)

The discretization of the beam in the radial, polar, and beam direction can be set using the method nbi.setDiscretization(), which allows setting the number of points used in each direction. The default values are 25 radial points, 25 polar points, and 50 points along the beam direction. .. code-block:: python

… # Set beam discretization parameters ds.eqsys.T_cold.nbi.setDiscretization(n_beam_radius = 30, n_beam_theta = 30, n_beam_s = 60)

Class documentation