Mutation++: MUlticomponent Thermodynamic And Transport properties for IONized gases in C++

The Mutation++ library provides accurate and efficient computation of physicochemical properties associated with partially ionized gases in various degrees of thermal nonequilibrium. With v1.0.0, users can compute thermodynamic and transport properties, multiphase linearly-constrained equilibria, chemical production rates, energy transfer rates, and gas-surface interactions. The framework is based on an object-oriented design in C++, allowing users to plug-and-play various models, algorithms, and data as necessary. Mutation++ is available open-source under the GNU Lesser General Public License v3.0.


Motivation and Significance
The evaluation of thermochemical nonequilibrium, partially ionized gas properties is essential for a wide range of applications, including hypersonic flows, solar physics and space weather, ion thrusters, medical plasmas, combustion processes, meteor phenomena, and biomass pyrolysis. For example, the prediction of hypersonic flow plays an important role in the development of thermal protection systems for atmospheric entry vehicles. Such flows span a broad range of temporal scales, from local thermodynamic equilibrium to thermochemical nonequilibrium. As such, myriads of physicochemical models, data, and algorithms are used in today's hypersonic Computational Fluid Dynamics (CFD) codes and represent a significant body of work in the scientific literature. The thermochemical models employed in these codes directly affect the evaluation of gas properties necessary to close the conservation laws governing the fluid. These include mixture thermodynamic and transport properties, species chemical production rates, and energy transfer rates. Each of these properties further depends on the selection of a variety of specialized algorithms and data, such as species partition functions, transport collision integrals, and reaction rate coefficients.
The implementation, testing, and maintenance of the models, algorithms, and data required to simulate thermal nonequilibrium flows represent a significant cost, in terms of human resources and time necessary to develop a simulation tool. As new models, algorithms, or data become available, additional effort is required to update existing codes, especially when models are "hardcoded." A number of commercial and academic software packages are available which provide gaseous properties, including CEA [1], EGlib [2], PEGASE [3], Chemkin [4], MUTATION [5], Cantera [6], and KAPPA [7], however, these libraries tend to focus on a specific application, a narrow range of collisional time-scales, or are specialized in providing only certain types of properties.
These observations have led to the desire to reduce the work necessary to implement new models and algorithms and centralize their development into a single software library which may be used by multiple CFD codes to maximize code reuse, testing, and open collaboration. This paper presents the Mutation ++ library, which has been developed to meet this objective. The library is designed with several goals in mind, including 1. provide accurate thermodynamic, transport, and chemical kinetic properties for multicomponent, partially ionized gases, 2. ensure the efficient evaluation of these properties using state-of-the-art, object-oriented algorithms and data structures in C++, 3. be easily extendable to incorporate new data or algorithms as they become available, 4. interface to any simulation tool based on the solution of conservation laws through a consistent and logical interface, 5. use self-documenting database formats to decrease data transcription errors and increase readability, and 6. be open source to promote code and data sharing among different research communities.
The latest version of Mutation ++ (v1.0.0) has recently been released open-source under the Lesser GNU Public License (LGPL v3) and is freely available on Github 3 . In the remainder of the paper, we present an overview of the library and its impact on the research community to-date. In particular, the four main modules of the library -thermodynamics, transport, chemical kinetics, and gas-surface interaction -are presented, with a few examples to illustrate the library's use.

Generalized Conservation Equations
While it is beyond the scope of this article to describe in detail all the various physicochemical models that are present in the literature, it is useful to briefly present a generalized model which has been used in the design of the library. For a more complete discussion, see the work of Scoggins [8]. We consider a generalized conservation law of the form, where U = ρ i ρu ρE ρẽ m T , is a vector of species mass, momentum, and total and internal energy densities, F (U , ∇ x U ) represents their flux, and S(U ) is a source function. The tilde over the indexed variables in the density vector denotes that these quantities must be expanded over their indices. The exact forms of U , F , and S depend on 1) the coordinate system, 2) the physical model (i.e.: Euler, Navier-Stokes), and 3) the thermochemical model of the gas (i.e.: equilibrium, reacting, multi-temperature, state-to-state). We define the thermochemical state-vector aŝ U = ρ i ρe ρẽ m T where ρe = ρE − ρu · u/2 is the static energy density of the gas. The flux and source functions are closed by constitutive relations for thermodynamic, transport, and chemical properties of the gas. These include quantities such as pressure, enthalpy, viscosity, thermal conductivity, diffusion coefficients, chemical production rates, and energy transfer source terms. In general, these properties are only functions of the local state-vectorÛ and possibly its gradient. This fact allows us to separate the solution of Eq. 1 into two separate domains with limited coupling controlled by the CFD solver and Mutation ++ , as shown in Fig. 1.

Software Architecture
Mutation ++ is designed with a strong focus on Object-Oriented Programming (OOP) patterns in C++. The library's Application Programming Interface (API) is thoroughly documented using the Doxygen format. A continuous integration strategy has been employed. Regression and black box testing are performed through a combination of the Catch2 header-only testing framework and CTest. The primary access to the library is through a Mixture object, which is implemented as a set of submodules encapsulating clearly separated physical quantities as depicted in the simplified Unified Modeling Language (UML) diagram in Fig. 2.

Software Functionalities
Each module in Fig. 2 is described in the following subsections. Specific examples of some of the outputs that the library can provide are given in the Sec. 3.

Thermodynamics
The thermodynamics module provides pure species and mixture thermodynamic quantities, such as enthalpy, entropy, specific heats, or Gibbs free energies. Mixture thermodynamic quantities are derived as sums of pure species properties, weighted by the composition of the mixture. Thermodynamic data for pure species can be found in several references [9][10][11][12][13][14][15][16][17]. Differences exist between each database, such as their format, temperature range of applicability, or degree of nonequilibrium supported. Such differences often sway simulation tool designers to select a single database format to support, or hard-code thermodynamic data directly into their models. This approach makes it difficult to update data as needed, or compare with other tools using different databases.
The Mutation ++ framework provides an abstraction layer which enforces a weak coupling between the concrete thermodynamic database, used for any given set of species, and the computation of mixture thermodynamic quantities. Such a design provides the flexibility to swap out different databases as needed, with minimal effort. The NASA 7-and 9-coefficient polynomial databases [18][19][20][21] and a custom XML format which implements a Rigid-Rotor/Harmonic-Oscillator (RRHO) model are currently implemented.
The NASA format is widely used and provides thermodynamic properties for pure species in thermal equilibrium. The RRHO model is suitable for thermal nonequilibrium calculations. In addition, a new database including more than 1200 neutral and ionized species containing C, H, O, and N is shipped with the library in the NASA 9-coefficient format. The details of this database have been published in [22]. The user can specify the concrete thermodynamic model when creating a mixture (Fig. 3a).
A closely related task to the calculation of thermodynamic properties is the solution of chemical equilibrium compositions. The efficient and robust computation of multiphase, constrained equilibrium compositions is an important topic in several fields, including combustion, aerospace and (bio)chemical engineering, metallurgy, paper processes, and the design of thermal protection systems for atmospheric entry vehicles (e.g., [23][24][25][26][27][28]). Several challenges associated with computing chemical equilibria make conventional methods hard to converge under some conditions [1,[29][30][31][32]. A new multiphase equilibrium solver, based on the single-phase Gibbs function continuation method [33,34], has been developed specifically for Mutation ++ . The Multiphase Gibbs Function Continuation (MPGFC) solver is robust for all well-posed constraints. More details about the solver can be found in [35].

Transport
Closure of transport fluxes is achieved through a multiscale Chapman-Enskog perturbative solution of the Boltzmann equation, yielding asymptotic expressions for the necessary transport coefficients, such as thermal conductivity, viscosity, and diffusion coefficients [36][37][38][39][40]. Explicit expressions for these coefficients are derived in terms of linear transport systems through a Laguerre-Sonine polynomial approximation of the Enskog expansion at increasing orders of accuracy [41][42][43]. These linear systems are functions of trans-  port collision integrals and the local state-vector, and may be solved through a variety of methods.
Collision integrals represent Maxwellian averages of collision cross-sections for each pair of species considered in a given mixture [36,44], weighted depending on the Laguerre-Sonine polynomial order used [45]. The preferred method to compute collision integrals is to numerically integrate from accurate ab initio potential energy surfaces. Such data is available for several important collision systems [46][47][48]. When potential energy surfaces are not available, collision integrals are integrated from model interaction potentials. The evaluation of these model integrals can be partitioned in several ways -neutral-neutral, ion-neutral, electron-neutral, and charged interactions, heavy, electron-heavy, and electron interactions -each with different functional dependencies on temperature and degree of ionization. Mutation ++ introduces a custom XML format for storing collision integral data (Fig. 3b)  Just-in-time loading and efficient evaluation of this data is handled by a CollisionDB object, as shown in Fig. 2.
The solution of the linear transport systems represents a significant CPU time for some CFD applications. Several algorithms have been proposed in the literature for reducing this cost [42,[49][50][51][52]. Mutation ++ provides plug-and-play transport algorithms through the use of selfregistering algorithm classes. For example, the abstract class ThermalConductivityAlgorithm, shown in Fig. 2, provides the necessary interface that all thermal conductivity algorithms must include, namely functions for computing the thermal conductivity and thermal diffusion ratios. Specific algorithms are then implemented by creating a concrete class which implements the interface. This pattern has been used for the calculation of the multicomponent diffusion matrix and shear viscosity as well.

Kinetics
The goal of the chemical kinetics module is the efficient and robust computation of species production rates due to finite-rate chemical reactions. For reaction set R involving species in S, we consider production rates of the forṁ where a full description of each term is given in [8].
Apart from the reaction rate temperatures, knowledge of reaction types is essential in some energy exchange mechanisms.
Manually inputting the type of every reaction in a mechanism of hundreds or thousands of reactions can be a tedious and error-prone process. Therefore, Mutation ++ provides a unique feature which determines the type of reaction automatically when a mechanism is loaded. The problem is formulated as classification tree [53], which can be constructed automatically using simple characteristics of each reaction. An example of such a classification tree is provided in Fig. 4.
In principle, the evaluation of Eq. 2 is straightforward, though great care is required to do it robustly and efficiently. A simplified class diagram of the kinetics module is presented in Fig. 2. The module contains a list of Reaction objects provided by the user through an XML reaction mechanism file (see Fig. 3c). The rest of the module is comprised of a set of computational "managers", which are responsible for the efficient evaluation of individual parts of Eq. 2. These include the evaluation of reaction rates, operations associated with the reaction stoichiometry (the sum and products in Eq. 2), and the evaluation of the third-body term, Θ r . An additional manager class is responsible for evaluating the Jacobian of species production rates, necessary for implicit time-stepping CFD algorithms. Finally, the Kinetics class orchestrates the use of each of these managers to evaluate Eq. 2 and its Jacobian with respect to species densities and temperatures.

Gas-Surface Interactions
The Gas-Surface Interaction (GSI) module provides surface boundary conditions for Eq. 1. They are obtained by applying the conservation of mass, momentum, and energy in a thin controlvolume on a surface at steady-state in the form where F g and F b are gas and bulk phase fluxes, n is the surface normal, and S s is the source term associated with surface processes. Mutation ++ provides several built-in terms that can be mixed to create a model through a custom XML format, as shown in Fig. 3d. Once a model is specified, balance equations are created dynamically through an object-oriented approach forming a Surface object (Fig. 2). The resulting nonlinear equations are functions of the surface thermodynamic state (stored asÛ = ρ i ρe ρẽ m T in SurfaceState), thermochemical properties of the interface, and the two connecting phases (e.g. SurfaceProperties and SolidProperties). This framework provides flexibility for the description of a variety of surfaces (e.g. chemically active, impermeable, porous with fixed outgassing, etc.).
For each type of surface, Mutation ++ provides the fluxes and source terms expressed in Eq. 3 to the client code. A very unique feature of the library is that it, on demand, solves the steadystate balances in a robust and efficient way, to obtain the boundary condition necessary for the CFD or material solvers. More information about the GSI models available can be found in [54].

Thermodynamic and Transport Properties
Fig. 5 presents a selection of thermodynamic and transport properties computed by Mutation ++ for an 11-species, isobaric air mixture in thermochemical equilibrium. The equilibrium mole fractions and thermodynamic properties are provided using both the NASA-9 and RRHO thermodynamic databases. Where possible, comparisons with equilibrium air curve-fits of D'Angola et al. [55] and thermal conductivity data from Murphy [56] and Azinovsky et al. [57] are also shown. Both "frozen" and "equilibrium" curves are shown for the specific heat at constant pressure and specific heat ratio. The frozen curves neglect the dependence of the species composition on the temperature through the equilibrium reactions, while equilibrium curves do not. This distinction is important, as it shows that these properties can vary substantially, depending on the thermochemical model employed by the user. For more information regarding the models, data, algorithms, and interpretation of these figures, please see the discussion in [58].

Equilibrium Ablation Rates
An important problem in the prediction of material response for thermal protection systems (TPS) of atmospheric entry vehicles is the solution of so-called "B-prime" tables. B-prime tables describe the equilibrium gas composition at the surface of an ablating TPS as well as its mass loss rate due to reactions at the surface, such as oxidation or nitridation. Assuming a thin control volume over an ablating TPS with equal diffusion coefficients for each species, conservation of elements inside the control volume yields where y is the elemental mass fraction for any element in the mixture, the subscripts w, c, g, and e refer to wall, char, pyrolysis gas, and boundary layer edge properties respectively, B ≡ m/(ρ e u e C M ) is a mass blowing rate, nondimensionalized by the boundary layer edge mass flux, and C M is the local Stanton number for mass transfer. When coupled with the minimization of Gibbs energy at a known surface condition, Eq. 4 can be solved to obtain species composition and char mass blowing rate B c . Fig. 6 shows such a calculation performed using Mutation ++ with the custom thermodynamic database discussed in Sec. 2.3.1. The B c results are compared with those obtained with the thermodynamic database from the NASA Chemical Equilibrium with Applications (CEA) code. More discussion on the differences in these databases can be found in [22].

Impact
The main goal of Mutation ++ is to promote open collaboration between different research groups and communities working in the broad areas of hypersonics, combustion, and plasma physics, by lowering the cost of developing new physicochemical models, data, and algorithms, for modeling gas and gas-surface phenomena. With an efficient and extensible frame-work, Mutation ++ can be easily coupled with existing CFD tools, allowing researchers to test effectively new thermodynamic, transport, or chemical models, physicochemical data, or numerical algorithms. In addition, users of the library can benefit from the work of others through collaborative testing, bug fixing, and maintenance which is supported by a continuous development and integration strategy.
Since its creation, Mutation ++ has been used in many diverse applications, branching out from the original motivation of hypersonic flows for atmospheric entry [59][60][61][62]. These include the study of biomass pyrolysis [63], solar physics, magnetized transport [43], and meteor phenomena [64,65]. Application of the library to fields other than originally intended serves to highlight the extensibility of the framework and the impact it can have on basic research. The library has also been used in limited commercial settings. Recently, Mutation ++ has been coupled to the material response codes Amaryllis of LMS Samcef (Siemens) and the Porous material Analysis Toolbox (PATO), developed jointly between C la Vie and NASA Ames Research Center.

Conclusions
The Mutation ++ library provides an OOP framework for computing thermodynamic, transport, and kinetic properties of non to fully ionized gas mixtures at all points on thermochemical nonequilibrium spectrum. The library leverages the simple dependence of thermochemical properties on the local thermodynamic state of the gas to implement a weak coupling between the computation of those properties and the simulation tools which need them, through a clean and consistent API.
The code is freely available on GitHub with an LGPL v3.0 license, following a continuous integration development strategy with periodic versioning to alleviate backward compatibility concerns for its users. This paper marks version v1.0.0 for the library. Future versions will aim to provide additional features and greater flexibility for the end-user.    Figure 6: Equilibrium char mass blowing rate and wall mole fractions for a carbon-phenolic ablator in air computed with Mutation ++ [22].