Theory

The problem of determining the equilibrium concentration of various species given a set of input elements can be formulated in terms of either entropy, Gibbs or Helmholtz function. For example, if the system is defined in terms of temperature and pressure, then minimisation of Gibbs function is appropriate since temperature and pressure are independent variables of Gibbs function [Kenneth1956] . If the system is defined in terms of temperature and volume, then minimisation of Helmholtz function is appropriate.

For a multi-phase system with \(\mathrm{(x_1, x_2,\ldots, x_{N_g},\ldots x_{N_g + N_s})}\) moles of \(\mathrm{(N_g + N_s)}\) species containing \(\mathrm{N_e}\) elements, the Helmholtz function to be minimized is as follows:

(1)\[\begin{split}\mathrm{A_{system}(T, V)} & = \mathrm{G - PV}\\ & = \sum_{i = 1}^{\mathrm{N_s + N_g}} \mu_i x_i - \mathrm{PV}\end{split}\]

Where G is the Gibbs function. The system pressure P is in bars. Dividing Eq. (1) by RT, we get,

(2)\[\mathrm{\tilde{A}(T, V)} = \frac{\mathrm{A_{system}}}{\mathrm{RT}} = \underbrace{\sum_{i = 1}^{\mathrm{N_s + N_g}}\tilde{\mu}_i x_i - \bar{x}^g}_{\mathrm{F(x)}}\]

Where \(\bar{x}^g\) is the total number of moles of gas in the system. \(\tilde{\mu}_i\) is the reduced (dimension less) chemical potential and can be given in terms of their standard chemical potential functions as,

(3)\[\begin{split}\tilde{\mu}_i = \begin{cases} {\left(\tilde{\mu}^o\right)}_i^g + \ln\mathrm{\frac{RT}{V}} + \ln x_i^g & \text{for }i = 1, 2 \ldots, \mathrm{N_g} \\ {\left(\tilde{\mu}^{o}\right)}_{i}^{c} & \text{for }i = 1,2\ldots,\mathrm{N_s} \end{cases}\end{split}\]

The superscript g and c are used for gaseous and condensed phase species. For example, \(x_{i}^{g}\) is the number of moles for \(i^{th}\) chemical species in gas phase. \({\left(\tilde{\mu}^o\right)}_{i}^{g}\) is chemical potential of \(i^{th}\) gaseous chemical species in gas phase at standard conditions. \(\mathrm{N_s}\) is the total number of the condensed species, \(\mathrm{N_g}\) is the total number of the species in the gaseous phase.

While the free energy Eq. (2) is minimized to solve for equilibrium, the species have to satisfy non-negativity constraint \(x_i\ge 0\) and the constraint for element conservation given as,

(4)\[\begin{align} \underbrace{\sum_{i = 1}^{\mathrm{N_g}} a_{ij}^{g} x_{i}^{g} + \sum _{i = 1}^{\mathrm{N_s}} a_{ij}^{c}x_{i}^{c} = b_j}_{\mathrm{C_j(x)}} &&j = 1, 2, 3\ldots, \mathrm{N_e} \end{align}\]

Where \(a_{ij} ^{g}\) is the number of atoms of \(j^{th}\) element in specie \(i\) in the gas phase and \(a_{ij}^{c}\) is the number of atoms of \(j^{th}\) element in \(i^{th}\) species with condensed phase. \(b_j\) is the total number of moles of element j, originally present in system mixture. The above free energy functions are minimised with two methods, viz., (1) Quadratic gradient descent method (2) sequential least square minimisation (SLSQP).

Quadratic gradient descent method:

To minimise the free energy (Eq. (2)) of the system containing \(\mathrm{(x_1, x_2,\ldots x_{N_g}, \ldots,x_{N_g + N_s})}\) moles of \(\mathrm{N_g}\) + \(\mathrm{N_s}\) species, with constraints described in Eq. (4), the method of Lagrange \('\) s underdetermined multipliers [Bertsekas2014] [Weisstein2020] is used. Here, the formulation is given for constant temperature and constant volume problem. Formulation for constant pressure and constant temperature can be found in literature [White1958] [McBride_and_Gordon1996a] [Eriksson1979] [Eriksson_and_Hack1990]. The Lagrangian function to be minimized can be written as,

(5)\[\mathrm{L} = \mathrm{\tilde{A}} - \sum_{j=1}^{N_e}\pi_j C_j(x)\]

Where, \(\pi_j\) are the undetermined multipliers. The partial derivatives of L with respect to \(i^{th}\) chemical species can be given after regrouping in terms of the gas species part and condensed species part as,

(6)\[\frac{\partial L}{\partial x_i} = \mathrm{\frac{\partial L^g} {\partial x_i} +\frac{\partial L^c}{\partial x_i}}\]

Where, the \(\mathrm{\frac{\partial L^g}{\partial x_i}}\) and \(\mathrm{\frac{\partial L^c}{\partial x_i}}\) are given by,

(7)\[\begin{aligned} \mathrm{\frac{\partial L^g}{\partial x_i}} = (c_i + \ln x_i^g) - \sum_{j=1}^{\mathrm{N_e}}a_{ij}^{g}\pi_j = 0 && i = 1,2\ldots,\mathrm{N_g} \end{aligned}\]

(8)\[\begin{aligned} \mathrm{\frac{\partial L^c}{\partial x_i}} = {\left(\tilde{\mu}^o\right)}_i^c - \sum_{j=1}^{\mathrm{N_e}}a_{ij}^{c} \pi_j = 0 &&i= 1,2\ldots,\mathrm{N_s} \end{aligned}\]

Where,

(9)\[c_i = {\left(\tilde{\mu}^o\right)}_i^g + \ln\mathrm{\frac{RT}{V}}\]

To linearize Eq. (7), Taylor series expansion about \(y_{i}^{g}\) is carried out, which results,

(10)\[\begin{align} f_i + \frac{x_i^g}{y_i^g} - 1 - \sum_{j=1}^{\mathrm{N_e}}a_{ij}^{g}\pi_j = 0 &&i = 1,2\ldots,\mathrm{N_g} \end{align}\]

Where, \(f_{i}\) can be given as,

(11)\[f_i = c_i + \ln y_i^g\]

From the Taylor expanded form, the improved mole numbers for the next iteration can be obtained from the rearranged Eq. (10)

(12)\[x_i^g = -f_i y_i^g + y_i^g \left(\sum_{i=1}^{\mathrm{N_g}} \pi _{j} a_{ij}^{g} + 1\right)\]

Substituting Eq. (12) in Eq. (4) we have,

(13)\[\begin{align} \sum_{k=1}^{\mathrm{N_e}} r_{jk} \pi_k + \sum_{i=1}^{\mathrm{N_s}} a_{ij}^{c}x_{i}^{c} & = b_j + \sum_{i=1}^{\mathrm{N_g}} a_{ij}^{g} f_i y_i^g - \sum_{i=1}^{\mathrm{N_g}}a_{ij}^{g}y_i^g &&j= 1,2\ldots, \mathrm{N_e} \end{align}\]

Where,

(14)\[\begin{align} r_{jk} & = r_{kj} = \sum_{i=1}^{\mathrm{N_g}} (a_{ij}^{g}a_{ik}^{g})y_{i}^{g} \text{ }j = 1, 2\ldots, \mathrm{N_e} \end{align}\]

The Eqs. (8) and (13) are solved simultaneously to get \(\pi _j\) and \(x_i^c\). Using \(\pi _j\), updated values of \(x^g\) are obtained using Eq. (12). Here, it should be noted that the formulation of condensed species is such that, the moles of each condensed phase species is directly obtained from the solution of Eq. (8) and (13), without applying any correction. If all \(x^g\) values are positive, they are considered as the guessed value for the next iteration. If not, then guessed values are corrected with the following Eq.as [Gunnar_Eriksson1971],

(15)\[y_{i, improved}^{p} = y^{p}_{i} + \lambda (x^{p}_{i} - y^{p}_{i})\]

Where, p is the phase of chemical species (gaseous species g or condensed species c). The \(\lambda\) is the correction factor, and can be given according to [Gunnar_Eriksson1971],

(16)\[\lambda = 0.99 \lambda ' (1 - 0.5 \lambda ')\]

Where, \(\lambda '\) is the value required for the next step to remain positive as given by,

(17)\[\lambda ' = \min_{i}\left(\frac{y^{g,c}_{i} }{(y^{g,c}_{i} - x^{g,c}_{i})} \right)\]

The corrected values of \(y^g\) are considered for the next iteration. Since, in the above formulation, the set of condensed species is not known beforehand, for the first iteration, only gaseous species equilibrium is obtained. In subsequent iterations, species are added such that they reduce overall system \('\) s free energy. This can be assured by [McBride_and_Gordon1996a] [Gunnar_Eriksson1971],

(18)\[{\left(\tilde{\mu}^o\right)}_i^c - \sum_{j=1}^{\mathrm{N_e}}\pi_{j}b_{j} \le 0\]

Sometimes with particular species, the combination can lead to a set of the dependent Eqs.in the formulation. Set of dependent species are identified with the help of reduced row echelon form. Subsequently, dependent species, as well as the matching combination of dependent species present in species set, are removed temporarily, and only species which decrease overall free energy of the system is included in species set. For example, species set containing, U (L), \(\mathrm{UO_{2}(L)}\) and \(\mathrm{U_4O_9(III)}\), leads to linear dependence and all three species are removed, and species which lower free energy of the system are included in species set.

The above Helmholtz function minimization subjected to mole number conservation constraint is implemented in python and the code is hereafter referred to as MINICHEM (MINImisation of CHEMical potentials). Once the equilibrium species moles in various phases are determined, the release fraction of the \(\mathrm{j^{th}}\) chemical species are determined as follows:

\[\begin{aligned} \text{Elemental release fraction = }\mathrm{{RF(e)}_j} = \frac{{\sum\limits}_{i} \text{Number of moles of \(j^{th}\) element in \(i^{th}\) gaseous species}}{\text{Total mole inventory of \(j^{th}\) element}} \end{aligned}\]

The isotopic release fraction is defined as,

\[\text{Isotopic release fraction }{RF(iso)}_k = {\mathrm{RF(e)}_j} \times \text{isotopic fraction of \(k^{th}\) isotope}\]