Calculation of Chemical Equilibria in Multi-Phase: Multicomponent Systems

Abstract

The paper describes methods for calculating chemical equilibria based on a constrained Gibbs free energy minimization. The methods allow the treatment of multicomponent systems with multiple phases, including gaseous phases, condensed phases, and stoichiometric phases. A special aspect is the detection and treatment of miscibility gaps. The underlying mathematical problem is described in detail together with the algorithmic approach for its solution. Results are presented for some test cases, including the computation of phase diagrams for ternary systems.

Appendix: The Gibbs Free Energy Function

Each phase specific Gibbs free energy function G ⁽ⁱ⁾ is the sum of a reference part G _ref ⁽ⁱ⁾, an ideal part G _id ⁽ⁱ⁾, and an excess part G _ex ⁽ⁱ⁾:

$$ \displaystyle\begin{array}{rcl} G^{(i)}(x^{(i)})&:=& G_{\mathrm{ ref}}^{(i)}(x^{(i)}) + G_{\mathrm{ id}}^{(i)}(x^{(i)}) + G_{\mathrm{ ex}}^{(i)}(x^{(i)}).{}\end{array}$$

(16)

The reference and the ideal part are given by:

$$\displaystyle\begin{array}{rcl} G_{\mathrm{ref}}^{(i)}(x^{(i)})&:=& \sum _{ k=1}^{K^{i} }G_{k}^{0,i}x_{ k}^{i},{}\end{array}$$

(17)

$$\displaystyle\begin{array}{rcl} G_{\mathrm{id}}^{(i)}(x^{(i)})&:=& RT\sum _{ k=1}^{K^{i} }x_{k}^{i}\ln x_{ k}^{i}.{}\end{array}$$

(18)

Here, R = 8.31451 J/(mol K) is the ideal gas constant and T the temperature in K. In the non-ideal part, binary interactions of phase constituents are described by the terms B ⁱ, ternary interactions by the terms T ⁱ and quaternary interactions by the terms Q ⁱ (for a phase i):

$$\displaystyle\begin{array}{rcl} G_{\mathrm{ex}}^{(i)}(x^{(i)})&:=& \sum _{ k_{1}<k_{2}}^{K^{i} }B^{i}(\{k_{ 1},k_{2}\},x^{(i)}) +\sum _{ k_{1}<k_{2}<k_{3}}^{K^{i} }T^{i}(\{k_{ 1},k_{2},k_{3}\},x^{(i)}) \\ & +& \sum _{k_{1}<k_{2}<k_{3}<k_{4}}^{K^{i} }Q^{i}(\{k_{ 1},k_{2},k_{3},k_{4}\},x^{(i)}), {}\end{array}$$

(19)

with

$$\displaystyle\begin{array}{rcl} B^{i}(\{k_{ 1},k_{2}\},x^{(i)})&:=& \prod _{ j=1}^{2}x_{ k_{j}}^{i}\sum _{ \nu =1}^{m^{i}(k_{ 1},k_{2})}L_{k_{ 1}k_{2}}^{\nu,i}\left (x_{ k_{1}}^{i} - x_{ k_{2}}^{i}\right )^{\nu -1},\ m^{i}(k_{ 1},k_{2}) \in \mathbb{N},{}\end{array}$$

(20)

$$\displaystyle\begin{array}{rcl} T^{i}(\{k_{ 1},k_{2},k_{3}\},x^{(i)})&:=& \prod _{ j=1}^{3}x_{ k_{j}}^{i}\left (\sum _{ j=1}^{3}L_{ k_{j}}^{i}x_{ k_{j}}^{i} + \frac{1} {3}\sum _{j=1}^{3}L_{ k_{j}}^{i}\left (1 -\sum _{ j=1}^{3}x_{ k_{j}}^{i}\right )\right ),{}\end{array}$$

(21)

$$\displaystyle\begin{array}{rcl} Q^{i}(\{k_{ 1},k_{2},k_{3},k_{4}\},x^{(i)})&:=& \prod _{ j=1}^{4}x_{ k_{j}}^{i}L_{ k_{1}k_{2}k_{3}k_{4}}^{i}.{}\end{array}$$

(22)

The binary excess part of the phase specific Gibbs function (20) is described here by the Redlich–Kister–Muggianu polynomials (RKMP), cf. [ 15 ]. Interactions of higher order could also be involved in the model but are neglected here.

If K ⁱ = 1 for some i ∈ {1, …, n}, it holds:

$$\displaystyle{G^{(i)}(x^{(i)}) = G_{ 1}^{0,i},}$$

(Gibbs energy for a stoichiometric phase).

The coefficients G _k ^0,i,  i = 1, …, n,  k = 1, …, K ⁱ are given by:

$$\displaystyle\begin{array}{rcl} G_{k}^{0,i}(T)& =& A_{ k}^{i} + B_{ k}^{i}T + C_{ k}^{i}T\ln T + D_{ k}^{i}T^{2} + E_{ k}^{i}T^{3} + F_{ k}^{i}T^{-1} \\ & +& G_{k}^{i}T^{i_{G} } + H_{k}^{i}T^{i_{H} } + I_{k}^{i}T^{i_{I} } + J_{k}^{i}T^{i_{J} }, {}\end{array}$$

(23)

where $A_{k}^{i},\ldots,M_{k}^{i};\ i_{G},\ldots,i_{M} \in \mathbb{R},\ i = 1,\ldots,n,\ k = 1,\ldots,K^{i}$ are coefficients taken from a material database [4]. The pressure dependence of the Gibbs energy can be involved by additional additive terms.

The RKMP coefficients $L_{(k_{1},\ldots,k_{I})}^{i},\ i = 1,\ldots,n,\ k_{1},\ldots,k_{I} \in \{ 1,\ldots,K^{i}\},k_{1} <\ldots <k_{I}$, where I is the interaction order, are given by:

$$\displaystyle\begin{array}{rcl} L_{(k_{1},\ldots,k_{N})}^{i}(T,P)& =& (A_{ k}^{i})_{ L} + (B_{k}^{i})_{ L}T + (C_{k}^{i})_{ L}T\ln T + (D_{k}^{i})_{ L}T^{2}{}\end{array}$$

(24)

$$\displaystyle\begin{array}{rcl} & +& (E_{k}^{i})_{ L}T^{3} + (F_{ k}^{i})_{ L}T^{-1} + (G_{ k}^{i})_{ L}P + (H_{k}^{i})_{ L}P^{2},{}\end{array}$$

(25)

where the coefficients $(A_{k}^{i})_{L},\ldots,(H_{k}^{i})_{L} \in \mathbb{R},\ i = 1,\ldots,n,\ k = 1,\ldots,K^{i}$ were taken from a material database [4] as well. The pressure dependence can be involved by additional additive terms in this case as well.

Please note that the term lnP has to be added in the case of the presence of a gaseous phase g:

$$\displaystyle{ G^{(g)}(x^{(g)}) = G_{\mathrm{ ref}}^{(g)}(x^{(g)}) + G_{\mathrm{ id}}^{(g)}(x^{(g)}) + RT\ln P. }$$

(26)