Examples
1D flux-tube
These simulations follow the dynamics of one or more species along the magnetic field. By putting a source at one end of the domain, and a sheath at the other, this can be a useful model of plasma dynamics in the Scrape-Off Layer (SOL) of a tokamak or other magnetised plasma.
1D periodic domain, Te and Ti
A fluid is evolved in 1D, imposing quasineutrality and zero net current. Both electron and ion pressures are evolved, but there is no exchange of energy between them, or heat conduction.
To run this example:
./hermes-3 -d examples/1D-te-ti
Which takes a few seconds to run on a single core. Then in the
examples/1D-te-ti
directory run the analysis script
python3 makeplot.py
That should generate png files and an animated gif if ImageMagick is
installed (the convert
program). If an error like
ModuleNotFoundError: No module named 'boutdata'
occurs, then
install the boutdata
package with python3 -m pip install
boutdata
.
The model components are ions (i) and electrons (e), and a component which uses the force on the electrons to calculate the parallel electric field, which transfers the force to the ions.
[hermes]
components = i, e, electron_force_balance
The ion density, pressure and momentum equations are evolved:
[i] # Ions
type = evolve_density, evolve_pressure, evolve_momentum
which solves the equations
The electron density is set to the ion density by quasineutrality, the parallel velocity is set by a zero current condition, and only the electron pressure is evolved.
[e] # Electrons
type = quasineutral, zero_current, evolve_pressure
which adds the equations:
The zero_current component sets:
1D Scrape-off Layer (SOL)
This simulates a similar setup to the SD1D code: A 1D domain, with a source of heat and particles on one side, and a sheath boundary on the other. Ions recycle into neutrals, which charge exchange and are ionised. A difference is that separate ion and electron temperatures are evolved here.
Due to the short length-scales near the sheath, the grid is packed close to the target, by setting the grid spacing to be a linear function of index:
[mesh]
dy = (length / ny) * (1 + (1-dymin)*(1-y/pi))
where dymin
is 0.1 here, and sets the smallest grid spacing (at the
target) as a fraction of the average grid spacing.
The components are ion species d+
, atoms d
, electrons e
:
[hermes]
components = (d+, d, e,
zero_current, sheath_boundary, collisions, recycling, reactions,
neutral_parallel_diffusion)
The electron velocity is set to the ion by specifying zero_current; A sheath boundary is included; Collisions are needed to be able to calculate heat conduction, as well as neutral diffusion rates; Recycling at the targets provides a source of atoms; neutral_parallel_diffusion simulates cross-field diffusion in a 1D system.
The sheath boundary is only imposed on the upper Y boundary:
[sheath_boundary]
lower_y = false
upper_y = true
The reactions component is a group, which lists the reactions included:
[reactions]
type = (
d + e -> d+ + 2e, # Deuterium ionisation
d + d+ -> d+ + d, # Charge exchange
)
To run this example:
nice -n 10 ./hermes-3 -d examples/1D-recycling
This should take 5-10 minutes to run. There is a makeplots.py
script in the
examples/1D-recycling
directory which will generate plots and a gif animation
(if ImageMagick is installed).
2D drift-plane
Simulations where the dynamics along the magnetic field is not included, or only included in a parameterised way as sources or sinks. These are useful for the study of the basic physics of plasma “blobs” / filaments, and tokamak edge turbulence.
Blob2d
A seeded plasma filament in 2D. This version is isothermal and cold ion, so only the electron density and vorticity are evolved. A sheath-connected closure is used for the parallel current.
The model components are
[hermes]
components = e, vorticity, sheath_closure
The electron component consists of two types:
[e] # Electrons
type = evolve_density, isothermal
The evolve_density component type evolves the electron density Ne
. This component
has several options, which are set in the same section e.g.
poloidal_flows = false # Y flows due to ExB
and so solves the equation:
The isothermal component type sets the temperature to be a constant, and using
the density then sets the pressure. The constant temperature is also
set in this [e]
section:
temperature = 5 # Temperature in eV
so that the equation solved is
where \(T_e\) is the fixed electron temperature (5eV).
The vorticity component uses the pressure to calculate the diamagnetic current,
so must come after the e
component. This component then calculates the potential.
Options to control the vorticity component are set in the [vorticity]
section.
The sheath_closure
component uses the potential, so must come after vorticity.
Options are also set as
[sheath_closure]
connection_length = 10 # meters
This adds the equation
where \(L_{||}\) is the connection length.
Blob2D-Te-Ti
A seeded plasma filament in 2D. This version evolves both electron and ion temperatures. A sheath-connected closure is used for the parallel current.
The model components are
[hermes]
components = e, h+, vorticity, sheath_closure
The electron component evolves density (saved as Ne
) and pressure
(Pe
), and from these the temperature is calculated.
[e]
type = evolve_density, evolve_pressure
The ion component sets the ion density from the electron density, by
using the quasineutrality of the plasma; the ion pressure (Ph+
) is evolved.
[h+]
type = quasineutral, evolve_pressure
The equations this solves are similar to the previous Blob2d case, except now there are pressure equations for both ions and electrons:
2D-drift-plane-turbulence-te-ti
A 2D turbulence simulation, similar to the Blob2D-Te-Ti case, but with extra source and sink terms, so that a statistical steady state of source-driven turbulence can be reached.
The model components are
[hermes]
components = e, h+, vorticity, sheath_closure
The electron component evolves density (saved as Ne
) and pressure
(Pe
), and from these the temperature is calculated.
[e]
type = evolve_density, evolve_pressure
The ion component sets the ion density from the electron density, by
using the quasineutrality of the plasma; the ion pressure (Ph+
) is evolved.
[h+]
type = quasineutral, evolve_pressure
The sheath closure now specifies that additional sink terms should be added
[sheath_closure]
connection_length = 50 # meters
potential_offset = 0.0 # Potential at which sheath current is zero
sinks = true
and radially localised sources are added in the [Ne]
, [Pe]
, and [Ph+]
sections.
The equations this solves are the same as the previous Blob2D-Te-Ti case, except wih extra source and sink terms. In SI units (except temperatures in eV) the equations are:
Where \(\overline{T}\) and \(\overline{n}\) are the reference temperature (units of eV) and density (in units of \(m^{-3}\)) used for normalisation. \(\overline{c_s} = \sqrt{e\overline{T} / m_p}\) is the reference sound speed, where \(m_p\) is the proton mass. The mean ion atomic mass \(\overline{A}\) is set to 1 here.
These reference values enter into the sheath current \(\mathbf{j}_{sh}\) because that is a simplified, linearised form of the full expression. Likewise the vorticity (\(\omega\)) equation used the Boussinesq approximation to simplify the polarisation current term, leading to constant reference values being used.
The sheath heat transmission coefficients default to \(\gamma_e = 6.5\) and \(\gamma_i = 2.0\) (\(\gamma_i\) as suggested in Stangeby’s textbook between equations (2.92) and (2.93)). Note the sinks in may not be correct or the best choices, especially for cases with multiple ion species; they were chosen as being simple to implement by John Omotani in May 2022.
2D axisymmetric tokamak
These are transport simulations, where the cross-field transport is given by diffusion, and fluid-like equations are used for the parallel dynamics (as in the 1D flux tube cases).
The input settings (in BOUT.inp) are set to read the grid from a file tokamak.nc
.
This is linked to a default file compass-36x48.grd.nc
, a COMPASS-like lower single
null tokamak equilibrium. Due to the way that BOUT++ uses communications between
processors to implement branch cuts, these simulations require a multiple of 6 processors.
You don’t usually need 6 physical cores to run these cases, if MPI over-subscription
is enabled.
heat-transport
In examples/tokamak/heat-transport
, this evolves only electron pressure with
a fixed density. It combines cross-field diffusion with parallel heat conduction
and a sheath boundary condition.
To run this simulation with the default inputs requires (at least) 6 processors because it is a single-null tokamak grid. From the build directory:
cd examples/tokamak
mpirun -np 6 ../../hermes-3 -d heat-transport
That will read the grid from tokamak.nc
, which by default links to
the compass-36x48.grd.nc
file.
The components of the model are given in heat-transport/BOUT.inp
:
[hermes]
components = e, h+, collisions, sheath_boundary_simple
We have two species, electrons and hydrogen ions, and add collisions between them and a simple sheath boundary condition.
The electrons have the following components to fix the density, evolve the pressure, and include anomalous cross-field diffusion:
[e]
type = fixed_density, evolve_pressure, anomalous_diffusion
The fixed_density
takes these options:
AA = 1/1836
charge = -1
density = 1e18 # Fixed density [m^-3]
so in this simulation the electron density is a uniform and constant value.
If desired, that density can be made a function of space (x
and y
coordinates).
The evolve_pressure
component has thermal conduction enabled, and outputs
extra diagnostics i.e. the temperature Te
:
thermal_conduction = true # Spitzer parallel heat conduction
diagnose = true # Output additional diagnostics
There are other options that can be set to modify the behavior,
such as setting kappa_limit_alpha
to a value between 0 and 1 to impose
a free-streaming heat flux limit.
Since we’re evolving the electron pressure we should set initial and
boundary conditions on Pe
:
[Pe]
function = 1
bndry_core = dirichlet(1.0) # Core boundary high pressure
bndry_all = neumann
That sets the pressure initially uniform, to a normalised value of 1,
and fixes the pressure at the core boundary. Other boundaries are set
to zero-gradient (neumann) so there is no cross-field diffusion of heat out of
the outer (SOL or PF) boundaries. Flow of heat through the sheath is
governed by the sheath_boundary_simple
top-level component.
The hydrogen ions need a density and temperature in order to calculate the collision frequencies. If the ion temperature is fixed to be the same as the electron temperature then there is no transfer of energy between ions and electrons:
[h+]
type = quasineutral, set_temperature
The quasineutral
component sets the ion density so that there is no net charge
in each cell. In this case that means the hydrogen ion density is set equal to
the electron density. To perform this calculation the component requires that the
ion atomic mass and charge are specified:
AA = 1
charge = 1
The set_temperature
component sets the ion temperature to the temperature of another
species. The name of that species is given by the temperature_from
option:
temperature_from = e # Set Th+ = Te
The collisions
component is described in the manual, and calculates both electron-electron
and electron-ion collisions. These can be disabled if desired, using individual options.
There are also ion-ion, electron-neutral, ion-neutral and neutral-neutral collisions that
are not used here.
The sheath_boundary_simple
component is a simplified Bohm-Chodura sheath boundary
condition, that allows the sheath heat transmission coefficient to be specified for
electrons and (where relevant) for ions.
The equations solved by this example are:
The calculation of the Coulomb logarithms follows the NRL formulary,
and the above expression is used for temperatures above 10eV. See
the collisions
manual section for the expressions used in other regimes.
recycling-dthene
The recycling-dthene
example includes cross-field diffusion,
parallel flow and heat conduction, collisions between species, sheath
boundary conditions and recycling. It simulates the density, parallel
flow and pressure of the electrons; ion species D+, T+, He+, Ne+; and
neutral species D, T, He, Ne.
The model components are a list of species, and then collective components which couple multiple species.
[hermes]
components = (d+, d, t+, t, he+, he, ne+, ne, e,
collisions, sheath_boundary, recycling, reactions)
Note that long lists like this can be split across multiple lines by using parentheses.
Each ion species has a set of components, to evolve the density, momentum and pressure. Anomalous diffusion adds diffusion of particles, momentum and energy. For example deuterium ions contain:
[d+]
type = evolve_density, evolve_momentum, evolve_pressure, anomalous_diffusion
AA = 2
charge = 1
Atomic reactions are specified as a list:
[reactions]
type = (
d + e -> d+ + 2e, # Deuterium ionisation
t + e -> t+ + 2e, # Tritium ionisation
he + e -> he+ + 2e, # Helium ionisation
he+ + e -> he, # Helium+ recombination
ne + e -> ne+ + 2e, # Neon ionisation
ne+ + e -> ne, # Neon+ recombination
)