Equations

Contents

Equations#

This section contains components which determine which equations are being solved in the code. There are two broad classes of components:

Whole equations

For example, fixed_temperature, evolve_pressure, evolve_energy allow the solution of energy in three levels of fidelity: constant temperature, a pressure equation and the conservative total energy equation. neutral_mixed contains both parallel and perpendicular transport of neutrals and has several equations included within.

Terms

For example, anomalous_diffusion adds cross-field transport to the density, energy and momentum equations if they are available while diamagnetic_drift and polarisation_drift add drift terms.

Please refer to the examples for common component configurations.

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 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 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.

fixed_fraction_ions#

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

struct FixedFractionIons : public Component

Set ion densities from electron densities

Public Functions

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

Inputs

  • <name>

    • fractions A comma-separated list of pairs separated by @ e.g. ‘d+ @ 0.5, t+ @ 0.5’

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.

struct Quasineutral : public Component

Calculate density from sum of other species densities * charge to ensure that net charge = 0

This is useful in simulations where multiple species are being evolved. Note that only one species’ density can be calculated this way, and it should be calculated last once all other densities are known.

Saves the density to the output (dump) files as N<name>

Public Functions

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

Inputs

Parameters:

  • name – Short name for species e.g. “e”

  • alloptions – Component configuration options

    • <name>

      • charge Required to have a particle charge

      • AA Atomic mass

virtual void finally(const Options &state) override

Get the final density for output including any boundary conditions applied

virtual void outputVars(Options &state) override

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

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 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 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<species> where <species> is the short name of the evolving species e.g. Pe.

Parallel conduction is included if the global braginskii_conduction component has been used.

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}S_E + \frac{2}{3}S_N\frac{1}{2}mNV^2\]

where \(S_E\) is the energy_source (thermal energy source), and \(S_N\) is the density source. If conduction has been calculated, it will be included in the energy_source term.

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

    • 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 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<species> (e.g Ed) and the pressure as P<species> (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

    • poloidal_flows Include poloidal ExB flows? Default is true

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

  • 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 finally(const Options &state) override

Optional inputs

  • species

    • <name>

      • velocity. Must have sound_speed or temperature

      • energy_source

  • 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 outputVars(Options &state) override

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

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/node55.html

braginskii_conduction#

This is a global component that calculates the parallel thermal conduction for all species that use evolve_pressure or evolve_energy, storing it in energy_source. If this is not desired for a particular species then it can be turned off by setting thermal_conduction = false in the input options for that species.

This component requires a collision time to have been calculated (i.e., with the Braginskii Collisions component component). It is recommended that this be one of the last component to run, to ensure density, pressure, and temperature have their final values. However, it must be run before Recycling, as that component will need to use the energy_flow_ylow value, to which conduction contributes.

The choice of collision frequency used for conduction is set by the flag conduction_collisions_mode: multispecies uses all available collision frequencies involving the chosen species, while braginskii uses only self-collisions .The default is multispecies and it is recommended for use if solving more than one ion. If you are solving for a single ion and want to recover Braginskii, use the braginskii mode.

\[\nabla\cdot\left(\kappa_{||}\mathbf{b}\mathbf{b}\cdot\nabla T\right)\]
struct BraginskiiConduction : public Component

Calculates parallel heat conduction due to collisions

NOTE: This is global as that is the only way to ensure it gets run after all collisions have been calculated. Logically this should really apply species-by-species, but we can’t do that until we have dynamic ordering to make sure collision rates get calculated first.

Public Functions

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

Inputs

  • <component name>

    • diagnose Output additional diagnostic fields?

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

    • kappa_limit_alpha Flux limiter, off by default.

    • conduction_collisions_mode Can be multispecies: all collisions, or braginskii: self collisions and ie

  • <species name>

    • type Checks whether energy or pressure are evolved

    • thermal_conduction Include parallel heat conduction? Default is true

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.

The velocity may be set in the mesh file as an array (2D or 3D), or in the options as an expression. The options value overrides the value in the mesh. If neither mesh array nor option are set then an exception will be thrown. Both mesh array and option should be specified in units of meters per second.

The name of the array in the mesh file is V<name>0 where <name> is the name of the species e.g. for species e (electrons), fixed_velocity will try to read Ve0 from the mesh file, and then the velocity option in the [e] section:

[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 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 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
struct ZeroCurrent : public Component

Set the velocity of a species so that there is no net current, by summing the current from other species.

This is most often used in the electron species, but does not need to be.

Public Functions

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

Inputs

Parameters:

  • name – Short name for species e.g. “e”

  • alloptions – Component configuration options

    • <name>

      • charge (must not be zero)

inline virtual void finally(const Options &state) override

Use the final simulation state to update internal state (e.g. time derivatives)

virtual void outputVars(Options &state) override

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

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 outputVars(Options &state) override

Save output diagnostics.

braginskii_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, ..., braginskii_collisions, braginskii_electron_viscosity
struct BraginskiiElectronViscosity : 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

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

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 outputVars(Options &state) override

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

braginskii_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 =  ..., braginskii_collisions, braginskii_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\]

The choice of collision frequency is set by the flag viscosity_collisions_mode: multispecies uses all available collision frequencies involving the chosen species, while braginskii uses only ii collisions. The default is multispecies and it is recommended when solving more than one ion. If you are solving for a single ion and want to recover Braginskii, use the braginskii mode.

If the perpendicular option is set:

[braginskii_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.61\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\]

Note

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 BraginskiiIonViscosity : 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

BraginskiiIonViscosity(const 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 outputVars(Options &state) override

Save variables to the output.

braginskii_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 BraginskiiThermalForce class:

struct BraginskiiThermalForce : public Component

Simple calculation of the thermal force for the Braginskii closure

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. The restrictions can be overridden by setting the override_ion_mass_restrictions=true.

Options used:

  • <name>

    • electron_ion : bool Include electron-ion collisions?

    • ion_ion : bool Include ion-ion elastic collisions?

Neutral gas models#

In 1D, neutral transport is currently done through the same components as for plasma, i.e. evolve_density, evolve_momentum and evolve_pressure with the additional, optional neutral_parallel_diffusion component. In 2D, all of this functionality is implemented in one component called neutral_mixed which additionally has cross-field transport. This discrepancy is due to historical reasons and will be refactored.

1D: 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 = ... , braginskii_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 \frac{1}{p_n}{\partial_{||}p_n}\right) \\ \frac{\partial p_n}{\partial t} =& \ldots + \nabla\cdot\left(\mathbf{b}D_n p_n \frac{1}{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(m_nn_nv_{||n}\right) =& \ldots + \nabla\cdot\left(\mathbf{b}D_n m_n n_nv_{||n} \frac{1}{p_n} \partial_{||}p_n\right) + \nabla\cdot\left(\mathbf{b}\eta_n \partial_{||}v_{||n}\right) \end{aligned}\end{split}\]

The diffusion coefficient is in \(m^2/s\) and is calculated as

\[D_n = \left(\frac{B}{B_{pol}}\right)^2 \frac{eT_n}{m_{n} \nu}\]

where m_{n} is the neutral species mass in kg and \(\nu\) is the collision frequency (by default, this sums up all of the enabled neutral collisions from the collisions component as well as the charge exchange rate). The factor \(B / B_{pol}\) is the projection of the cross-field direction on the parallel transport, and is the dneut input setting. Currently, the recommended use case for this component is to represent the neutrals diffusing orthogonal to the target wall, and it is recommended to set dneut according to the field line pitch at the target.

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 outputVars(Options &state) override

Save variables to the output.

2D/3D: neutral_mixed#

The below describes the neutral_mixed component used for 2D and 3D simulations. Note that all dimensionalities are compatible with the neutral_boundary component which facilitates energy losses to the wall through neutral reflection.

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 [M.V. Umansky, J.N.M (2003)]. The Hermes-3 approach is more advanced in having a separate neutral pressure equation, similar to the new AFN (Advanced Fluid Neutral) model in SOLPS-ITER [N. Horsten, N.F. (2017)].

\[ \begin{align}\begin{aligned}\begin{aligned}\\\begin{split}\frac{\partial n_n}{\partial t} =& -\nabla\cdot\left(n_n\mathbf{b}v_{\parallel, n} + n_n\mathbf{v}_{\perp n}\right) \\ & + S \\ \frac{\partial}{\partial t}\left(m_nn_nv_{\parallel, n}\right) =& -m_n \nabla\cdot\left(n_n v_{\parallel, n} (\mathbf{b}v_{||n} + \mathbf{v}_{\perp n}\right)) \\ & - \nabla_{\parallel}p_n \\ & + \nabla \cdot (\eta_{n} (\nabla_{\perp} v_{\parallel n} + \mathbf{b} \nabla_{\parallel} v_{\parallel n} )) \\ & + F \\ \frac{\partial p_n}{\partial t} =& -\frac{5}{3} \nabla\cdot\left(p_n\mathbf{b}v_{\parallel, n} + p_n\mathbf{v}_{\perp n}\right) \\ & + \frac{2}{3}v_{\parallel, n} \nabla_{\parallel} p_n \\ & + \frac{2}{3} \nabla\cdot\left(\kappa_n \left(\nabla_\perp T_n + \mathbf{b} \nabla_{\parallel} T_n\right)\right) \\ & - \frac{2}{3} v_{\parallel,n} \nabla \cdot \left(\eta_{n} \left( \nabla_{\perp} v_{\parallel n} + \mathbf{b} \nabla_{\parallel} v_{\parallel n} \right) \right) \\ & + \frac{2}{3}E \\\end{split}\\\end{aligned}\end{aligned}\end{align} \]

Where for the density equation, the first row of terms contains the parallel and perpendicular advection and the second row the particle sources. In the parallel momentum equation, the first row of terms features parallel and perpendicular advection of parallel momentum. This is followed by the compression term and the perpendicular and parallel viscosity (diffusion of parallel momentum) as well as the momentum source term. In the pressure equation, the first row contains the parallel and perpendicular advection of pressure. This is followed by the compression term, the perpendicular and parallel conduction (diffusion of temperature) and perpendicular and parallel viscous heating, finally followed by the energy sources.

While parallel momentum is evolved and is exchanged with the plasma parallel momentum, the advection of momentum is neglected in the perpendicular direction, resulting in the pressure diffusion model, where the pressure gradient is balanced by frictional forces. This is similar to Fickian diffusion with the pressure gradient replacing the density gradient as the flow driver, in an approach similar to that taken in nuclear fission neutronic transport modelling and several other edge codes.

The perpendicular velocity is calculated as:

\[\begin{aligned} v_{\perp} =& -D_n \frac{1}{P_n} \nabla_{\perp} p_n \end{aligned}\]

Where in the code, \(\frac{1}{P_n} \nabla_{\perp}P_n\) is represented as \(ln(P_n)\), which helps preserve pressure positivity.

The diffusion coefficients are defined as:

\[\begin{split}\begin{aligned} D_n =& \frac{v_{th,n}^{2}}{\nu_{n, tot}} = \frac{T_n}{m_n \nu_{n, tot}} \\ \kappa_{n} =& \frac{5}{2} D_n N_n \\ \eta_{n} =& \frac{2}{5} m_n \kappa_{n} \\ \end{aligned}\end{split}\]

Where \(v_{th,n}= \sqrt{\frac{T_n}{m_n}}\) is the thermal velocity of neutrals and \(\nu_{n, tot}\) is the total neutral collisionality. When the AFN diffusion collision mode is selected using the diffusion_collisions_mode setting, this collisionality is the sum of charge exchange, ionisation, neutral-neutral collisions and the pseudo-collisionality Rnn, which represents a mean-free path limit. When the multispecies mode is selected, all available collision frequencies are enabled instead of ionisation. AFN is recommended in all cases, with the multispecies mode representing a legacy approach. The Rnn pseudo-collisionality is based on the neutral_lmax parameter, currently hardcoded to 0.1m, which acts as an effective maximum neutral mean free path. It represents the distance that neutrals can travel before hitting a solid surface.

In an additional effort to limit the diffusivitiy to more physical values, a flux limiter has been implemented which clamps \(D_n\) to \(D_{n,max}\) defined as:

\[\begin{aligned} D_{n,max} =& f_l \frac{v_{th,n}}{abs(\nabla ln(P_n) + 1/l_{max})} \end{aligned}\]

This formulation is equivalent to defining a \(D_n\) with a free streaming velocity while accounting for the pseudo collisionality due to the maximum vessel mean free path \(l_{max}\). The flux limiter \(f_l\) is set to 1.0 by default.

struct NeutralMixed : public Component

Evolve density, parallel momentum and pressure for a neutral gas species with cross-field diffusion

Public Functions

NeutralMixed(const std::string &name, Options &options, Solver *solver)
Parameters:

  • name – The name of the species e.g. “h”

  • options – Top-level options. Settings will be taken from options[name]

  • solver – Time-integration solver to be used

virtual void finally(const Options &state) override

Use the final simulation state to update internal state (e.g. time derivatives)

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.

2D/3D: neutral_full_velocity#

This model evolves the equations for a neutral fluid, assuming axisymmetry (constant in \(Z\)), for the density \(n_n\), velocity \(\mathbf{v}_n\) and pressure \(p_n\).

\[\begin{split}\begin{aligned} \frac{\partial n_n}{\partial t} =& -\nabla\cdot\left(n_n\mathbf{v}_n\right) \nonumber \\ \frac{\partial \mathbf{v}_n}{\partial t} =& - \mathbf{v}_n\cdot\nabla\mathbf{v}_n -\frac{1}{n_n m_n}\nabla p_n + \nabla\cdot\left(\nu \nabla\mathbf{v}\right)\\ \frac{\partial p_n}{\partial t} =& -\gamma \nabla\cdot\left(p_n\mathbf{v}_n\right) + \left(\gamma - 1\right)\mathbf{v}_n\cdot\nabla p_n + \nabla\cdot\left(n_n \chi_n \nabla T_n\right) \nonumber \end{aligned}\end{split}\]

where the adiabatic index \(\gamma\) and dissipation parameters \(\nu\) (kinematic viscosity) and \(\chi\) (thermal conduction) are constants set in the options:

[d]
type = neutral_full_velocity

adiabatic_index = 5./3 # Ratio of specific heats
viscosity = 1.0   # Kinematic viscosity [m^2/s]
conduction = 1.0  # Heat conduction [m^2/s]

The contravariant components of \(\mathbf{v}_n\) are evolved in the same \(\left(x,y,z\right)\) field-aligned coordinate system as the plasma. To evaluate the nonlinear advection term, whilst avoiding the use of noisy Christoffel symbols coming from derivatives of basis vectors, these components are transformed into \(\left(R,Z,\phi\right)\) cylindrical coordinates, advected, then transformed back. This is done using matrices which are calculated in the initialisation stage by finite differences of the input mesh:

\[\begin{split}\begin{aligned} \left(\begin{array}{c} \nabla R \\ \nabla Z\end{array}\right) =& \left(\begin{array}{cc} \frac{\partial R}{\partial x} & \frac{\partial R}{\partial y} \\ \frac{\partial Z}{\partial x} & \frac{\partial Z}{\partial y}\end{array}\right)\left(\begin{array}{c} \nabla x \\ \nabla y\end{array}\right) \\ =& \left(\begin{array}{cc} \texttt{Urx} & \texttt{Ury} \\ \texttt{Uzx} & \texttt{Uzy} \end{array}\right)\left(\begin{array}{c} \nabla x \\ \nabla y\end{array}\right) \end{aligned}\end{split}\]

These components are calculated by finite differences of the Rxy and Zxy arrays in the input, then adjusted to match the given values of hthe and Bpxy:

\[\sqrt{\left(\frac{\partial R}{\partial y}\right)^2 + \left(\frac{\partial R}{\partial y}\right)^2} = h_\theta\]
\[\sqrt{\left(\frac{\partial R}{\partial x}\right)^2 + \left(\frac{\partial R}{\partial x}\right)^2} = 1 / \left(R B_\theta\right)\]

(Note that this second equality only works if \(x\) and \(y\) are orthogonal).

This matrix is then inverted, to give:

\[\begin{split}\begin{aligned} \left(\begin{array}{c} \nabla x \\ \nabla y\end{array}\right) =& \left(\begin{array}{cc} \texttt{Txr} & \texttt{Tyr} \\ \texttt{Txz} & \texttt{Tyz} \end{array}\right)\left(\begin{array}{c} \nabla R \\ \nabla Z\end{array}\right) \end{aligned}\end{split}\]

The components of \(\mathbf{v}_n\) are evolved in contravariant form:

\[\mathbf{v}_n = v^x \mathbf{e}_x + v^y \mathbf{e}_y + v^z \mathbf{e}_z\]

These components are stored in the output. In the RHS function the velocity is converted to covariant form:

\[\mathbf{v}_n = v_x \nabla x + v_y \nabla y + v_z \nabla z\]

which is then transformed to \(v_r\), \(v_Z\) and \(v_\phi\):

\[\begin{split}\begin{aligned} v_r =& \mathbf{v}_n \cdot \nabla R = \frac{\partial x}{\partial R} v_x + \frac{\partial y}{\partial R} v_y \\ v_Z =& \mathbf{v}_n \cdot \nabla Z = \frac{\partial x}{\partial Z} v_x + \frac{\partial y}{\partial Z} v_y \\ v_\phi =& \mathbf{v}_n \cdot \hat{\phi} = v_z / \left(\sigma_{Bpol} R\right) \end{aligned}\end{split}\]

which are implemented as

Field2D vr = Txr * Vn2D.x + Tyr * Vn2D.y;
Field2D vz = Txz * Vn2D.x + Tyz * Vn2D.y;
Field2D vphi = Vn2D.z / (sigma_Bp * Rxy);

These components are then advected as scalars for the \(\mathbf{v}_n\cdot\nabla\mathbf{v}_n\) term, and are diffused for the \(\nabla\cdot\left(\mu \nabla\mathbf{v}\right)\) kinematic viscosity.

The advection of momentum \(\mathbf{v}\cdot\nabla\mathbf{v}\) is controlled by these settings:

  1. momentum_advection is false by default, disabling this nonlinear advection term. This keeps the inertia in the time derivative, but neglects the neutral dynamic pressure in the momentum balance.

  2. toroidal_flow is true by default, which includes the toroidal (\(z\)) component of the neutral flow. Importantly, this allows the parallel and poloidal flows to evolve independently: The parallel flow can follow the plasma towards the target, while the poloidal flow can be away from the target.

  3. curved_torus is true by default, and is only active when both momentum_advection and toroidal_flow are enabled. Neutrals travel in straight lines in real space, so toroidal flow is converted to radial flow. This appears in the \(v_r\) and \(v_\phi\) equations due to a combination of the radial centrifugal force and conservation of toroidal angular momentum.

Flow perpendicular to the magnetic field is damped by collisions e.g. CX reactions with the plasma. The steady-state flow is therefore a balance between the pressure gradient (including dynamic pressure if momentum_advection is enabled), and this friction. The neutral velocity perpendicular to the magnetic field is:

\[\begin{split}\begin{aligned} \mathbf{v}_{n\perp} =& \mathbf{v}_{n} - \mathbf{b}\mathbf{b}\cdot\mathbf{v}_{n} \\ =& \mathbf{v}_{n} - \mathbf{e}_y\frac{v_{ny}}{g_{yy}} \\ =& \mathbf{v}_{n} - \left(\nabla y + \frac{g_{yz}}{g_{yy}}\nabla z\right)v_{ny} \\ =& v_{nx}\nabla x + \left(v_{nz} - \frac{g_{yz}}{g_{yy}}v_{ny}\right)\nabla z \end{aligned}\end{split}\]

At boundaries neutral thermal energy is lost at a rate controlled by the option

neutral_gamma = 5./4

This sets the flux of power to the wall to:

\[q = \gamma n_n T_n c_s\]

Currently this is only done at target boundaries, not radial boundaries.

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.

Two forms are available, which are implemented differently for density, momentum, and pressure equations. In the density equation, form 0 uses the diamagnetic velocity perpendicular to b and the gradient of P; at the boundaries this velocity is perpendicular to the boundary. Form 1 uses the magnetic gyro-center drifts, which are mostly vertical; at the boundaries this form produces a flow through the boundary. Forms 0 and 1 are analytically equivalent and should give the same result away from boundaries, but form 0 doesn’t produce flows through boundaries. This is an approach that UEDGE uses to avoid unphysical boundary flows.

However, Form 1 is nice because the flow velocity depends on the temperature, not the pressure gradient. This usually makes it better behaved numerically. To make the most of both, the diamagnetic_drift component allows the forms to be mixed using the diamag_form setting. For example, the tcv-x21 example blends it such that form 0 is at the boundary:

[diamagnetic_drift]
diamag_form = x * (1 - x)  # 0 = gradient; 1 = divergence

A table of the two forms used in Hermes-3, and the corresponding terms in Simakov & Catto is shown below, where \(\mathbf{C}=\nabla\times\left(\frac{\mathbf{b}}{B}\right)\) is the curvature vector. Instead of the diamagnetic velocity, the whole terms associated are shown. The difference among the forms is the divergence of a curl, which vanishes. The diamagnetic velocity \(\mathbf{v}_{dia}\) is defined above. Notice that Simakov & Catto used Gaussian units, but Hermes-3 uses SI units.

Form 0

Form 1

Simakov & Catto

Density

\(\mathbf{C} \cdot \nabla\left(\dfrac{p}{q}\right)\)

\(\nabla \cdot (n \mathbf{v}_{dia})\)

Eq. (51): \(\dfrac{c}{q}\left(\nabla \times \dfrac{\mathbf{b}}{B}\right) \cdot \nabla p\)

Momentum

\(\mathbf{C} \cdot \nabla\left(\dfrac{mnv_\parallel T}{q}\right)\)

\(\nabla\cdot (mnv_\parallel \mathbf{v}_{dia})\)

Eq. (64): \(\nabla\cdot \left(\dfrac{1}{\Omega} \mathbf{b} \times \nabla(p v_\parallel)\right)\)

Pressure

\(\dfrac{5}{2}\mathbf{C} \cdot \nabla\left(\dfrac{pT}{q}\right)\)

\(\dfrac{5}{2}\nabla\cdot (p \mathbf{v}_{dia})\)

Eq. (56): \(\nabla\cdot\left(\dfrac{5}{2 m \Omega} \mathbf{b} \times \nabla(p T)\right)\)

* Eq.(64) in Simakov & Catto is derived for ion parallel momentum, but it is also applicable to electrons since it comes from the gyro-viscosity and the mass factors of \(m_i\) or \(m_e\) cancel out.

struct DiamagneticDrift : public Component

Calculate diamagnetic flows.

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 outputVars(Options &state) override

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

Field3D calcDivJ(GuardedOptions &state)

Calculate divergence of all currents except polarisation.

void diamagneticCompression(GuardedOptions &state, Field3D DivJ)

Calculate energy exchange term nonlinear in pressure due to compression of polarisation drift (3 / 2) ddt(Pi) += (Pi * m_i / n0 / Z) * DivJ

Adds energy_source for all species that have charge, mass and pressure Throws a BoutException if boussinesq=true

Field3D calcMassDensity(GuardedOptions &state)

Solve for time derivative of potential Using Div(mass_density / B^2 Grad_perp(dphi/dt)) = DivJ

Field3D calcPolFlowPotential(Field3D mass_density, Field3D DivJ)

Calculate poloidal drift potential-flow approximation.

void polarisationAdvection(GuardedOptions &state, Field3D phi_pol)

Polarisation drift approximated by a potential flow

v_p = - (m_i / (Z_i * B^2)) * Grad(phi_pol)

Sets density_source, energy_source and momentum_source for all species with mass and charge.

Stellarator cross-field transport: 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 outputVars(Options &state) override

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

Tokamak cross-field transport: 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 outputVars(Options &state) override

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

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.

Kinematic viscosity can be included by setting e.g. viscosity = 0.1 in SI units (m^2/s). This adds a diffusion of vorticity and corresponding ion heating. The viscous friction force in this simplified operator is

\[\mathbf{F}_\nu = - \nu B \mathbf{b} \times \nabla \Omega\]

This gives rise to a drift and current with divergence:

\[\begin{split}\begin{aligned}\nabla\cdot\mathbf{J_{\nu}} =& \nabla\cdot\left[\frac{\mathbf{b}\times\mathbf{F}_\nu}{B}\right] \\ =& \nabla\cdot\left[\nu \nabla_\perp \Omega\right]\end{aligned}\end{split}\]

Viscous heating is calculated using the work done by a fluid velocity consistent with the Boussinesq approximation:

\[\mathbf{u} = \frac{\mathbf{b}\times\nabla\Phi}{B}\]

where the generalized potential is

\[\Phi = \phi + \hat{p} / \overline{n}\]

The work done is

\[\mathbf{F}_\nu\cdot\mathbf{u} = -\nu \nabla_\perp\Omega \cdot \nabla_\perp\Phi\]

This heating is distributed between charged species in proportion to their local mass density. The properties of the work done can be analysed by writing in terms of a vector \(\mathbf{g}\):

\[\mathbf{g} = \frac{\overline{A}\overline{n}}{B^2}\nabla_\perp\Phi\]

to write:

\[\begin{split}\begin{aligned}\mathbf{F}_\nu\cdot\mathbf{u} =& -\nu\frac{B^2}{\overline{A}\overline{n}} \nabla_\perp\left(\nabla\cdot\mathbf{g}\right)\cdot\mathbf{g} \\ =& \nu\frac{B^2}{\overline{A}\overline{n}} \left[\left(\nabla_\perp\mathbf{g} : \nabla_\perp\mathbf{g}\right) - \nabla\cdot\left(\mathbf{g}\cdot\nabla\mathbf{g}\right)\right]\end{aligned}\end{split}\]

The last term is not in general positive definite, so this simple form of viscosity could in some cases lead to cooling.

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_core_averagey: bool, default false Average phi core boundary in Y? (if phi_boundary_relax)

    • 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 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.

Field3D calculatePihat(GuardedOptions allspecies)

Diamagnetic term in vorticity. Note this is weighted by the mass This includes all species, including electrons

Field3D calculateDivJdia(Field3D phi, GuardedOptions allspecies)

Calculates Div(Jdia) and sets energy_source for all charged species with pressure.

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>

    • 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)?

    • 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_core_averagey: bool, default false Average phi core boundary in Y? (if phi_boundary_relax)

    • 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?

    • viscosity_perp: Field2D, default 0.0 Kinematic viscosity in perpendicular diraction [m^2/s]

    • viscosity_par: Field2D, default 0.0 Kinematic viscosity in parallel diraction [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 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#

Notes: When using this module,

  1. Set sound_speed:alfven_wave=true so that the shear Alfven wave speed is included in the calculation of the fastest parallel wave speed for numerical dissipation.

  2. For tokamak simulations use Neumann boundary condition on the core and Dirichlet on SOL and PF boundaries by setting electromagnetic:apar_core_neumann=true (this is the default).

  3. Set the potential core boundary to be constant in Y by setting vorticity:phi_core_averagey = true

  4. Magnetic flutter terms must be enabled to be active (electromagnetic:magnetic_flutter=true). They use an Apar_flutter field, not the Apar field that is calculated from the induction terms.

  5. If using vorticity:phi_boundary_relax to evolve the radial boundary of the electrostatic potential, the timescale phi_boundary_timescale should be set much longer than the Alfven wave period or unphysical instabilities may grow from the boundaries.

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 requires an additional term in the momentum equation:

\[\frac{\partial p_s}{\partial t} = \cdots + Z_s e A_{||} \frac{\partial n_s}{\partial t}\]

This is implemented so that the density time-derivative is calculated using the lowest order terms (parallel flow, ExB drift and a low density numerical diffusion term).

The above equations are 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\]

The toroidal variation of density \(n_s\) must be kept in this equation. By default the iterative “Naulin” solver is used to do this: A fast FFT-based method is used in a fixed point iteration, correcting for the density variation.

Magnetic flutter terms are disabled by default, and can be enabled by setting

[electromagnetic]
magnetic_flutter = true

This writes an Apar_flutter field to the state, which then enables perturbed parallel derivative terms in the evolve_density, evolve_pressure, evolve_energy and evolve_momentum components. Parallel flow terms are modified, and parallel heat conduction.

\[\begin{split}\begin{aligned}\mathbf{b}\cdot\nabla f =& \mathbf{b}_0\cdot\nabla f + \delta\mathbf{b}\cdot\nabla f \\ =& \mathbf{b}_0\cdot\nabla f + \frac{1}{B}\nabla\times\left(\mathbf{b}A_{||}\right)\cdot\nabla f \\ \simeq& \mathbf{b}_0\cdot\nabla f + \frac{1}{B_0}\left[A_{||}\nabla\times\mathbf{b} + \left(\nabla A_{||}\right)\times\mathbf{b}_0\right]\cdot\nabla f \\ \simeq& \mathbf{b}_0\cdot\nabla f + \frac{1}{B_0}\mathbf{b}_0\times \nabla A_{||} \cdot \nabla f\end{aligned}\end{split}\]
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 restartVars(Options &state) override

Add extra fields to restart files.

virtual void outputVars(Options &state) override

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

Contents