Calculation of Reaction Extent as a Function of Pressure and Temperature: Issue with Activity-Based Approach

Question

I'm working on modeling reactions involving combinations of hydrogen (H), nitrogen (N), and strontium (Sr) elements. I've obtained the temperature-dependent Gibbs free energy of formation, enthalpy, and entropy for all the species involved in these reactions. My goal is to determine the extent of the reaction (denoted as x) as a function of pressure (p) and temperature (T) using the activity of these species.

In my model, all species containing Sr* are considered solid, and I've set their activity to one. For gas-phase species, I'm calculating the activity using the following expression:

$$ \text{activity(p, x)}_\text{gas} = \left( \frac{p}{p_{\text{ref}}} \cdot \frac{n_{i,0} + \text{exponent_sign} \cdot v_i \cdot x}{\sum n_i + x \cdot \sum v_i} \right)^{v_i \cdot \text{exponent_sign}} $$

Here, $v_i$ represents the stoichiometric coefficient of species, $n_{i,0}$ is the initial mole amount of species, $\sum n_i$ is the total initial mole amount, $\sum v_i$ is the total stoichiometric coefficient, and exponent_sign is negative for reactants and positive for products.

To calculate the Gibbs energy of the reaction: $$ G(T, p, x) = G_{\text{formation}}(T) + R \cdot T \cdot \ln{\left( \prod_{i} a_i(x, p) \right)} $$

However, when I attempt to minimize this expression, I'm obtaining reaction extents of either 1 or 0, but not values in between. I suspect there might be an error in my approach, or perhaps I'm overlooking something. Can anyone offer insights into what might be causing this issue? Alternatively, if my approach is incorrect, I would appreciate suggestions on how to accurately calculate the reaction extent as a function of pressure and temperature using the available data.

Any assistance or guidance on this matter would be greatly appreciated! Thank you in advance for your help :)

Edit:

• I'm currently confused about calculating the total mole amount using: $$ \sum {n_i}, $$ Should this summation include both species in solid and gas phases, or only species in gas phase? Similarly, I'm uncertain about:$$ \sum {v_i}, $$ In my understanding, both terms are meant to reflect "Le Chatelier's Principle." Do solids have an effect on this summation, or should it only consider gas species?
• If the function values of $G(T, p, x)$ for all reaction extend $x$ are negative respective positive, what does that mean?

Answer 1

Since you asked for general assistance and I have dealt with these type of calculations, I will give you some general points as the question evolved these last days:

• The expression for the activity of a gaseous species is strange. I imagine that a single species participates in various reactions, not one. Thus, there should be several extents of reactions in the numerator. Let $i$ denote the chemical species for which we have $i = 1,2,...,M$, and $j$ the chemical reaction for which we have $j = 1,2,...,N$. Then, the molar fraction of species $i$ in the gas phase is \begin{equation} y_i = \frac{n_i}{n} = \frac{n_{i0} + \sum_{j = 1}^N \nu_{ij}\xi_{j}} {\sum_{i = 1}^M \left[n_{i0} + \sum_{j = 1}^N \nu_{ij}\xi_{j}\right]} \tag1 \end{equation} where the brackets are not really needed, but are added for clarity.
• The gas composition $y_i$ is a variable that belongs to the gas phase. Thus, the total amount in the denominator of Eq. (1) only includes the amounts of the species that 'live' in the gaseous phase.
• The activity for the solid species are one, but care must be taken once you solve the system and hopefully reach a solution. Once you have all the extents of reactions $\xi_1,\xi_2,...,\xi_N$, you calculate the final amounts for all the solid species. If any of those are zero or negative, then something is wrong in the calculation or the chosen solution is wrong. Remember that a set of non-linear equations yields many solutions, but only one in this case will be physically meaningful.
• The approach that you are taking which is to minimise the total Gibbs energy is not correct, and the equation you are minimising is wrong. Try to solve the $N$ equations which appear with the $N$ laws of mass action \begin{equation} K_j = \prod a_{ij}^{\nu_{ij}} \quad j = 1,2,...,N \tag2 \end{equation} If you want to have the value of the total Gibbs energy, then this value is \begin{equation} G = \sum_{i = 1}^M n_i G_i = \sum_{i = 1}^M n_i [\Delta_f G_i + RT\ln a_i] \tag3 \end{equation}
• The value of Eq. (3) is the minimum value and has units of $\ce{J}$. It can be positive or negative.
• If you insist on the approach of minimising $G$, i.e. Eq. (3), you can extend this method using the Lagrange multipliers. It can be proved that you will have to solve even more equations, a total of $M$ (number of species) and $w$ (number of material balances). The number $w$ coincides with the number of Lagrange multipliers $\lambda$'s which are not known. Search online and you will get many papers that tried this idea. However, keep it simple and just do it like in class like I said.