Computational Plasma Models

Global Models

• Volume Averaged and Simplified Kinetics

• Spatially independent simulations provide:

  • Volume-averaged plasma parameters.

  • Energy-dependent coefficients.

Fuild Models (Zapdos)

• Local Field Approximation (LFA).

• Local Mean Energy Approximation (LMEA).

• Monte Carlo Fluids (MCF).

• Typically utilizes the first two or three moments of the Boltzmann equation

• Particle species are represented as continuous fluids.

  • Various assumptions can be used for energy-dependent phenomena.

Higher Order and Hybrid Models

• Higher order fluid models

• Hybrid particle-fluid models

• Higher order models utilize higher order moments of the Boltzmann equation.

• Hybrid models represent a fraction of a species, a single species, or multiple species in a fully-kinetic manner.

Kinetic Particle Methods

• Direct Kinetic Methods

• Molecular Dynamics (MD)

• Particle In Cell (PIC) with Monte Carlo Collisions (MCC)

• Individual particles or groups of particles are tracked.

  • Direct Kinetic Models explicity solve the Boltzmann equation.

  • In MD models, inter-molecular forces between particles are calculated and used to update paricle data.

  • In PIC-MCC models, individual particles are traced and are subject to macroscopic fields and particle-particle Monte Carlo collisions.

MOOSE Introduction

History and Purpose

• Development started in 2008

• Open-sourced in 2014

• Designed to solve computational engineering problems and reduce the expense and time required to develop new applications by:

  • being easily extended and maintained

  • working efficiently on a few and many processors

  • providing an object-oriented, extensible system for creating all aspects of a simulation tool

MOOSE By The Numbers

• 208 contributors

• 46,000 commits

• 5000 unique visitors per month

• ~6 new Discussion participants per week

• 1500 citations for the MOOSE papers

  • Most cited paper in Elsevier Software-X

  • More than 500 publications using MOOSE

• 30M tests per week

General Capabilities

• Continuous and Discontinuous Galerkin FEM

• Finite Volume

• Supports fully coupled or segregated systems, fully implicit and explicit time integration

• Automatic differentiation (AD)

• Unstructured mesh with FEM shapes

• Higher order geometry

• Mesh adaptivity (refinement and coarsening)

• Massively parallel (MPI and threads)

• User code agnostic of dimension, parallelism, shape functions, etc.

• Operating Systems:

  • macOS

  • Linux

  • Windows (WSL)

Object-oriented, pluggable system

Example Code

Software Quality

• MOOSE follows an Nuclear Quality Assurance Level 1 (NQA-1) development process

• All commits undergo review using GitHub Pull Requests and must pass a set of application regression tests before they are available to our users

• MOOSE includes a test suite and documentation system to allow for agile development while maintaining a NQA-1 process

• Utilizes the Continuous Integration Environment for Verification, Enhancement, and Testing (CIVET)

• External contributions are guided through the process by the team, and are very welcome!

Development Process

Community

https://github.com/idaholab/moose/discussions

License

• LGPL 2.1

• Does not limit what you can do with your application

  • Can license/sell your application as closed source

• Modifications to the library itself (or the modules) are open source

• New contributions are automatically LGPL 2.1

Zapdos Introduction

Zapdos Overview

• Zapdos is an open source multi-fluid plasma code, that:

  • Utilizes the finite element method, and

  • Solves for species densities in log-molar form.

• As a MOOSE-based application, Zapdos can take advantage of:

  • Scalable parallel computing, and

  • Flexible multiphysics coupling.

Current Research Efforts

• Validation and Verification

  • Modeling physical systems and comparing to experimental results.

  • Utilizing the Method of Manufactured Solutions (MMS) to ensure correct implementation.

Experimental Comparison with results from Greenberg and Hebner (1993).

1D MMS Convergence study results for a single ionic species, electrons and mean electron energy.

Code Development

• Improving physics

• Electromagnetic coupling

• Improving the user interface

• Finite volume implementation

Some selected results

Electron Oscillation

2D Electron Oscilation in the GEC reference cell

Experimental Validation: Metastable Profiles

Argon Ground To Plate McMillin and Zachariah (1995)

75 [V] peak to peak at 100 [mTorr]

75 [V] peak to peak at 100 [Torr]

MOOSE Electromagnetics

Cost Jet Gas Mixing

Zapdos Multi-fluid Equations

Logarithmic and Logarithmic-Molar Form

In Zapdos, the logarithmic form or the logarithmic-molar form is used to calculate the species density. This is done for numerical stability and to prevent the calculation of nonphysical, negative densities by the numerical solver. The logarithmic form is given by

and the logarithmic-molar form is given by

Where is the traditional species density, is the modified species density, and is Avogadro's number

Governing Equations

Drift-Diffusion Approximation

where and denote some particle species, denotes the number of bodies in the source term, and denotes sources of species from chemical reactions.

where and are the mobility and diffusivity of species respectively, and is the electric potential.

where denotes charged particles and denotes neutral particles.

Source Terms (from Crane)

Where is the stoichiometric coefficient of the reaction, is the reaction rate, and is the density of species .

Electromagnetic Field Calculations

When evaluating an electrostatic system, Poisson's equation is used to calculate the electrostatic potential,

where is the elementary charge, is the permitivity of free space, and is the charge density in the system.

Additionally, Zapdos is capable of calculating an effective electric field using the form

In cases where electromagnetic systems are considered, MOOSE's Electromagnetics Module is used to perform field calculations.

Energy Conservation

Boundary Conditions (BCs)

Electrostatic Potential Boundary Conditions

Flux Boundary Conditions: Type 1

Flux Boundary Conditions: Type 2

Flux Boundary Conditions: Type 3

Brief Visualization Overivew

Visualizing 1D results

The results shown here are from Tutorial 1

Visualizing results from a 0D model

The results shown here are from Tutorial 2

Tutorial 1 - Diffusion + Constant Source

Problem Statement

Consider a one dimensional single species, , plasma on the domain . Additionally, consider that this is a purely diffusive system with a constant source term due to a constant background gas density . This system takes the following form

Where is a constant diffusion coefficient, is the species density, and is a constant first order reaction rate. Additionally, a zero-density boundary condition will be imposed

Explore what happens as the diffusion coefficient is varied.

[Materials]
  [gas_species_0]
    type = ADHeavySpecies
    heavy_species_name = Ar+
    heavy_species_mass = 6.64e-26
    heavy_species_charge = 1.0
    diffusivity = 1
    #diffusivity = 2
  []
[]

Explore what happens as the reaction rate is varied.

[Materials]
  [FirstOrder_Reaction]
    type = GenericRateConstant
    reaction = FirstOrder
    reaction_rate_value = 9.8696
    #reaction_rate_value = 4.9348
  []
[]

Analytic Solution

Diffusive Species Solution.

Tutorial 2 - Reaction Network

Problem Statement

Consider a system of 3 species, .

• decays at a rate of into species

• decays at a rate of into species

• Species is a stable species

The initial condition for this system is given by

System of Equations

How would the reaction rates need to change to decrease the decay of and increase the growth of ?

[Reactions]
  [Gas]
    #Name of each variable on the reactant side
    species = 'nA nB nC'
    #Define type of coefficient (rate or townsend)
    reaction_coefficient_format = 'rate'
    #Define if using log form
    use_log = false
    #Define if using automatic differentiation
    use_ad = true
    #Define which material block the reactions take place.
    # For undefine blocks, naming starts at 0
    block = 0
    #Define reactions and coefficients
    reactions = 'nA -> nB  : 1
                 nB -> nC  : 5'
  []
[]

Analytic Solution

Species evolution over time.

Tutorial 3 - Electric Field Driven Ion Loss

Problem Statement

Consider a single species plasma, an ionic Argon plasma, in a one-dimensional system on the domain . In this system

• There is an initial uniform distribution of singly charged Argon ions;

• The ion-only plasma is subject to the electric field generated by its own charge density; and

• The walls of this system destroy (or remove) ions.

This is the same problem as described on page 26-27 of Lieberman and Lichtenberg (1994).

Sytem of Equations

What would happen to the rate of the potential decay as pressure changes?

[Materials]
  [GasBasics]
    type = GasElectronMoments
    interp_trans_coeffs = false
    interp_elastic_coeff = false
    ramp_trans_coeffs = false
    user_p_gas = 13.3322
    #user_p_gas = 1.33322
    potential = potential
    property_tables_file = rate_coefficients/electron_moments.txt
  []
[]

Solution Demonstration

Tutorial 3 - Kinetic Implementation using FENIX

FENIX Intro

As commercial interest in fusion energy has increased, the need for high-fidelity, multiphysics simulations of fusion devices has also increased. To address this, Idaho National Laboratory, North Carolina State University, the University of Illinois at Urbana-Champaign, and the United Kingdom Atomic Energy Authority are developing an open source application, the Fusion ENergy Integrated multiphys-X (FENIX) framework, for modeling plasma facing components. FENIX is built on the MOOSE framework and couples MOOSE ecosystem capabilities as well as external codes:

• TMAP8 (Tritium migration)

• Cardinal

  • Monte-Carlo Neutronics via OpenMC

  • Computational Fluid Dynamics via NekRS

• Electromagnetics Module

• Ray Tracing Module (Kinetic Plasma Foundation)

• Heat Transfer

• Solid Mechanics

• Thermal Hydraulics

PIC in FENIX

In Lieberman and Lichtenberg (1994), this tutorial problem is solved with a PIC code, a core capability of FENIX still in-development. As a very basic benchmarking problem, the PIC example presented in Lieberman and Lichtenberg (1994) has been replicated in FENIX.

• 100 Macroparticles are placed on a 100 element mesh

• Each particle is initially placed exactly in the center of each cell

• Leap Frog particle stepping is used

• =

• Total simulation time of

• The potential is represented with first order nodal basis functions

Results Comparison

All figures show FENIX results superimposed on figures 2.2.a, 2.2.b, and 2.2.c from Lieberman and Lichtenberg (1994). Results from Lieberman and Lichtenberg (1994) are in black and grey while FENIX results are in color.

Particle Population Evolution

Electrostatic Potential Evolution

Particle Phase Space Evolution

Current PIC Status

Verified Capabilities

• Leap frog and Boris particle stepping

• Electrostatic capabilities using all first order finite elements in libMesh.

• Uniform random particle initialization for all first order finite elements in libMesh.

  • Bounding box particle initialization.

Under Development

• Particle based current density source calculations.

• Particle collision capabilities.

Future PIC Work

• Continued verification efforts by replicating Kinetic instabilities

  • Landau Damping

  • Two Stream

  • Dory-Ghast-Haris Instability

• Replicating more complex analytic plasma solutions

• Benchmarking against other codes

• Validation by replicating experimental conditions and data

Tutorial 4 - Pressure versus Electron Temperature Parametric Study

Problem Statement

Consider an Argon capacitively coupled plasma (CCP). Utilize the MOOSE Postprocessor system to calculate the volume-averaged electron temperature.

How is the average electron temperature affected as the pressure is varied?

Background Gas Density

[AuxKernels]
  [Ar_val]
    type = FunctionAux
    variable = Ar
    function = 'log(3.22e22/6.022e23)'
    execute_on = INITIAL
  []
[]

Pressure in the electron material block, and ion mobility and diffusivity coefficients.

  [GasBasics]
    type = GasElectronMoments
    #True means variable electron coeff, defined by user
    interp_trans_coeffs = true
    #Leave as false (CRANE accounts of elastic coeff.)
    interp_elastic_coeff = false
    #Leave as false, unless computational error is due to rapid coeff. changes
    ramp_trans_coeffs = false
    #Name for electrons (usually 'em')
    em = em
    #Name for potential (usually 'potential')
    potential = potential
    #Name for the electron mean energy density (usually 'mean_en')
    mean_en = mean_en
    #User difine pressure in pa
    user_p_gas = 133.322
    #True if pressure dependent coeff.
    pressure_dependent_electron_coeff = true
    #Name of text file with electron properties
    property_tables_file = rate_coefficients/electron_moments.txt
  []
  #The material properties of the ion
  [gas_species_0]
    type = ADHeavySpecies
    heavy_species_name = Ar+
    heavy_species_mass = 6.64e-26
    heavy_species_charge = 1.0
    mobility = 0.144409938
    diffusivity = 6.428571e-3
  []

Expected Output

Tutorial 5 - Plasma-Water Interface

Problem Statement

Consider an Argon DC discharge over water. Within Zapdos, users can model material surface interfaces of plasma processes (in this case, electrons diffusing into a water’s surface). This is done using what is called the Interface System within Zapdos/MOOSE.

  [em_physical_right]
    type = HagelaarElectronBC
    variable = em
    boundary = 'master0_interface'
    potential = potential
    electron_energy = mean_en
    r = 0.00
    position_units = ${dom0Scale}
  []
  [Arp_physical_right_diffusion]
    type = HagelaarIonDiffusionBC
    variable = Arp
    boundary = 'master0_interface'
    r = 0
    position_units = ${dom0Scale}
  []

Expected Output