Components

This section describes the model components currently available.

Species density

The density of a species can be calculated in several different ways, and are usually needed by other components.

fixed_density

Set the density to a value which does not change in time. For example:

[d]
type = fixed_density, ...

AA = 2 # Atomic mass
charge = 0
density = 1e17 # In m^-3

Note that the density can be a function of x, y and z coordinates for spatial variation.

The implementation is in the FixedDensity class:

struct FixedDensity : public Component

Set ion density to a fixed value

Public Functions

inline FixedDensity(std::string name, Options &alloptions, Solver *solver)

Inputs

  • <name>

    • AA

    • charge

    • density value (expression) in units of m^-3

inline virtual void transform(Options &state) override

Sets in the state the density, mass and charge of the species

  • species

    • <name>

      • AA

      • charge

      • density

inline virtual void outputVars(Options &state) override

Add extra fields for output, or set attributes e.g docstrings.

evolve_density

This component evolves the species density in time, using the BOUT++ time integration solver. The species charge and atomic mass must be set, and the initial density should be specified in its own section:

[d]
type = evolve_density, ...

AA = 2 # Atomic mass
charge = 0

[Nd]
function = 1 - 0.5x # Initial condition, normalised to Nnorm

The equation solved is:

\[\frac{\partial n}{\partial t} = -\nabla\cdot\left[n \left(\frac{1}{B}\mathbf{b}\times\nabla\phi + v_{||}\mathbf{b}\right)\right] + S_n\]

where the source \(S_n\) is a combination of external source, and other processes that nay be included, including drift terms (e.g. magnetic drift) or atomic processes (e.g. ionization).

Notes:

  1. The density will be saved in the output file as N + species label, e.g Nd in the above example.

  2. If diagnose=true is set in the species options then the net source \(S_n\) is saved as SN + species, e.g. SNd; the external source is saved as S + species + _src e.g. Sd_src. The time derivative of density is saved as ddt(N + species + ) e.g. ddt(Nd).

  3. The density source can be set in the input mesh file as a field S + species + _src e.g. Sd_src. This can be overridden by specifying the source in the input options.

  4. The poloidal_flows switch controls whether the X-Y components of the ExB flow are included (default is true). If set to false then ExB flows are only in the X-Z plane.

The implementation is in the EvolveDensity class:

struct EvolveDensity : public Component

Evolve species density in time

Mesh inputs

N<name>_src A source of particles, per cubic meter per second. This can be over-ridden by the source option setting.

Public Functions

EvolveDensity(std::string name, Options &options, Solver *solver)

Inputs

  • <name>

    • charge Particle charge e.g. hydrogen = 1

    • AA Atomic mass number e.g. hydrogen = 1

    • bndry_flux Allow flow through radial boundaries? Default is true.

    • poloidal_flows Include poloidal ExB flows? Default is true.

    • density_floor Minimum density floor. Default is 1e-5 normalised units

    • low_n_diffuse Enhance parallel diffusion at low density? Default false

    • hyper_z Hyper-diffusion in Z. Default off.

    • evolve_log Evolve logarithm of density? Default false.

    • diagnose Output additional diagnostics?

  • N<name> e.g. “Ne”, “Nd+”

    • source Source of particles [/m^3/s] NOTE: This overrides mesh input N<name>_src

    • source_only_in_core Zero the source outside the closed field-line region?

    • neumann_boundary_average_z Apply Neumann boundaries with Z average?

virtual void transform(Options &state) override

This sets in the state

  • species

    • <name>

      • AA

      • charge

      • density

virtual void finally(const Options &state) override

Calculate ddt(N).

Requires state components

  • species

    • <name>

      • density

Optional components

  • species

    • <name>

      • velocity If included, requires sound_speed or temperature

      • density_source

  • fields

    • phi If included, ExB drift is calculated

virtual void outputVars(Options &state) override

Add extra fields for output, or set attributes e.g docstrings.

upstream_density_feedback

This is intended for 1D simulations, where the density at \(y=0\) is set by adjusting an input source. This component uses a PI controller method to scale the density source up and down, to maintain the specified upstream density. The source, e.g. Sd+_feedback, is calculated as a product of the control signal density_source_multiplier, and the array density_source_shape which defines the source region. The signal is non-dimensional and the controller depends on the value of density_source_shape to have a good initial guess of the source. It should be set to a reasonable value in the units of [m-3s-1]. A good reasonable value is the expected steady state domain particle loss (for example due to unrecycled ions at the target).

For example:

[d+]
type = ..., upstream_density_feedback

density_upstream = 1e19      # Density in m^-3
density_controller_p = 1e-2  # Feedback controller proportional (p) parameter
density_controller_i = 1e-3  # Feedback controller integral (i) parameter

[Nd+]
source_shape = h(pi - y) * 1e20  # Source shape

There are two additional settings which can make the controller more robust without excessive tuning:

density_source_positive ensures the controller never takes particles away, which can prevent oscillatory behaviour. Note that this requires some other domain particle sink to ensure control, or else the particle count can never reduce.

density_integral_positive This makes sure the integral component only adds particles. The integral component takes a long time to change value, which can result in large overshoots if the initial guess was too small. This setting mitigates this by disabling the integral term if the density is above the desired value.

Notes:

  • The example cases have their PI parameters tuned properly without the need of the above two settings.

  • Under certain conditions, the use of the PI controller can make the upstream density enter a very small oscillation (~0.05% of upstream value).

  • There is a separate source setting that includes a fixed (non varying) density source.

The implementation is in the UpstreamDensityFeedback class:

struct UpstreamDensityFeedback : public Component

Adds a time-varying density source, depending on the difference between the upstream density at y=0 and the specified value

Public Functions

inline UpstreamDensityFeedback(std::string name, Options &alloptions, Solver*)

Inputs

  • <name> (e.g. “d+”)

    • density_upstream Upstream density (y=0) in m^-3

    • density_controller_p Feedback proportional to error

    • density_controller_i Feedback proportional to error integral

    • density_integral_positive Force integral term to be positive? (default: false)

    • density_source_positive Force density source to be positive? (default: true)

    • diagnose Output diagnostic information?

  • N<name> (e.g. “Nd+”)

    • source_shape The initial source that is scaled by a time-varying factor

virtual void transform(Options &state) override

Inputs

  • <name>

    • density

Outputs

  • <name>

    • density_source

inline virtual void outputVars(Options &state) override

Add extra fields for output, or set attributes e.g docstrings.

inline virtual void restartVars(Options &state) override

Add extra fields to restart files.

fixed_fraction_ions

This sets the density of a species to a fraction of the electron density.

quasineutral

This component sets the density of one species, so that the overall charge density is zero everywhere. This must therefore be done after all other charged species densities have been calculated. It only makes sense to use this component for species with a non-zero charge.

Species pressure and temperature

isothermal

Sets the temperature of a species to a fixed value which is constant in space and time. If the species density is set then this component also calculates the pressure.

By default only saves the temperature once as a non-evolving variable. If diagnose is set then pressure is also saved as a time-evolving variable.

[e]
type = ..., isothermal

temperature = 10   # Constant temperature [eV]
struct Isothermal : public Component

Set temperature to a fixed value

Public Functions

virtual void transform(Options &state) override

Inputs

  • species

    • <name>

      • density (optional)

Sets in the state

  • species

    • <name>

      • temperature

      • pressure (if density is set)

virtual void outputVars(Options &state) override

Add extra fields for output, or set attributes e.g docstrings.

fixed_temperature

Sets the temperature of a species to a fixed value which is constant in time but can vary in space. If the species density is set then this component also calculates the pressure.

By default only saves the temperature once as a non-evolving variable. If diagnose is set then pressure is also saved as a time-evolving variable.

[e]
type = ..., fixed_temperature

temperature = 10 - x   # Spatially dependent temperature [eV]
struct FixedTemperature : public Component

Set species temperature to a fixed value

Public Functions

inline FixedTemperature(std::string name, Options &alloptions, Solver *solver)

Inputs

  • <name>

    • temperature value (expression) in units of eV

inline virtual void transform(Options &state) override

Sets in the state the temperature and pressure of the species

Inputs

  • species

    • <name>

      • density (optional)

Sets in the state

  • species

    • <name>

      • temperature

      • pressure (if density is set)

inline virtual void outputVars(Options &state) override

Add extra fields for output, or set attributes e.g docstrings.

evolve_pressure

Evolves the pressure in time. This pressure is named P where <species> is the short name of the evolving species e.g. Pe.

By default parallel thermal conduction is included, which requires a collision time. If collisions are not calculated, then thermal conduction should be turned off by setting thermal_conduction = false in the input options.

If the component option diagnose = true then additional fields will be saved to the dump files: The species temperature T + name (e.g. Td+ or Te), the time derivative ddt(P + name) (e.g. ddt(Pd+) or ddt(Pe)), and the source of pressure from other components is saved as SP + name (e.g. SPd+ or SPe). The pressure source is the energy density source multiplied by 2/3 (i.e. assumes a monatomic species).

\[\frac{\partial P}{\partial t} = -\nabla\cdot\left(P\mathbf{v}\right) - \frac{2}{3} P \nabla\cdot\mathbf{b}v_{||} + \frac{2}{3}\nabla\cdot\left(\kappa_{||}\mathbf{b}\mathbf{b}\cdot\nabla T\right) + \frac{2}{3}S_E + S_N\frac{1}{2}mNV^2\]

where \(S_E\) is the energy_source (thermal energy source), and \(S_N\) is the density source.

Notes:

  • Heat conduction through the boundary is turned off currently. This is because heat losses are usually calculated at the sheath, so any additional heat conduction would be in addition to the sheath heat transmission already included.

The implementation is in EvolvePressure:

struct EvolvePressure : public Component

Evolves species pressure in time

Mesh inputs

P<name>_src A source of pressure, in Pascals per second This can be over-ridden by the source option setting.

Public Functions

EvolvePressure(std::string name, Options &options, Solver *solver)

Inputs

  • <name>

    • bndry_flux Allow flows through radial boundaries? Default is true

    • density_floor Minimum density floor. Default 1e-5 normalised units.

    • diagnose Output additional diagnostic fields?

    • evolve_log Evolve logarithm of pressure? Default is false

    • hyper_z Hyper-diffusion in Z

    • kappa_coefficient Heat conduction constant. Default is 3.16 for electrons, 3.9 otherwise

    • kappa_limit_alpha Flux limiter, off by default.

    • poloidal_flows Include poloidal ExB flows? Default is true

    • precon Enable preconditioner? Note: solver may not use it even if enabled.

    • p_div_v Use p * Div(v) form? Default is v * Grad(p) form

    • thermal_conduction Include parallel heat conduction? Default is true

  • P<name> e.g. “Pe”, “Pd+”

    • source Source of pressure [Pa / s]. NOTE: This overrides mesh input P<name>_src

    • source_only_in_core Zero the source outside the closed field-line region?

    • neumann_boundary_average_z Apply Neumann boundaries with Z average?

virtual void transform(Options &state) override

Inputs

  • species

    • <name>

      • density

Sets

  • species

    • <name>

      • pressure

      • temperature Requires density

virtual void finally(const Options &state) override

Optional inputs

  • species

    • <name>

      • velocity. Must have sound_speed or temperature

      • energy_source

      • collision_rate (needed if thermal_conduction on)

  • fields

    • phi Electrostatic potential -> ExB drift

virtual void outputVars(Options &state) override

Add extra fields for output, or set attributes e.g docstrings.

virtual void precon(const Options &state, BoutReal gamma) override

Preconditioner

evolve_energy

Note This is currently under development and has some unresolved issues with boundary conditions. Only for testing purposes.

This evolves the sum of species internal energy and parallel kinetic energy, \(\mathcal{E}\):

\[\mathcal{E} = \frac{1}{\gamma - 1} P + \frac{1}{2}m nv_{||}^2\]

Note that this component requires the parallel velocity \(v_{||}\) to calculate the pressure. It must therefore be listed after a component that sets the velocity, such as evolve_momentum:

[d]
type = ..., evolve_momentum, evolve_energy

The energy density will be saved as E (e.g Ed) and the pressure as P (e.g. Pd). Additional diagnostics, such as the temperature, can be saved by setting the option diagnose = true.

struct EvolveEnergy : public Component

Evolves species internal energy in time

Mesh inputs

P<name>_src A source of pressure, in Pascals per second This can be over-ridden by the source option setting.

Public Functions

EvolveEnergy(std::string name, Options &options, Solver *solver)

Inputs

  • <name>

    • bndry_flux Allow flows through radial boundaries? Default is true

    • density_floor Minimum density floor. Default 1e-5 normalised units.

    • diagnose Output additional diagnostic fields?

    • evolve_log Evolve logarithm of pressure? Default is false

    • hyper_z Hyper-diffusion in Z

    • kappa_coefficient Heat conduction constant. Default is 3.16 for electrons, 3.9 otherwise

    • kappa_limit_alpha Flux limiter, off by default.

    • poloidal_flows Include poloidal ExB flows? Default is true

    • precon Enable preconditioner? Note: solver may not use it even if enabled.

    • thermal_conduction Include parallel heat conduction? Default is true

  • E<name> e.g. “Ee”, “Ed+”

    • source Source of energy [W / s]. NOTE: This overrides mesh input P<name>_src

    • source_only_in_core Zero the source outside the closed field-line region?

    • neumann_boundary_average_z Apply Neumann boundaries with Z average?

virtual void transform(Options &state) override

Inputs

  • species

    • <name>

      • density

      • velocity

Sets

  • species

    • <name>

      • pressure

      • temperature

virtual void finally(const Options &state) override

Optional inputs

  • species

    • <name>

      • velocity. Must have sound_speed or temperature

      • energy_source

      • collision_rate (needed if thermal_conduction on)

  • fields

    • phi Electrostatic potential -> ExB drift

virtual void outputVars(Options &state) override

Add extra fields for output, or set attributes e.g docstrings.

virtual void precon(const Options &state, BoutReal gamma) override

Preconditioner

SNB nonlocal heat flux

Calculates the divergence of the electron heat flux using the Shurtz-Nicolai-Busquet (SNB) model. Uses the BOUT++ implementation which is documented here.

struct SNBConduction : public Component

Calculate electron heat flux using the Shurtz-Nicolai-Busquet (SNB) model

This component will only calculate divergence of heat flux for the electron (e) species.

Usage

Add as a top-level component after both electron temperature and collision times have been calculated.

Important: If evolving electron pressure, disable thermal conduction or that will continue to add Spitzer heat conduction.

[hermes]
components = e, ..., collisions, snb_conduction

[e]
type = evolve_pressure, ...
thermal_conduction = false # For evolve_pressure

[snb_conduction]
diagnose = true # Saves heat flux diagnostics

Useful references:

Public Functions

inline SNBConduction(std::string name, Options &alloptions, Solver*)

Inputs

  • <name>

    • diagnose Saves Div_Q_SH and Div_Q_SNB

virtual void transform(Options &state) override

Inputs

  • species

    • e

      • density

      • collision_frequency

Sets

  • species

    • e

      • energy_source

virtual void outputVars(Options &state) override

Add extra fields for output, or set attributes e.g docstrings.

Species parallel dynamics

fixed_velocity

Sets the velocity of a species to a fixed value which is constant in time but can vary in space. If the species density is set then this component also calculates the momentum.

Saves the temperature once as a non-evolving variable.

[e]
type = ..., fixed_velocity

velocity = 10 + sin(z)   # Spatially dependent velocity [m/s]
struct FixedVelocity : public Component

Set parallel velocity to a fixed value

Public Functions

inline virtual void transform(Options &state) override

This sets in the state

  • species

    • <name>

      • velocity

      • momentum

inline virtual void outputVars(Options &state) override

Add extra fields for output, or set attributes e.g docstrings.

evolve_momentum

Evolves the momentum NV in time. The evolving quantity includes the atomic mass number, so should be divided by AA to obtain the particle flux.

If the component option diagnose = true then additional fields will be saved to the dump files: The velocity V + name (e.g. Vd+ or Ve), the time derivative ddt(NV + name) (e.g. ddt(NVd+) or ddt(NVe)), and the source of momentum density (i.e force density) from other components is saved as SNV + name (e.g. SNVd+ or SNVe).

The implementation is in EvolveMomentum:

struct EvolveMomentum : public Component

Evolve parallel momentum.

Public Functions

virtual void transform(Options &state) override

This sets in the state

  • species

    • <name>

      • momentum

      • velocity if density is defined

virtual void finally(const Options &state) override

Calculate ddt(NV).

Inputs

  • species

    • <name>

      • density

      • velocity

      • pressure (optional)

      • momentum_source (optional)

      • sound_speed (optional, used for numerical dissipation)

      • temperature (only needed if sound_speed not provided)

  • fields

    • phi (optional)

virtual void outputVars(Options &state) override

Add extra fields for output, or set attributes e.g docstrings.

zero_current

This calculates the parallel flow of one charged species so that there is no net current, using flows already calculated for other species. It is used like quasineutral:

[hermes]
components = h+, ..., e, ...   # Note: e after all other species

[e]
type = ..., zero_current,... # Set e:velocity

charge = -1 # Species must have a charge

electron_force_balance

This calculates a parallel electric field which balances the electron pressure gradient and other forces on the electrons (including collisional friction, thermal forces):

\[E_{||} = \left(-\nabla p_e + F\right) / n_e\]

where \(F\) is the momentum_source for the electrons. This electric field is then used to calculate a force on the other species:

\[F_z = Z n_z E_{||}\]

which is added to the ion’s momentum_source.

The implementation is in ElectronForceBalance:

struct ElectronForceBalance : public Component

Balance the parallel electron pressure gradient against the electric field. Use this electric field to calculate a force on the other species

E = (-∇p_e + F) / n_e

where F is the momentum source for the electrons.

Then uses this electric field to calculate a force on all ion species.

Note: This needs to be put after collisions and other components which impose forces on electrons

Public Functions

virtual void transform(Options &state) override

Required inputs

  • species

    • e

      • pressure

      • density

      • momentum_source [optional] Asserts that charge = -1

Sets in the input

  • species

    • <all except e> if both density and charge are set

      • momentum_source

electron_viscosity

Calculates the Braginskii electron parallel viscosity, adding a force (momentum source) to the electron momentum equation:

\[F = \sqrt{B}\nabla\cdot\left[\frac{\eta_e}{B}\mathbf{b}\mathbf{b}\cdot\nabla\left(\sqrt{B}V_{||e}\right)\right]\]

The electron parallel viscosity is

\[\eta_e = \frac{4}{3} 0.73 p_e \tau_e\]

where \(\tau_e\) is the electron collision time. The collisions between electrons and all other species therefore need to be calculated before this component is run:

[hermes]
components = ..., e, ..., collisions, electron_viscosity
struct ElectronViscosity : public Component

Electron viscosity

Adds Braginskii parallel electron viscosity, with SOLPS-style viscosity flux limiter

Needs to be calculated after collisions, because collision frequency is used to calculate parallel viscosity

References

Public Functions

ElectronViscosity(std::string name, Options &alloptions, Solver*)

Braginskii electron viscosity.

Inputs

  • <name>

    • diagnose: bool, default false Output diagnostic SNVe_viscosity?

    • eta_limit_alpha: float, default -1.0 Flux limiter coefficient. < 0 means no limiter

virtual void transform(Options &state) override

Inputs

  • species

    • e

      • pressure (skips if not present)

      • velocity (skips if not present)

      • collision_frequency

Sets in the state

  • species

    • e

      • momentum_source

virtual void outputVars(Options &state) override

Add extra fields for output, or set attributes e.g docstrings.

ion_viscosity

Adds ion viscosity terms to all charged species that are not electrons. The collision frequency is required so this is a top-level component that must be calculated after collisions:

[hermes]
components =  ..., collisions, ion_viscosity

By default only the parallel diffusion of momentum is included, adding a force to each ion’s momentum equation:

\[F = \sqrt{B}\nabla\cdot\left[\frac{\eta_i}{B}\mathbf{b}\mathbf{b}\cdot\nabla\left(\sqrt{B}V_{||i}\right)\right]\]

The ion parallel viscosity is

\[\eta_i = \frac{4}{3} 0.96 p_i \tau_i\]

If the perpendicular option is set:

[ion_viscosity]
perpendicular = true # Include perpendicular flows

Then the ion scalar viscous pressure is calculated as:

\[\Pi_{ci} = \Pi_{ci||} + \Pi_{ci\perp}\]

where \(\Pi_{ci||}\) corresponds to the parallel diffusion of momentum above.

\[\Pi_{ci||} = - 0.96 \frac{2p_i\tau_i}{\sqrt{B}} \partial_{||}\left(\sqrt{B} V_{||i}\right)\]

The perpendicular part is calculated from:

\[\begin{split}\begin{aligned}\Pi_{ci\perp} =& 0.96 p_i\tau_i \kappa \cdot \left[\mathbf{V}_E + \mathbf{V}_{di} + 1.16\frac{\mathbf{b}\times\nabla T_i}{B} \right] \\ =& -0.96 p_i\tau_i\frac{1}{B}\left(\mathbf{b}\times\kappa\right)\cdot\left[\nabla\phi + \frac{\nabla p_i}{en_i} + 1.61\nabla T_i \right]\end{aligned}\end{split}\]

A parallel force term is added, in addition to the parallel viscosity above:

\[F = -\frac{2}{3}B^{3/2}\partial_{||}\left(\frac{\Pi_{ci\perp}}{B^{3/2}}\right)\]

In the vorticity equation the viscosity appears as a divergence of a current:

\[\mathbf{J}_{ci} = \frac{\Pi_{ci}}{2}\nabla\times\frac{\mathbf{b}}{B} - \frac{1}{3}\frac{\mathbf{b}\times\nabla\Pi_{ci}}{B}\]

that transfers energy between ion internal energy and \(E\times B\) energy:

\[\begin{split}\begin{aligned}\frac{\partial \omega}{\partial t} =& \ldots + \nabla\cdot\mathbf{J}_{ci} \\ \frac{\partial p_i}{\partial t} =& \ldots - \mathbf{J}_{ci}\cdot\nabla\left(\phi + \frac{p_i}{n_0}\right)\end{aligned}\end{split}\]

Note that the sum of the perpendicular and parallel contributions to the ion viscosity act to damp the net poloidal flow. This can be seen by assuming that \(\phi\), \(p_i\) and \(T_i\) are flux functions. We can then write:

\[\Pi_{ci\perp} = -0.96 p_i\tau_i \frac{1}{B}\left(\mathbf{b}\times\kappa\right)\cdot\nabla\psi F\left(\psi\right)\]

where

\[F\left(\psi\right) = \frac{\partial\phi}{\partial\psi} + \frac{1}{en}\frac{\partial p_i}{\partial\psi} + 1.61\frac{\partial T_i}{\partial\psi}\]

Using the approximation

\[\left(\mathbf{b}\times\kappa\right)\cdot\nabla\psi \simeq -RB_\zeta \partial_{||}\ln B\]

expanding:

\[\frac{2}{\sqrt{B}}\partial_{||}\left(\sqrt{B}V_{||i}\right) = 2\partial_{||}V_{||i} + V_{||i}\partial_{||}\ln B\]

and neglecting parallel gradients of velocity gives:

\[\Pi_{ci} \simeq 0.96 p_i\tau_i \left[ \frac{RB_{\zeta}}{B}F\left(\psi\right) - V_{||i} \right]\partial_{||}\ln B\]

Notes and implementation details: - The magnitude of \(\Pi_{ci\perp}\) and \(\Pi_{ci||}\) are individually

limited to be less than or equal to the scalar pressure \(Pi\) (though can have opposite sign). The reasoning is that if these off-diagonal terms become large then the model is likely breaking down. Occasionally happens in low-density regions.

struct IonViscosity : public Component

Ion viscosity terms

Adds a viscosity to all species which are not electrons

Uses Braginskii collisional form, combined with a SOLPS-like flux limiter.

Needs to be calculated after collisions, because collision frequency is used to calculate parallel viscosity

The ion stress tensor Pi_ci is split into perpendicular and parallel pieces:

Pi_ci = Pi_ciperp + Pi_cipar

In the parallel ion momentum equation the Pi_cipar term is solved as a parallel diffusion, so is treated separately All other terms are added to Pi_ciperp, even if they are not really parallel parts

Public Functions

IonViscosity(std::string name, Options &alloptions, Solver*)

Inputs

  • <name>

    • eta_limit_alpha: float, default -1 Flux limiter coefficient. < 0 means off.

    • perpendicular: bool, default false Include perpendicular flows? Requires curvature vector and phi potential

virtual void transform(Options &state) override

Inputs

  • species

    • <name> (skips “e”)

      • pressure (skips if not present)

      • velocity (skips if not present)

      • collision_frequency

Sets in the state

  • species

    • <name>

      • momentum_source

virtual void outputVars(Options &state) override

Save variables to the output.

simple_conduction

This is a simplified parallel heat conduction model that can be used when a linearised model is needed. If used, the thermal conduction term in evolve_pressure component should be disabled.

[hermes]
components = e, ...

[e]
type = evolve_pressure, simple_conduction

thermal_conduction = false  # Disable term in evolve_pressure

To linearise the heat conduction the temperature and density used in calculating the Coulomb logarithm and heat conduction coefficient can be fixed by specifying conduction_temperature and conduction_density.

Note: For hydrogenic plasmas this produces very similar parallel electron heat conduction as the evolve_pressure term with electron-electron collisions disabled.

struct SimpleConduction : public Component

Simplified models of parallel heat conduction

Intended mainly for testing.

Expressions taken from: https://farside.ph.utexas.edu/teaching/plasma/lectures1/node35.html

Public Functions

inline virtual void transform(Options &state) override

Modify the given simulation state All components must implement this function

Drifts and transport

The ExB drift is included in the density, momentum and pressure evolution equations if potential is calculated. Other drifts can be added with the following components.

diamagnetic_drift

Adds diamagnetic drift terms to all species’ density, pressure and parallel momentum equations. Calculates the diamagnetic drift velocity as

\[\mathbf{v}_{dia} = \frac{T}{q} \nabla\times\left(\frac{\mathbf{b}}{B}\right)\]

where the curvature vector \(\nabla\times\left(\frac{\mathbf{b}}{B}\right)\) is read from the bxcv mesh input variable.

struct DiamagneticDrift : public Component

Calculate diamagnetic flows.

Public Functions

virtual void transform(Options &state) override

For every species, if it has:

  • temperature

  • charge

Modifies:

  • density_source

  • energy_source

  • momentum_source

polarisation_drift

This calculates the polarisation drift of all charged species, including ions and electrons. It works by approximating the drift as a potential flow:

\[\mathbf{v}_{pol} = - \frac{m}{q B^2} \nabla_\perp\phi_{pol}\]

where \(\phi_{pol}\) is approximately the time derivative of the electrostatic potential \(\phi\) in the frame of the fluid, with an ion diamagnetic contribution. This is calculated by inverting a Laplacian equation similar to that solved in the vorticity equation.

This component needs to be run after all other currents have been calculated. It marks currents as used, so out-of-order modifications should raise errors.

See the examples/blob2d-vpol example, which contains:

[hermes]
components = e, vorticity, sheath_closure, polarisation_drift

[polarisation_drift]
diagnose = true

Setting diagnose = true saves DivJ to the dump files with the divergence of all currents except polarisation, and phi_pol which is the polarisation flow potential.

struct PolarisationDrift : public Component

Calculates polarisation drift terms for all charged species, both ions and electrons.

Approximates the polarisation drift by a generalised flow potential phi_pol

v_pol = - (A / (Z * B^2)) * Grad_perp(phi_pol)

phi_pol is approximately the time derivative of the electric potential in the frame of the flow, plus an ion diamagnetic contribution

phi_pol is calculated using:

Div(mass_density / B^2 * Grad_perp(phi_pol)) = Div(Jpar) + Div(Jdia) + …

Where the divergence of currents on the right is calculated from:

  • species[…][“momentum”] The parallel momentum of charged species

  • DivJdia, diamagnetic current, calculated in vorticity component

  • DivJcol collisional current, calculated in vorticity component

  • DivJextra Other currents, eg. 2D parallel closures

The mass_density quantity is the sum of density * atomic mass for all charged species (ions and electrons)

Public Functions

virtual void transform(Options &state) override

Inputs

  • species

    • … All species with both charge and mass

      • AA

      • charge

      • density

      • momentum (optional)

  • fields

    • DivJextra (optional)

    • DivJdia (optional)

    • DivJcol (optional)

Sets

  • species

    • … All species with both charge and mass

      • density_source

      • energy_source (if pressure set)

      • momentum_source (if momentum set)

virtual void outputVars(Options &state) override

Add extra fields for output, or set attributes e.g docstrings.

anomalous_diffusion

Adds cross-field diffusion of particles, momentum and energy to a species.

[hermes]
components = e, ...

[e]
type = evolve_density, evolve_momentum, evolve_pressure, anomalous_diffusion

anomalous_D = 1.0   # Density diffusion [m^2/s]
anomalous_chi = 0,5 # Thermal diffusion [m^2/s]
anomalous_nu = 0.5  # Kinematic viscosity [m^2/s]

Anomalous diffusion coefficients can be functions of x and y. The coefficients can also be read from the mesh input file: If the mesh file contains D_ + the name of the species, for example D_e then it will be read and used as the density diffusion coefficient. Similarly, chi_e is the thermal conduction coefficient, and nu_e is the kinematic viscosity. All quantities should be in SI units of m^2/s. Values that are set in the options (as above) override those in the mesh file.

The sources of particles \(S\), momentum \(F\) and energy \(E\) are calculated from species density \(N\), parallel velocity \(V\) and temperature \(T\) using diffusion coefficients \(D\), \(\chi\) and \(\nu\) as follows:

\[\begin{split}\begin{aligned} S =& \nabla\cdot\left(D \nabla_\perp N\right) \\ F =& \nabla\cdot\left(m V D \nabla_\perp N\right) + \nabla\cdot\left(m N \nu \nabla_\perp V\right)\\ E =& \nabla\cdot\left(\frac{3}{2}T D \nabla_\perp N\right) + \nabla\cdot\left(N \chi \nabla_\perp T\right) \end{aligned}\end{split}\]

Note that particle diffusion is treated as a density gradient-driven flow with velocity \(v_D = -D \nabla_\perp N / N\).

struct AnomalousDiffusion : public Component

Add anomalous diffusion of density, momentum and energy

Mesh inputs

D_<name>, chi_<name>, nu_<name> e.g D_e, chi_e, nu_e

in units of m^2/s

Public Functions

AnomalousDiffusion(std::string name, Options &alloptions, Solver*)

Inputs

  • <name>

    • anomalous_D This overrides D_<name> mesh input

    • anomalous_chi This overrides chi_<name>

    • anomalous_nu Overrides nu_<name>

    • anomalous_sheath_flux Allow anomalous flux into sheath?

virtual void transform(Options &state) override

Inputs

  • species

    • <name>

      • density

      • temperature (optional)

      • velocity (optional)

Sets in the state

  • species

    • <name>

      • density_source

      • momentum_source

      • energy_source

virtual void outputVars(Options &state) override

Add extra fields for output, or set attributes e.g docstrings.

Neutral gas models

The neutral_mixed component solves fluid equations along \(y\) (parallel to the magnetic field), and uses diffusive transport in \(x\) and \(z\). It was adopted from the approach used in UEDGE and this paper [Journal of Nuclear Materials, vol. 313-316, pp. 559-563 (2003)].

\[\begin{split}\begin{aligned}\frac{\partial n_n}{\partial t} =& -\nabla\cdot\left(n_n\mathbf{b}v_{||n} + n_n\mathbf{v}_{\perp n}\right) + S\\ \frac{\partial}{\partial t}\left(n_nv_{||n}\right) =& -\nabla\cdot\left(n_nv_{||n} \mathbf{b}v_{||n} + n_nv_{||n}\mathbf{v}_{\perp n}\right) - \partial_{||}p_n + \nabla_{||}\left(D_{nn}n_n\partial_{||}v_{||n}\right) + F \\ \frac{\partial p_n}{\partial t} =& -\nabla\cdot\left(p_n\mathbf{b}v_{||n} + p_n\mathbf{v}_{\perp n}\right) - \frac{2}{3}p_n\nabla\cdot\left(\mathbf{b}v_{||n}\right) + \nabla\cdot\left(D_{nn}n_n\nabla_\perp T_n\right) + \frac{2}{3}Q \end{aligned}\end{split}\]

The parallel momentum is evolved, so that it can be exchanged with the plasma parallel momentum, but the mass is neglected for perpendicular motion. In the perpendicular direction, therefore, the motion is a balance between the friction (primarily with the plasma through charge exchange) and the pressure gradient:

\[\mathbf{v}_{\perp n} = -D_{nn}\frac{1}{p_n}\nabla_\perp p_n\]

At the moment there is no attempt to limit these velocities, which has been found necessary in UEDGE to get physical results in better agreement with kinetic neutral models [Discussion, T.Rognlien].

Sources

Applying sources using the input file

The simplest way to implement a source in one of the Hermes-3 equations is through the input file. This is done by defining an array representing values of the source across the entire domain using the BOUT++ input file syntax (see BOUT++ documentation).

Sources are available for the density, pressure and momentum equations, and are prescribed under a header corresponding to the chosen equation and species.

For example, this is how a pressure source is prescribed in the 1D-recycling example. First the domain and grid are defined using input file functions. This creates a 400 element 1D grid with a length of 30m and an X-point at the 10m mark. The grid increases in resolution towards the target, with a minimum grid spacing of 0.1 times the average grid spacing:

[mesh]
# 1D simulation, use "y" as the dimension along the fieldline
nx = 1
ny = 400   # Resolution along field-line
nz = 1
length = 30           # Length of the domain in meters
length_xpt = 10   # Length from midplane to X-point [m] (i.e. this is where the source ends)

dymin = 0.1  # Minimum grid spacing near target, as fraction of average. Must be > 0 and < 1

# Parallel grid spacing — grid refinement near the divertor target (which is where the interesting
# stuff happens)
dy = (length / ny) * (1 + (1-dymin)*(1-y/pi))

# Calculate where the source ends in grid index (i.e. at the X-point)
source = length_xpt / length
y_xpt = pi * ( 2 - dymin - sqrt( (2-dymin)^2 - 4*(1-dymin)*source ) ) / (1 - dymin)

And here is how the calculated geometric information is used to prepare a pressure source. First, the required total ion power flux is converted to a pressure according to \(E = 3/2P\), then it is divided by the length of the heating region to obtain the power flux required in each cell. Note that this assumes that \(dx = dz = J = 0\) and that the volume upstream of the X-point is simply an integral of \(dy = mesh:length\_xpt\). If you are imposing a full B-field profile in your 1D simulation, you will need to account for the fact that \(J\) is no longer constant. In order to limit the pressure source to just the region above the X-point, it is multiplied by a Heaviside function which returns 1 upstream of \(y=mesh:y\_xpt\) and 0 downstream of it.

[Pd+]

# Initial condition for ion pressure (in terms of hermes:Nnorm * hermes:Tnorm)
function = 1

# Input power flux to ions in W/m^2
powerflux = 2.5e7

source = (powerflux*2/3 / (mesh:length_xpt))*H(mesh:y_xpt - y)  # Input power as function of y

[Pe]

# Input power flux to electrons in W/m^2
function = `Pd+:function`  # Same as ion pressure initially

source = `Pd+:source`  # Same as ion pressure source

Applying sources using the grid file

The input file has limitations, and sometimes it is useful to prepare an arbitrary profile outside of BOUT++ and import it through the grid file. In 2D, this can be done by adding an appropriate Field3D or Field2D to the grid netCDF file with the sources in the appropriate units.

Time-dependent sources

Any source can be made time-dependent by adding a flag and providing a prefactor function in the input file. The already defined source will be multiplied by the prefactor, which is defined by a time-dependent input file function.

Here is the implementation in the 1D-time-dependent-sources example, where the electrons and ions are set to receive 8MW of mean power flux each with a +/-10% sinusoidal fluctuation of a period of 50us. The density source has a mean of zero and oscillates between \(-1\times10^{22}\) and \(1\times10^{22}\), also with a period of 50us.

Note that if you have the density controller enabled, it will work to counteract the imposed density source oscillation.

[Nd+]
function = 5e19 / hermes:Nnorm # Initial conditions
source_time_dependent = true
source = 1e22 * H(mesh:y_xpt - y)
source_prefactor = sin((2/50)*pi*1e6*t)   #  Oscillation between -1 and 1, period 50us

[Pe]
function = 0.01
powerflux = 16e6  # Input power flux in W/m^2
source = 0.5 * (powerflux*2/3 / (mesh:length_xpt))*H(mesh:y_xpt - y)  # Input power as function of y
source_time_dependent = true
source_prefactor = 1 + 0.1 * sin((2/50)*pi*1e6*t)   #  10% fluctuation on on  top of background source, period 50us

[Pd+]
function = 0.01
source = Pe:source
source_time_dependent = true
source_prefactor = Pe:source_prefactor

Boundary conditions

Simple boundary conditions

BOUT++ simple boundary conditions

BOUT++ provides a number of fundamental boundary conditions including: - dirichlet(x): boundary set to constant value of x - neumann: boundary set to zero gradient - free_o2: boundary set by linear extrapolation (using 2 points) - free_o3: boundary set by quadratic extrapolation (using 3 points)

These can be set on different parts of the domain using the keywords core, sol, pf, lower_target, upper_target, xin, xout, yup, ydown and bndry_all.

The boundary conditions can also be applied over a finite width as well as relaxed over a specified timescale.

These boundary conditions are implemented in BOUT++, and therefore have no access to the normalisations within Hermes-3 and so must be used in normalised units. Please see the BOUT++ documentation for more detail, including the full list of boundary conditions and more guidance on their use. In case the documentation is incomplete or insufficient, please refer to the BOUT++ boundary condition code .

Hermes-3 simple boundary conditions

Currently, there is only one additional simple boundary condition implemented in Hermes-3. decaylength(x) sets the boundary according to a user-set radial decay length. This is a commonly used setting for plasma density and pressure in the tokamak SOL boundary in 2D and 3D but is not applicable in 1D. Note that this must be provided in normalised units just like the BOUT++ simple boundary conditions.

Simple boundary condition examples

The below example for a 2D tokamak simulation sets the electron density to a constant value of 1e20 m:sup:-3 in the core and sets a decay length of 3mm in the SOL and PFR regions, while setting the remaining boundaries to neumann. Example settings of the fundamental normalisation factors and the calculation of the derived ones is provided in the hermes component which can be accessed by using the hermes: prefix in any other component in the input file.

[hermes]
Nnorm = 1e17  # Reference density [m^-3]
Bnorm = 1   # Reference magnetic field [T]
Tnorm = 100   # Reference temperature [eV]
qe = 1.60218e-19   # Electron charge
Mp = 1.67262e-27   # Proton mass
Cs0 = sqrt(qe * Tnorm / Mp)   # Reference speed [m/s]
Omega_ci = qe * Bnorm / Mp   # Reference frequency [1/s]
rho_s0 = Cs0 / Omega_ci   # Refence length [m]

[Ne]
bndry_core = dirichlet(1e20 / hermes:Nnorm)
bndry_sol = decaylength(0.003 / hermes:rho_s0)
bndry_pf = decaylength(0.003 / hermes:rho_s0)
bndry_all = neumann()

Component boundary conditions

Hermes-3 includes additional boundary conditions whose complexity requires their implementation as components. They may overwrite simple boundary conditions and must be set in the same way as other components.

noflow_boundary

This is a species component which imposes a no-flow boundary condition on y (parallel) boundaries.

  • Zero-gradient boundary conditions are applied to density, temperature and pressure fields, if they are set.

  • Zero-value boundary conditions are applied to velocity and momentum if they are set.

By default both yup and ydown boundaries are set, but can be turned off by setting noflow_lower_y or noflow_upper_y to false.

Example: To set no-flow boundary condition on an ion d+ at the lower y boundary, with a sheath boundary at the upper y boundary:

[hermes]
components = d+, sheath_boundary

[d+]
type = noflow_boundary

noflow_lower_y = true   # This is the default
noflow_upper_y = false  # Turn off no-flow at upper y for d+ species

[sheath_boundary]
lower_y = false         # Turn off sheath lower boundary for all species
upper_y = true

Note that currently noflow_boundary is set per-species, whereas sheath_boundary is applied to all species. This is because sheath boundary conditions couple all charged species together, and doesn’t affect neutral species.

The implementation is in NoFlowBoundary:

struct NoFlowBoundary : public Component

Public Functions

virtual void transform(Options &state) override

Inputs

  • species

    • <name>

      • density [Optional]

      • temperature [Optional]

      • pressure [Optional]

      • velocity [Optional]

      • momentum [Optional]

neutral_boundary

Sets Y (sheath/target) boundary conditions on neutral particle density, temperature and pressure. A no-flow boundary condition is set on parallel velocity and momentum. It is a species-specific component and so goes in the list of components for the species that the boundary condition should be applied to.

Just like ions can undergo fast and thermal recycling, neutrals can undergo fast or thermal reflection at the wall. In edge codes using the kinetic neutral code EIRENE, this is typically controlled by the TRIM database. Hermes-3 features a simpler implementation for a constant, user-set fast reflection fraction \(R_{f}\) and energy reflection coefficient \(\alpha_{n}\) based on the approach in the thesis of D.Power 2023.

The two types of reflection are as follows:

  • Fast reflection, where a neutral atom hits the wall and reflects having lost some energy,

  • Thermal reflection, where a neutral atom hits the wall, recombines into a molecule, and then is assumed to immediately dissociate at the Franck Condon dissociation temperature of 3eV.

They are both implemented as a neutral energy sink calculated from the cooling heat flux \(Q_{cool}\):

\[\begin{split}\begin{aligned} Q_{cool} &= Q_{inc} - Q_{fast_refl} - Q_{th_refl} \\ Q_{incident} &= 2n_{n} T_{n} v_{th}^{x} \\ Q_{fast} &= 2n_{n} T_{n} v_{th}^{x} (R_{f} \alpha_{n}) \\ Q_{thermal} &= T_{FC} n_{n} v_{th}^{x} (1 - R_{f}) \\ v_{th}^{x} &= \frac{1}{4}\sqrt{\frac{8k_{B}T_{n}}{\pi m_{n}}} \end{aligned}\end{split}\]

Where \(Q_{incident}\) is the neutral heat flux incident on the wall, \(Q_{fast}\) is the returning heat flux from fast reflection, \(Q_{thermal}\) is the returning heat flux from thermal reflection and \(T_{FC}\) is the Franck-Condon dissociation temperature, currently hardcoded to 3eV. Note that the fast and incident heat flux are both of a Maxwellian distribution, and so their formula corresponds to the 1 dimensional static Maxwellian heat flux and \(v_{th}^{x}\) the corresponding 1D static Maxwellian thermal velocity (Stangeby p.69). The thermal heat flux represents a monoenergetic distribution at \(T_{n}=T_{FC}\) and is therefore calculated with a simpler formula.

Since different regions of the tokamak feature different incidence angles and may feature different materials, the energy reflection coefficient and the fast reflection fraction can be set individually for the target, PFR and SOL walls. The default values are 0.75 for \(\alpha_{n}\) and 0.8 for \(R_{r}\) and correspond to approximate values for tungsten for incidence angles seen at the target. (Power, 2023)

Here are the options set to their defaults. Note that the SOL and PFR are set to have no reflection by default so that it is compatible with a model of any dimensionality which has a target.

[hermes]
components = d

[d]
type = ... , neutral_boundary

neutral_boundary_sol = true
neutral_boundary_pfr = true
neutral_boundary_upper_y = true
neutral_boundary_lower_y = true

target_energy_refl_factor = 0.75
sol_energy_refl_factor = 0.75
pfr_energy_refl_factor = 0.75

target_fast_refl_fraction = 0.80
sol_fast_refl_fraction = 0.80
pfr_fast_refl_fraction = 0.80
struct NeutralBoundary : public Component

Per-species boundary condition for neutral particles at sheath (Y) boundaries.

Sets boundary conditions:

  • Free boundary conditions on logarithm of density, temperature and pressure

  • No-flow boundary conditions on velocity and momentum.

Adds an energy sink corresponding to a flux of heat to the walls.

Heat flux into the wall is q = gamma_heat * n * T * v_th

where v_th = sqrt(eT/m) is the thermal speed

Public Functions

virtual void transform(Options &state) override

state

  • species

    • <name>

      • density Free boundary

      • temperature Free boundary

      • pressure Free boundary

      • velocity [if set] Zero boundary

      • momentum [if set] Zero boundary

      • energy_source Adds wall losses

virtual void outputVars(Options &state) override

Add extra fields for output, or set attributes e.g docstrings.

Others

See sheath_boundary and simple_sheath_boundary.

Collective quantities

These components combine multiple species together. They are typically listed after all the species groups in the component list, so that all the species are present in the state.

One of the most important is the collisions component. This sets collision times for all species, which are then used

sound_speed

Calculates the collective sound speed, by summing the pressure of all species, and dividing by the sum of the mass density of all species:

\[c_s = \sqrt{\sum_i P_i / \sum_i m_in_i}\]

This is set in the state as sound_speed, and is used for the numerical diffusion terms in the parallel advection.

neutral_parallel_diffusion

This adds diffusion to all neutral species (those with no or zero charge), because it needs to be calculated after the collision frequencies are known.

[hermes]
components = ... , collisions, neutral_parallel_diffusion

[neutral_parallel_diffusion]
dneut = 1         # Diffusion multiplication factor
diagnose = true   # This enables diagnostic output for each species

It is intended mainly for 1D simulations, to provide effective parallel diffusion of particles, momentum and energy due to the projection of cross-field diffusion:

\[\begin{split}\begin{aligned} \frac{\partial n_n}{\partial t} =& \ldots + \nabla\cdot\left(\mathbf{b}D_n n_n\partial_{||}p_n\right) \\ \frac{\partial p_n}{\partial t} =& \ldots + \nabla\cdot\left(\mathbf{b}D_n p_n\partial_{||}p_n\right) + \frac{2}{3}\nabla\cdot\left(\mathbf{b}\kappa_n \partial_{||}T_n\right) \\ \frac{\partial}{\partial t}\left(n_nv_{||n}\right) =& \ldots + \nabla\cdot\left(\mathbf{b}D_n n_nv_{||n} \partial_{||}p_n\right) + \nabla\cdot\left(\mathbf{b}\eta_n \partial_{||}T_n\right) \end{aligned}\end{split}\]

The diffusion coefficient is calculated as

\[D_n = \left(\frac{B}{B_{pol}}\right)^2 \frac{T_n}{A \nu}\]

where A is the atomic mass number; \(\nu\) is the collision frequency. The factor \(B / B_{pol}\) is the projection of the cross-field direction on the parallel transport, and is the dneut input setting.

struct NeutralParallelDiffusion : public Component

Add effective diffusion of neutrals in a 1D system, by projecting cross-field diffusion onto parallel distance.

Note: This needs to be calculated after the collision frequency, so is a collective component. This therefore applies diffusion to all neutral species i.e. those with no (or zero) charge

If diagnose = true then the following outputs are saved for each neutral species

  • D<name>_Dpar Parallel diffusion coefficient e.g. Dhe_Dpar

  • S<name>_Dpar Density source due to diffusion

  • E<name>_Dpar Energy source due to diffusion

  • F<name>_Dpar Momentum source due to diffusion

Public Functions

virtual void transform(Options &state) override

Inputs

  • species

    • <all neutrals> # Applies to all neutral species

      • AA

      • collision_frequency

      • density

      • temperature

      • pressure [optional, or density * temperature]

      • velocity [optional]

      • momentum [if velocity set]

Sets

  • species

    • <name>

      • density_source

      • energy_source

      • momentum_source [if velocity set]

virtual void outputVars(Options &state) override

Save variables to the output.

collisions

For collisions between charged particles. In the following all quantities are in SI units except the temperatures: \(T\) is in eV, so \(eT\) has units of Joules.

Debye length \(\lambda_D\)

\[\lambda_D = \sqrt{\frac{\epsilon_0 T_e}{n_e e}}\]

Coulomb logarithm, from [NRL formulary 2019], adapted to SI units

  • For thermal electron-electron collisions

    \[\ln \lambda_{ee} = 30.4 - \frac{1}{2} \ln\left(n_e\right) + \frac{5}{4}\ln\left(T_e\right) - \sqrt{10^{-5} + \left(\ln T_e - 2\right)^2 / 16}\]

    where the coefficient (30.4) differs from the NRL value due to converting density from cgs to SI units (\(30.4 = 23.5 - 0.5\ln\left(10^{-6}\right)\)).

  • Electron-ion collisions

    \[\begin{split}\ln \lambda_{ei} = \left\{\begin{array}{ll} 10 & \textrm{if } T_e < 0.1 \textrm{eV or } n_e < 10^{10}m^{-3} \\ 30 - \frac{1}{2}\ln\left(n_e\right) - \ln(Z) + \frac{3}{2}\ln\left(T_e\right) & \textrm{if } T_im_e/m_i < T_e < 10Z^2 \\ 31 - \frac{1}{2}\ln\left(n_e\right) + \ln\left(T_e\right) & \textrm{if } T_im_e/m_i < 10Z^2 < T_e \\ 23 - \frac{1}{2}\ln\left(n_i\right) + \frac{3}{2}\ln\left(T_i\right) - \ln\left(Z^2\mu\right) & \textrm{if } T_e < T_im_e/m_i \\ \end{array}\right.\end{split}\]

  • Mixed ion-ion collisions

    \[\ln \lambda_{ii'} = 29.91 - ln\left[\frac{ZZ'\left(\mu + \mu'\right)}{\mu T_{i'} + \mu'T_i}\left(\frac{n_iZ^2}{T_i} + \frac{n_{i'} Z'^2}{T_{i'}}\right)^{1/2}\right]\]

    where like the other expressions the different constant is due to converting from cgs to SI units: \(29.91 = 23 - 0.5\ln\left(10^{-6}\right)\).

The frequency of charged species a colliding with charged species b is

\[\nu_{ab} = \frac{1}{3\pi^{3/2}\epsilon_0^2}\frac{Z_a^2 Z_b^2 n_b \ln\Lambda}{\left(v_a^2 + v_b^2\right)^{3/2}}\frac{\left(1 + m_a / m_b\right)}{m_a^2}\]

Note that the cgs expression in Hinton is divided by \(\left(4\pi\epsilon_0\right)^2\) to get the expression in SI units. The thermal speeds in this expression are defined as:

\[v_a^2 = 2 e T_a / m_a\]

Note that with this definition we recover the Braginskii expressions for e-i and i-i collision times.

For conservation of momentum, the collision frequencies \(\nu_{ab}\) and \(\nu_{ba}\) are related by:

\[m_a n_a \nu_{ab} = m_b n_b \nu_{ba}\]

Momentum exchange, force on species a due to collisions with species b:

\[F_{ab} = C_m \nu_{ab} m_a n_a \left( u_b - u_a \right)\]

Where the coefficient \(C_m\) for parallel flows depends on the species: For most combinations of species this is set to 1, but for electron-ion collisions the Braginskii coefficients are used: \(C_m = 0.51\) if ion charge \(Z_i = 1\); 0.44 for \(Z_i = 2\); 0.40 for \(Z_i = 3\); and 0.38 is used for \(Z_i \ge 4\). Note that this coefficient should decline further with increasing ion charge, tending to 0.29 as \(Z_i \rightarrow \infty\).

Frictional heating is included by default, but can be disabled by setting the frictional_heating option to false. When enabled it adds a source of thermal energy corresponding to the resistive heating term:

\[Q_{ab,F} = \frac{m_b}{m_a + m_b} \left( u_b - u_a \right) F_{ab}\]

This term has some important properties:

  1. It is always positive: Collisions of two species with the same temperature never leads to cooling.

  2. It is Galilean invariant: Shifting both species’ velocity by the same amount leaves \(Q_{ab,F}\) unchanged.

  3. If both species have the same mass, the thermal energy change due to slowing down is shared equally between them.

  4. If one species is much heavier than the other, for example electron-ion collisions, the lighter species is preferentially heated. This recovers e.g. Braginskii expressions for \(Q_{ei}\) and \(Q_{ie}\).

This can be derived by considering the exchange of energy \(W_{ab,F}\) between two species at the same temperature but different velocities. If the pressure is evolved then it contains a term that balances the change in kinetic energy due to changes in velocity:

\[\begin{split}\begin{aligned} \frac{\partial}{\partial t}\left(m_a n_a u_a\right) =& \ldots + F_{ab} \\ \frac{\partial}{\partial t}\left(\frac{3}{2}p_a\right) =& \ldots - F_{ab} u_a + W_{ab, F} \end{aligned}\end{split}\]

For momentum and energy conservation we must have \(F_{ab}=-F_{ba}\) and \(W_{ab,F} = -W_{ba,F}\). Comparing the above to the Braginskii expression we see that for ion-electron collisions the term \(- F_{ab}u_a + W_{ab, F}\) goes to zero, so \(W_{ab, F} \sim u_aF_{ab}\) for \(m_a \gg m_b\). An expression that has all these desired properties is

\[W_{ab,F} = \left(\frac{m_a u_a + m_b u_a}{m_a + m_b}\right)F_{ab}\]

which is not Galilean invariant but when combined with the \(- F_{ab} u_a\) term gives a change in pressure that is invariant, as required.

Thermal energy exchange, heat transferred to species \(a\) from species \(b\) due to temperature differences, is given by:

\[Q_{ab,T} = \nu_{ab}\frac{3n_a m_a\left(T_b - T_a\right)}{m_a + m_b}\]

  • Ion-neutral and electron-neutral collisions

    The cross-section for elastic collisions between charged and neutral particles can vary significantly. Here for simplicity we just take a value of \(5\times 10^{-19}m^2\) from the NRL formulary.

  • Neutral-neutral collisions

    The cross-section is given by

\[\sigma = \pi \left(\frac{d_1 + d_2}{2}\right)^2\]

where \(d_1\) and \(d_2\) are the kinetic diameters of the two species. Typical values are [Wikipedia] for H2 2.89e-10m, He 2.60e-10m, Ne 2.75e-10m.

The mean relative velocity of the two species is

\[v_{rel} = \sqrt{\frac{eT_1}{m_1} + \frac{eT_2}{m_2}}\]

and so the collision rate of species 1 on species 2 is:

\[\nu_{12} = v_{rel} n_2 \sigma\]

The implementation is in Collisions:

struct Collisions : public Component

Calculates the collision rate of each species with all other species

Important: Be careful when including both ion_neutral collisions and reactions such as charge exchange, since that may result in double counting. Similarly for electron_neutral collisions and ionization reactions.

Public Functions

Collisions(std::string name, Options &alloptions, Solver*)

The following boolean options under alloptions[name] control which collisions are calculated:

  • electron_electron

  • electron_ion

  • electron_neutral

  • ion_ion

  • ion_neutral

  • neutral_neutral

There are also switches for other terms:

  • frictional_heating Include R dot v heating term as energy source? (includes Ohmic heating)

Parameters:

alloptions – Settings, which should include:

  • units

    • eV

    • inv_meters_cubed

    • meters

    • seconds

virtual void transform(Options &state) override

Modify the given simulation state All components must implement this function

virtual void outputVars(Options &state) override

Add extra fields for output, or set attributes e.g docstrings.

thermal_force

This implements simple expressions for the thermal force. If the electron_ion option is true (which is the default), then a momentum source is added to all ions:

\[F_z = 0.71 n_z Z^2 \nabla_{||}T_e\]

where \(n_z\) is the density of the ions of charge \(Z\). There is an equal and opposite force on the electrons.

If the ion_ion option is true (the default), then forces are calculated between light species (atomic mass < 4) and heavy species (atomic mass > 10). If any combinations of ions are omitted, then a warning will be printed once. The force on the heavy ion is:

\[\begin{split}\begin{aligned} F_z =& \beta \nabla_{||}T_i \\ \beta =& \frac{3\left(\mu + 5\sqrt{2}Z^2\left(1.1\mu^{5/2} - 0.35\mu^{3/2}\right) - 1\right)}{2.6 - 2\mu + 5.4\mu^2} \\ \mu =& m_z / \left(m_z + m_i\right) \end{aligned}\end{split}\]

where subscripts \(z\) refer to the heavy ion, and \(i\) refers to the light ion. The force on the light ion fluid is equal and opposite: \(F_i = -F_z\).

The implementation is in the ThermalForce class:

struct ThermalForce : public Component

Simple calculation of the thermal force

Important: This implements a quite crude approximation, which is intended for initial development and testing. The expressions used are only valid for trace heavy ions and light main ion species, and would not be valid for Helium impurities in a D-T plasma, for example. For this reason only collisions where one ion has an atomic mass < 4, and the other an atomic mass > 10 are considered. Warning messages will be logged for species combinations which are not calculated.

Options used:

  • <name>

    • electron_ion : bool Include electron-ion collisions?

    • ion_ion : bool Include ion-ion elastic collisions?

Public Functions

virtual void transform(Options &state) override

Inputs

  • species

    • e [ if electron_ion true ]

      • charge

      • density

      • temperature

    • <species>

      • charge [ Checks, skips species if not set ]

      • AA

      • temperature [ If AA < 4 i.e. “light” species ]

Outputs

  • species

    • e

      • momentum_source [ if electron_ion true ]

    • <species> [ if AA < 4 (“light”) or AA > 10 (“heavy”) ]

      • momentum_source

recycling

This component calculates the flux of a species into a boundary due to recycling of flow out of the boundary of another species.

The boundary fluxes might be set by sheath boundary conditions, which potentially depend on the density and temperature of all species. Recycling therefore can’t be calculated until all species boundary conditions have been set. It is therefore expected that this component is a top-level component (i.e. in the Hermes section) which comes after boundary conditions are set.

Recycling has been implemented at the target, the SOL edge and the PFR edge. Each is off by default and must be activated with a separate flag. Each can be assigned a separate recycle multiplier and recycle energy.

Configuring thermal recycling

A simple and commonly used way to model recycling is to assume it is fully thermal, i.e. that every incident ion recombines into a neutral molecule and thermalises with the surface before becoming re-emitted. Hermes-3 does not yet have a hydrogenic molecule model, and so the molecules are assumed to instantly dissociate at the Franck-Condon dissociation temperature of 3.5eV.

In order to set this up, the chosen species must feature an outflow through the boundary - any cells with an inflow have their recycling source set to zero. If a sheath boundary condition is enabled, then this is automatically satisfied at the target through the Bohm condition. If it is not enabled, then the target boundary must be set to free_o2, free_o3 or decaylength to allow an outflow.

The recycling component has a species option, that is a list of species to recycle. For each of the species in that list, recycling will look in the corresponding section for the options recycle_as, recycle_multiplier and recycle_energy for each of the three implemented boundaries. Note that the resulting recycling source is a simple multiplication of the outgoing species flow and the multiplier factor. This means that recycling d+ ions into d2 molecules would require a multiplier of 0.5 to maintain a particle balance in the simulation.

For example, recycling d+ ions into d atoms with a recycling fraction of 0.95 at the target and 1.0 at the SOL and PFR edges. Each returning atom has an energy of 3.5eV:

[hermes]
components = d+, d, sheath_boundary, recycling

[recycling]
species = d+   # Comma-separated list of species to recycle

[d+]
recycle_as = d         # Species to recycle as

target_recycle = true
target_recycle_multiplier = 0.95 # Recycling fraction
target_recycle_energy = 3.5   # Energy of recycled particles [eV]

sol_recycle = true
sol_recycle_multiplier = 1 # Recycling fraction
sol_recycle_energy = 3.5   # Energy of recycled particles [eV]

pfr_recycle = true
pfr_recycle_multiplier = 1 # Recycling fraction
pfr_recycle_energy = 3.5   # Energy of recycled particles [eV]

Allowing for fast recycling

In reality, a fraction of incident ions will undergo specular reflection off the surface and preserve a fraction of their energy. In the popular Monte-Carlo neutral code EIRENE, the fast recycling fraction and the energy reflection factor are provided by the TRIM database as a function of incident angle, surface material and incident particle energy. Studies found that sheath acceleration can make the ion angle relatively consistent, e.g. 60 degrees; in (Jae-Sun Park et al 2021 Nucl. Fusion 61 016021).

The recycled heat flux is:

\[\begin{split}\begin{aligned} \Gamma_{E_{n}} &= R \times (R_{f} \alpha_{E} \Gamma_{E_{i}}^{sheath} + (1 - R_{f}) T_{R} \Gamma_{N_{i}})) \\ \end{aligned}\end{split}\]

Where \(R\) is the recycle multiplier, \(R_{f}\) is the fast reflection fraction, \(\alpha_{E}\) is the energy reflection factor, \(\Gamma_{E_{i}}^{sheath}\) is the incident heat flux from the sheath boundary condition, \(T_{R}\) is the recycle energy and \(\Gamma_{N_{i}}\) is the incident ion flux.

\(R_{f}\) and \(\alpha_{E}\) can be set as in the below example. They can also be set to different values for the SOL and PFR by replacing the word “target” with either “sol” or “pfr”.

[d+]
recycle_as = d         # Species to recycle as

target_recycle = true
target_recycle_multiplier = 0.95 # Recycling fraction
target_recycle_energy = 3.5   # Energy of recycled particles [eV]
target_fast_recycle_energy_factor = 0.70
target_fast_recycle_fraction = 0.80

Neutral pump

The recycling component also features a neutral pump which is currently implemented for the SOL and PFR edges only, and so is not available in 1D. The pump is a region of the wall which facilitates particle loss by incomplete recycling and neutral absorption.

The pump requires wall recycling to be enabled on the relevant wall region.

The particle loss rate \(\Gamma_{N_{n}}\) is the sum of the incident ions that are not recycled and the incident neutrals which are not reflected, both of which are controlled by the pump multiplier \(M_{p}\) which is set by the pump_multiplier option in the input file. The unrecycled ion flux \(\Gamma_{N_{i}}^{unrecycled}\) is calculated using the recycling model and allows for either thermal or fast recycling, but with the difference that the pump_multiplier replaces the recycle_multiplier.

\[\begin{split}\begin{aligned} \Gamma_{N_{n}} &= \Gamma_{N_{i}}^{unrecycled} + M_{p} \times \Gamma_{N_{n}}^{incident} \\ \Gamma_{N_{n}}^{incident} &= N_{n} v_{th} = N_{n} \frac{1}{4} \sqrt{\frac{8 T_{n}}{\pi m_{n}}} \\ \end{aligned}\end{split}\]

Where the thermal velocity formulation is for a static maxwellian in 1D (see Stangeby p.64, eqns 2.21, 2.24) and the temperature is in eV.

The heat loss rate \(\Gamma_{E_{n}}\) is calculated as:

\[\begin{split}\begin{aligned} \Gamma_{E_{n}} &= \Gamma_{E_{i}}^{unrecycled} + M_{p} \times \Gamma_{E_{n}}^{incident} \\ \Gamma_{E_{n}}^{incident} &= \gamma T_{n} N_{n} v_{th} = 2 T_{n} N_{n} \frac{1}{4} \sqrt{\frac{8 T_{n}}{\pi m_{n}}} \\ \end{aligned}\end{split}\]

Where the incident heat flux is for a static maxwellian in 1D (see Stangeby p.69, eqn 2.30).

The pump will be placed in any cell that

  1. Is the final domain cell before the guard cells

  2. Is on the SOL or PFR edge

  3. Has a is_pump value of 1

The field is_pump must be created by the user and added to the grid file as a Field2D.

Diagnostic variables

Diagnostic variables for the recycled particle and energy fluxes are provided separately for the targets, the pump as well as the SOL and PFR which are grouped together as wall. as well as the pump. In addition, the field is_pump is saved to help in plotting the pump location.

struct Recycling : public Component

Convert fluxes of species at boundaries

Since this must be calculated after boundary fluxes (e.g. sheath), it is included as a top-level component

Public Functions

Recycling(std::string name, Options &alloptions, Solver*)

Inputs

  • <name>

    • species A comma-separated list of species to recycle

  • <species>

    • recycle_as The species to recycle into

    • recycle_multiplier The recycled flux multiplier, between 0 and 1

    • recycle_energy The energy of the recycled particles [eV]

virtual void transform(Options &state) override

Inputs

  • species

    • <species>

      • density

      • velocity

Outputs

  • species

    • <species>

      • density_source

virtual void outputVars(Options &state) override

Add extra fields for output, or set attributes e.g docstrings.

binormal_stpm

This adds a term to all species which includes the effects of cross-field drifts following the stellarator two point model: Y. Feng et al., Plasma Phys. Control. Fusion 53 (2011) 024009

[hermes]
components = ... , binormal_stpm

[binormal_stpm]
D = 1         # [m^2/s]  Density diffusion coefficient
chi = 3       # [m^2/s]  Thermal diffusion coefficient
nu = 1        # [m^2/s]  Momentum diffusion coefficient

Theta = 1e-3  # Field line pitch

It is intended only for 1D simulations, to provide effective parallel diffusion of particles, momentum and energy due to the projection of cross-field diffusion:

\[\begin{split}\begin{aligned} \frac{\partial N}{\partial t} =& \ldots + \nabla\cdot\left(\mathbf{b}\frac{D}{\Theta}\partial_{||}N\right) \\ \frac{\partial P}{\partial t} =& \ldots + \frac{2}{3}\nabla\cdot\left(\mathbf{b}\frac{\chi}{\Theta} N\partial_{||}T\right) \\ \frac{\partial}{\partial t}\left(NV\right) =& \ldots + \nabla\cdot\left(\mathbf{b}\frac{\nu}{\Theta} \partial_{||}NV\right) \end{aligned}\end{split}\]

The diffusion coefficients D, chi and nu and field line pitch Theta are prescribed in the input file.

struct BinormalSTPM : public Component

Adds terms to Density, Pressure and Momentum equations following the stellarator 2-point model from Yuhe Feng et al., PPCF 53 (2011) 024009 The terms add the effective parallel contribution of perpendicular transport which is of significance in long connection length scenarios. B Shanahan 2023 brendan.shanahan@ipp.mpg.de

Public Functions

BinormalSTPM(std::string name, Options &options, Solver *solver)

Inputs

  • <name>

    • D Perpendicular density diffusion coefficient

    • chi Perpendicular heat diffusion coefficient

    • nu Perpendicular momentum diffusion coefficient

    • Theta Field line pitch as described by Feng et al.

virtual void transform(Options &state) override

Sets

  • species

    • <name>

      • pressure correction

      • momentum correction

      • density correction

virtual void outputVars(Options &state) override

Add extra fields for output, or set attributes e.g docstrings.

Atomic and molecular reactions

The formula for the reaction is used as the name of the component. This makes writing the input file harder, since the formula must be in the exact same format (e.g. h + e and e + h won’t be recognised as being the same thing), but makes reading and understanding the file easier.

To include a set of reactions, it is probably easiest to group them, and then include the group name in the components list

[hermes]
components = ..., reactions

[reactions]
type = (
        h + e -> h+ + 2e,  # ionisation
        h+ + e -> h,    # Radiative + 3-body recombination
       )

Note that brackets can be used to split the list of reactions over multiple lines, and trailing commas are ignored. Comments can be used if needed to add explanation. The name of the section does not need to be reactions, and multiple components could be created with different reaction sets. Be careful not to include the same reaction twice.

When reactions are added, all the species involved must be included, or an exception should be thrown.

Diagnostic variables

Diagnostic variables are provided for each reaction channel of density, momentum and energy transfer. Additionally, charge exchange features a diagnostic for the reaction rate (in ionisation and recombination, the reaction rate K is simply the density transfer rate S divided by the ion density). The sign convention is always in terms of a plasma source, so that a source of plasma density, energy or momentum is positive, and a sink is negative. Radiative energy transfer is provided separately as E is a transfer of energy between two species, while R is a net loss of energy from the system due to the plasma being transparent.

Variable prefix

Units

Description

K

\(s^{-1}\)

Reaction rate

S

\(m^{-3}s^{-1}\)

Density transfer rate

E

\(Wm^{-3}\)

Energy transfer rate

R

\(Wm^{-3}\)

Radiation

Notes:

  1. Charge exchange channel diagnostics: For two species a and b, the channel Fab_cx is a source of momentum for species a due to charge exchange with species b. There are corresponding sinks for the products of the charge exchange reaction which are not saved.

    For example,reaction d + t+ -> d+ + t will save the following forces (momentum sources): - Fdt+_cx is a source of momentum for deuterium atoms d and sink of momentum for deuterium ions d+. - Ft+d_cx is a source of momentum for tritium ions t+ and sink of momentum for tritium atoms t

    The reason for this convention is the existence of the inverse reactions: t + d+ -> t+ + d outputs diagnostics Ftd+_cx and Fd+t_cx.

  2. Reactions typically convert species from one to another, leading to a transfer of mass momentum and energy. For a reaction converting species \(a\) to species \(b\) at rate \(R\) (units of events per second per volume) we have transfers:

    \[\begin{split}\begin{aligned} \frac{\partial}{\partial t} n_a =& \ldots - R \\ \frac{\partial}{\partial t} n_b =& \ldots + R \\ \frac{\partial}{\partial t}\left( m n_a u_a\right) =& \ldots + F_{ab} \\ \frac{\partial}{\partial t}\left( m n_a u_a\right) =& \ldots + F_{ba} \\ \frac{\partial}{\partial t}\left( \frac{3}{2} p_a \right) =& \ldots - F_{ab}u_a + W_{ab} - \frac{1}{2}mRu_a^2 \\ \frac{\partial}{\partial t}\left( \frac{3}{2} p_b \right) =& \ldots - F_{ba}u_b + W_{ba} + \frac{1}{2}mRu_b^2 \end{aligned}\end{split}\]

where both species have the same mass: \(m_a = m_b = m\). In the pressure equations the \(-F_{ab}u_a\) comes from splitting the kinetic and thermal energies; \(W_{ab}=-W_{ba}\) is the energy transfer term that we need to find; The final term balances the loss of kinetic energy at fixed momentum due to a particle source or sink.

The momentum transfer \(F_{ab}=-F{ba}\) is the momentum carried by the converted ions: \(F_{ab}=-m R u_a\). To find \(W_{ab}\) we note that for \(p_a = 0\) the change in pressure must go to zero: \(-F_{ab}u_a + W_{ab} -\frac{1}{2}mRu_a^2 = 0\).

\[\begin{split}\begin{aligned} W_{ab} =& F_{ab}u_a + \frac{1}{2}mRu_a^2 \\ =& - mR u_a^2 + \frac{1}{2}mRu_a^2\\ =& -\frac{1}{2}mRu_a^2 \end{aligned}\end{split}\]

Substituting into the above gives:

\[\begin{split}\begin{aligned} \frac{\partial}{\partial t}\left( \frac{3}{2} p_b \right) =& \ldots - F_{ba}u_b + W_{ba} + \frac{1}{2}mRu_b^2 \\ =& \ldots - mRu_au_b + \frac{1}{2}mRu_a^2 + \frac{1}{2}mRu_a^2 \\ =& \ldots + \frac{1}{2}mR\left(u_a - u_b\right)^2 \end{aligned}\end{split}\]

This has the property that the change in pressure of both species is Galilean invariant. This transfer term is included in the Amjuel reactions and hydrogen charge exchange.

Hydrogen

Multiple isotopes of hydrogen can be evolved, so to keep track of this the species labels h, d and t are all handled by the same hydrogen atomic rates calculation. The following might therefore be used

[hermes]
components = d, t, reactions

[reactions]
type = (
        d + e -> d+ + 2e,  # Deuterium ionisation
        t + e -> t+ + 2e,  # Tritium ionisation
       )

Reaction

Description

h + e -> h+ + 2e

Hydrogen ionisation (Amjuel H.4 2.1.5)

d + e -> d+ + 2e

Deuterium ionisation (Amjuel H.4 2.1.5)

t + e -> t+ + 2e

Tritium ionisation (Amjuel H.4 2.1.5)

h + h+ -> h+ + h

Hydrogen charge exchange (Amjuel H.3 3.1.8)

d + d+ -> d+ + d

Deuterium charge exchange (Amjuel H.3 3.1.8)

t + t+ -> t+ + t

Tritium charge exchange (Amjuel H.3 3.1.8)

h + d+ -> h+ + d

Mixed hydrogen isotope CX (Amjuel H.3 3.1.8)

d + h+ -> d+ + h

h + t+ -> h+ + t

t + h+ -> t+ + h

d + t+ -> d+ + t

t + d+ -> t+ + d

h+ + e -> h

Hydrogen recombination (Amjuel H.4 2.1.8)

d+ + e -> d

Deuterium recombination (Amjuel H.4 2.1.8)

t+ + e -> t

Tritium recombination (Amjuel H.4 2.1.8)

In addition, the energy loss associated with the ionisation potential energy cost as well as the photon emission during excitation and de-excitation during multi-step ionisation is calculated using the AMJUEL rate H.10 2.1.5. The equivalent rate for recombination is H.10 2.1.8.

The code to calculate the charge exchange rates is in hydrogen_charge_exchange.[ch]xx. This implements reaction H.3 3.1.8 from Amjuel (p43), scaled to different isotope masses and finite neutral particle temperatures by using the effective temperature (Amjuel p43):

\[T_{eff} = \frac{M}{M_1}T_1 + \frac{M}{M_2}T_2\]

The effective hydrogenic ionisation rates are calculated using Amjuel reaction H.4 2.1.5, by D.Reiter, K.Sawada and T.Fujimoto (2016). Effective recombination rates, which combine radiative and 3-body contributions, are calculated using Amjuel reaction 2.1.8.

struct HydrogenChargeExchange : public Component

Hydrogen charge exchange total rate coefficient

p + H(1s) -> H(1s) + p

Reaction 3.1.8 from Amjuel (p43)

Scaled to different isotope masses and finite neutral particle temperatures by using the effective temperature (Amjuel p43)

T_eff = (M/M_1)T_1 + (M/M_2)T_2

Important: If this is included then ion_neutral collisions should probably be disabled in the collisions component, to avoid double-counting.

Subclassed by HydrogenChargeExchangeIsotope< Isotope1, Isotope2 >

Public Functions

inline HydrogenChargeExchange(std::string name, Options &alloptions, Solver*)
Parameters:

alloptions – Settings, which should include:

  • units

    • eV

    • inv_meters_cubed

    • seconds

Helium

Reaction

Description

he + e -> he+ + 2e

He ionisation, unresolved metastables (Amjuel 2.3.9a)

he+ + e -> he

He+ recombination, unresolved metastables (Amjuel 2.3.13a)

The implementation of these rates are in the AmjuelHeIonisation01 and AmjuelHeRecombination10 classes:

struct AmjuelHeIonisation01 : public AmjuelReaction

e + he -> he+ + 2e Amjuel reaction 2.3.9a, page 161 Not resolving metastables, only transporting ground state

Public Functions

inline virtual void transform(Options &state) override

Modify the given simulation state All components must implement this function

struct AmjuelHeRecombination10 : public AmjuelReaction

e + he+ -> he Amjuel reaction 2.3.13a Not resolving metastables. Includes radiative + threebody + dielectronic. Fujimoto Formulation II

Public Functions

inline virtual void transform(Options &state) override

Modify the given simulation state All components must implement this function

Lithium

These rates are taken from ADAS (‘96 and ‘89)

Reaction

Description

li + e -> li+ + 2e

Lithium ionisation

li+ + e -> li+2 + 2e

li+2 + e -> li+3 + 2e

li+ + e -> li

Lithium recombination

li+2 + e -> li+

li+3 + e -> li+2

li+ + h -> li + h+

Charge exchange with hydrogen

li+2 + h -> li+ + h+

li+3 + h -> li+2 + h+

li+ + d -> li + d+

Charge exchange with deuterium

li+2 + d -> li+ + d+

li+3 + d -> li+2 + d+

li+ + t -> li + t+

Charge exchange with tritium

li+2 + t -> li+ + t+

li+3 + t -> li+2 + t+

The implementation of these rates is in ADASLithiumIonisation, ADASLithiumRecombination and ADASLithiumCX template classes:

template<int level>
struct ADASLithiumIonisation : public OpenADAS

ADAS effective ionisation (ADF11)

Template Parameters:

level – The ionisation level of the ion on the left of the reaction

Public Functions

inline virtual void transform(Options &state) override

Modify the given simulation state All components must implement this function

template<int level>
struct ADASLithiumRecombination : public OpenADAS

ADAS effective recombination coefficients (ADF11)

Template Parameters:

level – The ionisation level of the ion on the right of the reaction

Public Functions

inline ADASLithiumRecombination(std::string, Options &alloptions, Solver*)
Parameters:

alloptions – The top-level options. Only uses the [“units”] subsection.

inline virtual void transform(Options &state) override

Modify the given simulation state All components must implement this function

template<int level, char Hisotope>
struct ADASLithiumCX : public OpenADASChargeExchange
Template Parameters:

  • level – The ionisation level of the ion on the right of the reaction

  • Hisotope – The hydrogen isotope (‘h’, ‘d’ or ‘t’)

Public Functions

inline ADASLithiumCX(std::string, Options &alloptions, Solver*)
Parameters:

alloptions – The top-level options. Only uses the [“units”] subsection.

inline virtual void transform(Options &state) override

Modify the given simulation state All components must implement this function

Neon

These rates are taken from ADAS (96): SCD and PLT are used for the ionisation rate and radiation energy loss; ACD and PRB for the recombination rate and radiation energy loss; and CCD (89) for the charge exchange coupling to hydrogen. The ionisation potential is also included as a source or sink of energy for the electrons.

Reaction

Description

ne + e -> ne+ + 2e

Neon ionisation

ne+ + e -> ne+2 + 2e

ne+2 + e -> ne+3 + 2e

ne+3 + e -> ne+4 + 2e

ne+4 + e -> ne+5 + 2e

ne+5 + e -> ne+6 + 2e

ne+6 + e -> ne+7 + 2e

ne+7 + e -> ne+8 + 2e

ne+8 + e -> ne+9 + 2e

ne+9 + e -> ne+10 + 2e

ne+ + e -> ne

Neon recombination

ne+2 + e -> ne+

ne+3 + e -> ne+2

ne+4 + e -> ne+3

ne+5 + e -> ne+4

ne+6 + e -> ne+5

ne+7 + e -> ne+6

ne+8 + e -> ne+7

ne+9 + e -> ne+8

ne+10 + e -> ne+9

ne+ + h -> ne + h+

Charge exchange with hydrogen

ne+2 + h -> ne+ + h+

ne+3 + h -> ne+2 + h+

ne+4 + h -> ne+3 + h+

ne+5 + h -> ne+4 + h+

ne+6 + h -> ne+5 + h+

ne+7 + h -> ne+6 + h+

ne+8 + h -> ne+7 + h+

ne+9 + h -> ne+8 + h+

ne+10 + h -> ne+9 + h+

ne+ + d -> ne + d+

Charge exchange with deuterium

ne+2 + d -> ne+ + d+

ne+3 + d -> ne+2 + d+

ne+4 + d -> ne+3 + d+

ne+5 + d -> ne+4 + d+

ne+6 + d -> ne+5 + d+

ne+7 + d -> ne+6 + d+

ne+8 + d -> ne+7 + d+

ne+9 + d -> ne+8 + d+

ne+10 + d -> ne+9 + d+

ne+ + t -> ne + t+

Charge exchange with tritium

ne+2 + t -> ne+ + t+

ne+3 + t -> ne+2 + t+

ne+4 + t -> ne+3 + t+

ne+5 + t -> ne+4 + t+

ne+6 + t -> ne+5 + t+

ne+7 + t -> ne+6 + t+

ne+8 + t -> ne+7 + t+

ne+9 + t -> ne+8 + t+

ne+10 + t -> ne+9 + t+

The implementation of these rates is in ADASNeonIonisation, ADASNeonRecombination and ADASNeonCX template classes:

template<int level>
struct ADASNeonIonisation : public OpenADAS

ADAS effective ionisation (ADF11)

Template Parameters:

level – The ionisation level of the ion on the left of the reaction

Public Functions

inline virtual void transform(Options &state) override

Modify the given simulation state All components must implement this function

template<int level>
struct ADASNeonRecombination : public OpenADAS

ADAS effective recombination coefficients (ADF11)

Template Parameters:

level – The ionisation level of the ion on the right of the reaction

Public Functions

inline ADASNeonRecombination(std::string, Options &alloptions, Solver*)
Parameters:

alloptions – The top-level options. Only uses the [“units”] subsection.

inline virtual void transform(Options &state) override

Modify the given simulation state All components must implement this function

template<int level, char Hisotope>
struct ADASNeonCX : public OpenADASChargeExchange
Template Parameters:

  • level – The ionisation level of the ion on the right of the reaction

  • Hisotope – The hydrogen isotope (‘h’, ‘d’ or ‘t’)

Public Functions

inline ADASNeonCX(std::string, Options &alloptions, Solver*)
Parameters:

alloptions – The top-level options. Only uses the [“units”] subsection.

inline virtual void transform(Options &state) override

Modify the given simulation state All components must implement this function

Fixed fraction radiation

These components produce volumetric electron energy losses, but don’t otherwise modify the plasma solution: Their charge and mass density are not calculated, and there are no interactions with other species or boundary conditions.

The fixed_fraction_hutchinson_carbon component calculates radiation due to carbon in coronal equilibrium, using a simple formula from I.H.Hutchinson Nucl. Fusion 34 (10) 1337 - 1348 (1994):

\[L\left(T_e\right) = 2\times 10^{-31} \frac{\left(T_e/10\right)^3}{1 + \left(T_e / 10\right)^{4.5}}\]

which has units of \(Wm^3\) with \(T_e\) in eV.

To use this component you can just add it to the list of components and then configure the impurity fraction:

[hermes]
components = ..., fixed_fraction_hutchinson_carbon, ...

[fixed_fraction_hutchinson_carbon]
fraction = 0.05   # 5% of electron density
diagnose = true   # Saves Rfixed_fraction_carbon to output

Or to customise the name of the radiation output diagnostic a section can be defined like this:

[hermes]
components = ..., c, ...

[c]
type = fixed_fraction_hutchinson_carbon
fraction = 0.05   # 5% of electron density
diagnose = true   # Saves Rc (R + section name)

Carbon is also provided as an ADAS rate along with nitrogen, neon and argon. The component names are fixed_fraction_carbon, fixed_fraction_nitrogen, fixed_fraction_neon and fixed_fraction_argon.

These can be used in the same way as fixed_fraction_hutchinson_carbon. Each rate is in the form of a 10 coefficient log-log polynomial fit of data obtained using the open source tool radas. The \(n {\tau}\) parameter representing the density and residence time assumed in the radas collisional-radiative model has been set to \(1\times 10^{20} \times 0.5ms\) based on David Moulton et al 2017 Plasma Phys. Control. Fusion 59(6).

Each rate has an upper and lower bound beyond which the rate remains constant. Please refer to the source code in fixed_fraction_radiation.hxx for the coefficients and bounds used for each rate.

In addition to the above rates, there are three simplified cooling curves for Argon: fixed_fraction_argon_simplified1, fixed_fraction_argon_simplified2 and fixed_fraction_argon_simplified3. They progressively reduce the nonlinearity in the rate by taking out the curvature from the slopes, taking out the RHS shoulder and taking out the LHS-RHS asymmetry, respectively. These rates may be useful in investigating the impact of the different kinds of curve nonlinearities on the solution.

Adjusting reactions

The reaction rates can be adjusted by a user-specified arbitrary multiplier. This can be useful for the analysis of the impact of individual reactions. The multiplier setting must be placed under the neutral species corresponding to the reaction, e.g. under [d] when adjusting deuterium ionisation, recombination or charge exchange. The multiplier for the fixed fraction impurity radiation must be placed under the impurity species header, e.g. under [ar] for argon. This functionality is not yet currently implemented for helium or neon reactions.

Setting

Specified under

Reaction

K_iz_multiplier

Neutral species

Ionisation rate

R_ex_multiplier

Neutral species

Ionisation (excitation) radiation rate

K_rec_multiplier

Neutral species

Recombination rate

R_rec_multiplier

Neutral species

Recombination radiation rate

K_cx_multiplier

Neutral species

Charge exchange rate

R_multiplier

Impurity species

Fixed frac. impurity radiation rate

Electromagnetic fields

These are components which calculate the electric and/or magnetic fields.

vorticity

Evolves a vorticity equation, and at each call to transform() uses a matrix inversion to calculate potential from vorticity.

In this component the Boussinesq approximation is made, so the vorticity equation solved is

\[\nabla\cdot\left(\frac{\overline{A}\overline{n}}{B^2}\nabla_\perp \phi\right) \underbrace{+ \nabla\cdot\left(\sum_i\frac{A_i}{Z_i B^2}\nabla_\perp p_i\right)}_{\mathrm{if diamagnetic\_polarisation}} = \Omega\]

Where the sum is over species, \(\overline{A}\) is the average ion atomic number, and \(\overline{n}\) is the normalisation density (i.e. goes to 1 in the normalised equations). The ion diamagnetic flow terms in this Boussinesq approximation can be written in terms of an effective ion pressure \(\hat{p}\):

\[\hat{p} \equiv \sum_i \frac{A_i}{\overline{A} Z_i} p_i\]

as

\[\nabla\cdot\left[\frac{\overline{A}\overline{n}}{B^2}\nabla_\perp \left(\phi + \frac{\hat{p}}{\overline{n}}\right) \right] = \Omega\]

Note that if diamagnetic_polarisation = false then the ion pressure terms are removed from the vorticity, and also from other ion pressure terms coming from the polarisation current (i.e. \(\hat{p}\rightarrow 0\).

This is a simplified version of the full vorticity definition which is:

\[\nabla\cdot\left(\sum_i \frac{A_i n_i}{B^2}\nabla_\perp \phi + \sum_i \frac{A_i}{Z_i B^2}\nabla_\perp p_i\right) = \Omega\]

and is derived by replacing

\[\sum_i A_i n_i \rightarrow \overline{A}\overline{n}\]

In the case of multiple species, this Boussinesq approximation means that the ion diamagnetic flow terms

The vorticity equation that is integrated in time is

\[\begin{split}\begin{aligned}\frac{\partial \Omega}{\partial t} =& \nabla\cdot\left(\mathbf{b}\sum_s Z_s n_sV_{||s}\right) \\ &+ \underbrace{\nabla\cdot\left(\nabla\times\frac{\mathbf{b}}{B}\sum_s p_s\right)}_{\textrm{if diamagnetic}} + \underbrace{\nabla\cdot\mathbf{J_{exb}}}_{\mathrm{if exb\_advection}} \\ &+ \nabla\cdot\left(\mathbf{b}J_{extra}\right)\end{aligned}\end{split}\]

The nonlinearity \(\nabla\cdot\mathbf{J_{exb}}\) is part of the divergence of polarisation current. In its simplified form when exb_advection_simplified = true, this is the \(E\times B\) advection of vorticity:

\[\nabla\cdot\mathbf{J_{exb}} = -\nabla\cdot\left(\Omega \mathbf{V}_{E\times B}\right)\]

When exb_advection_simplified = false then the more complete (Boussinesq approximation) form is used:

\[\nabla\cdot\mathbf{J_{exb}} = -\nabla\cdot\left[\frac{\overline{A}}{2B^2}\nabla_\perp\left(\mathbf{V}_{E\times B}\cdot\nabla \hat{p}\right) + \frac{\Omega}{2} \mathbf{V}_{E\times B} + \frac{\overline{A}\overline{n}}{2B^2}\nabla_\perp^2\phi\left(\mathbf{V}_{E\times B} + \frac{\mathbf{b}}{B}\times\nabla\hat{p}\right) \right]\]

The form of the vorticity equation is based on Simakov & Catto (corrected in erratum 2004), in the Boussinesq limit and with the first term modified to conserve energy. In the limit of zero ion pressure and constant \(B\) it reduces to the simplified form.

struct Vorticity : public Component

Evolve electron density in time

Public Functions

Vorticity(std::string name, Options &options, Solver *solver)

Options

  • <name>

    • average_atomic_mass: float, default 2.0 Weighted average ion atomic mass for polarisation current

    • bndry_flux: bool, default true Allow flows through radial (X) boundaries?

    • collisional_friction: bool, default false Damp vorticity based on mass-weighted collision frequency?

    • diagnose: bool, false Output additional diagnostics?

    • diamagnetic: bool, default true Include diamagnetic current, using curvature vector?

    • diamagnetic_polarisation: bool, default true Include ion diamagnetic drift in polarisation current?

    • exb_advection: bool, default true Include ExB advection (nonlinear term)?

    • hyper_z: float, default -1.0 Hyper-viscosity in Z. < 0 means off

    • laplacian: subsection Options for the Laplacian phi solver

    • phi_boundary_relax: bool, default false Relax radial phi boundaries towards zero-gradient?

    • phi_boundary_timescale: float, 1e-4 Timescale for phi boundary relaxation [seconds]

    • phi_dissipation: bool, default true Parallel dissipation of potential (Recommended)

    • poloidal_flows: bool, default true Include poloidal ExB flow?

    • sheath_boundary: bool, default false If phi_boundary_relax is false, set the radial boundary to the sheath potential?

    • split_n0: bool, default false Split phi into n=0 and n!=0 components?

    • viscosity: Field2D, default 0.0 Kinematic viscosity [m^2/s]

    • vort_dissipation: bool, default false Parallel dissipation of vorticity?

    • damp_core_vorticity: bool, default false Damp axisymmetric component of vorticity in cell next to core boundary

virtual void transform(Options &state) override

Optional inputs

  • species

    • pressure and charge => Calculates diamagnetic terms [if diamagnetic=true]

    • pressure, charge and mass => Calculates polarisation current terms [if diamagnetic_polarisation=true]

Sets in the state

  • species

    • [if has pressure and charge]

      • energy_source

  • fields

    • vorticity

    • phi Electrostatic potential

    • DivJdia Divergence of diamagnetic current [if diamagnetic=true]

Note: Diamagnetic current calculated here, but could be moved to a component with the diamagnetic drift advection terms

virtual void finally(const Options &state) override

Optional inputs

  • fields

    • DivJextra Divergence of current, including parallel current Not including diamagnetic or polarisation currents

virtual void outputVars(Options &state) override

Add extra fields for output, or set attributes e.g docstrings.

inline virtual void restartVars(Options &state) override

Add extra fields to restart files.

relax_potential

This component evolves a vorticity equation, similar to the vorticity component. Rather than inverting an elliptic equation at every timestep, this component evolves the potential in time as a diffusion equation.

struct RelaxPotential : public Component

Evolve vorticity and potential in time.

Uses a relaxation method for the potential, which is valid for steady state, but not for timescales shorter than the relaxation timescale.

Public Functions

RelaxPotential(std::string name, Options &options, Solver *solver)

Options

  • <name>

    • diamagnetic

    • diamagnetic_polarisation

    • average_atomic_mass

    • bndry_flux

    • poloidal_flows

    • split_n0

    • laplacian Options for the Laplacian phi solver

virtual void transform(Options &state) override

Optional inputs

  • species

    • pressure and charge => Calculates diamagnetic terms [if diamagnetic=true]

    • pressure, charge and mass => Calculates polarisation current terms [if diamagnetic_polarisation=true]

Sets in the state

  • species

    • [if has pressure and charge]

      • energy_source

  • fields

    • vorticity

    • phi Electrostatic potential

    • DivJdia Divergence of diamagnetic current [if diamagnetic=true]

Note: Diamagnetic current calculated here, but could be moved to a component with the diamagnetic drift advection terms

virtual void finally(const Options &state) override

Optional inputs

  • fields

    • DivJextra Divergence of current, including parallel current Not including diamagnetic or polarisation currents

virtual void outputVars(Options &state) override

Add extra fields for output, or set attributes e.g docstrings.

electromagnetic

This component modifies the definition of momentum of all species, to include the contribution from the electromagnetic potential \(A_{||}\).

Assumes that “momentum” \(p_s\) calculated for all species \(s\) is

\[p_s = m_s n_s v_{||s} + Z_s e n_s A_{||}\]

which arises once the electromagnetic contribution to the force on each species is included in the momentum equation. This is normalised so that in dimensionless quantities

\[p_s = A n v_{||} + Z n A_{||}\]

where \(A\) and \(Z\) are the atomic number and charge of the species.

The current density \(j_{||}\) in SI units is

\[j_{||} = -\frac{1}{\mu_0}\nabla_\perp^2 A_{||}\]

which when normalised in Bohm units becomes

\[j_{||} = - \frac{1}{\beta_{em}}\nabla_\perp^2 A_{||}\]

where \(\beta_{em}\) is a normalisation parameter which is half the plasma electron beta as normally defined:

\[\beta_{em} = \frac{\mu_0 e \overline{n} \overline{T}}{\overline{B}^2}\]

To convert the species momenta into a current, we take the sum of \(p_s Z_s e / m_s\). In terms of normalised quantities this gives:

\[- \frac{1}{\beta_{em}} \nabla_\perp^2 A_{||} + \sum_s \frac{Z^2 n_s}{A}A_{||} = \sum_s \frac{Z}{A} p_s\]
struct Electromagnetic : public Component

Electromagnetic potential A||

Reinterprets all species’ parallel momentum as a combination of a parallel flow and a magnetic contribution, i.e. canonical momentum. 

m n v_{||} + Z e n A_{||}
Changes the “momentum” of each species so that after this component the momentuum of each species is just 
m n v_{||}
This component should be run after all species have set their momentum, but before the momentum is used e.g to set boundary conditions.

Calculates the electromagnetic potential A_{||} using

Laplace(Apar) - alpha_em * Apar = -Ajpar

By default outputs Apar every timestep. When diagnose = true in also saves alpha_em and Ajpar.

Public Functions

Electromagnetic(std::string name, Options &options, Solver *solver)

Options

  • units

  • <name>

    • diagnose Saves Ajpar and alpha_em time-dependent values

virtual void transform(Options &state) override

Inputs

  • species

    • <..> All species with charge and parallel momentum

      • charge

      • momentum

      • density

      • AA

Sets

  • species

    • <..> All species with charge and parallel momentum

      • momentum (modifies) to m n v||

      • velocity (modifies) to v||

  • fields

virtual void outputVars(Options &state) override

Add extra fields for output, or set attributes e.g docstrings.