Chemical equilibrium calculations using SCILAB

Introduction

Many physico-chemical phenomena occurring in transformations of raw materials, includes chemical reactions. Some of them are reversible, which means that, they can reach an equilibrium (in thermodynamical way). The knowledge of the mechanisms and principles of the chemical equilibria is very important for any chemical engineer.

One of the most common chemical reaction of industrial importance is the gas phase synthesis of methanol from syngas (carbon monoxide and hydrogen) on catalyst. This reaction is reversible and essentially exothermic since its reaction enthalpy is negative along a broad range of temperature.

In this article we will solve a problem related to chemical reaction equilibrium. The solution includes finding a root of a nonlinear algebraic equation, which is formulated using Guldberg-Waage’s and Van’t Hoff laws. Due to nonlinear nature of the equation, the application of computational software is necessary. In this case, computer program based on Newton-Raphson method in SCILAB was used.

Definition of the problem

We will use a problem from “Chemical reactor analysis and design” by J.B. Rawlings, namely the exercise 3.8 on page 66. This problem, which concerns the methanol synthesis, is defined below.

Methanol can be manufactured by the gas-phase reaction

$CO+2H_{2}\leftrightharpoons CH_3OH$

A batch reactor is charged with H₂ and CO mixture at 1.0 atm and the pressure is maintained at 1.0 atm. The initial mole fraction of H₂ is 0.6 and CO is 0.4. At 400K,

$\Delta G^o=-333\text{ }cal/mol$

$\Delta H^o=-22,580\text{ }cal/mol$

Determine the temperature T at which the equilibrium mole fraction of H₂is 0.05.

Solution

Our goal is to formulate an algebraic equation with one unknown, which is the temperature required to obtain assumed molar fraction of hydrogen. In calculations we have accepted ideal gas law and constant reaction enthalpy through the temperature range.

Step 1. We need to find the limiting reactant in the reaction between carbon monoxide and hydrogen. To do that, we will compare ratios:

$y_{i0}/ v_i$

where:

y_i₀ – initial molar fraction of i-th reactant,

v_i – stoichiometric coefficient in the chemical reaction of i-th reactant.

Since 0.4/1 > 0.6/2, the limiting reactant is hydrogen.

Step. 2. Find expressions for final mole fractions of reactants

Let’s use symbol x as the molar advancement of the reaction, and n₁, n₂, n₃as the numbers of moles of carbon monoxide, hydrogen and methanol, respectively. Using such notation, the number of moles of reactants will be described by the system of equations:

(1) $\begin{equation*} n_1=n_{10}-x \end{equation*}$

(2) $\begin{equation*} n_2=n_{20}-2x \end{equation*}$

(3) $\begin{equation*} n_3=x \end{equation*}$

The total number of moles is :

(4) $\begin{equation*} n=n_{0}-2x \end{equation*}$

Since the limiting reactant is hydrogen, we will introduce a new variable called the fractional conversion X of hydrogen, as :

(5) $\begin{equation*} X=2x/n_{20} \end{equation*}$

According to the definition, fractional conversion is the number of moles of a compound that reacted divided by the amount of the moles that were fed.

After introducing the fractional conversion X into Eq. (1), Eq. (2), Eq. (3) and Eq. (4), we get:

(6) $\begin{equation*} n_1=n_{10}-(X\cdot n_{20})/2 \end{equation*}$

(7) $\begin{equation*} n_2=n_{20}-X\cdot n_{20} \end{equation*}$

(8) $\begin{equation*} n_3=(X\cdot {n_20})/2 \end{equation*}$

(9) $\begin{equation*} n=n_0-X\cdot n_{20}=n_{20}\cdot (1/y_{20} -X) \end{equation*}$

with $y_{20}=0.6$

The mole fractions of components 1, 2, 3 are respectively:

(10) $\begin{equation*} y_1=(y_{10}/y_{20} -X)/(1/y_{20} -X) \end{equation*}$

(11) $\begin{equation*} y_2=(1-X)/(1/y_{20} -X) \end{equation*}$

(12) $\begin{equation*} y_3=(0.5\cdot X)/(1/y_{20} -X) \end{equation*}$

Step 3. Find the equilibrium fractional conversion X_eq matching the final mole fraction of hydrogen y₂=0.05.

(13) $\begin{equation*} \frac{(1-X)}{(1/y_{20} -X)}=0.05\Rightarrow X_{eq}=\left(1-\frac{1}{y_{20}}\cdot y_2\right)/(1-y_2 ) \end{equation*}$

Step 4. Find the expressions of equilibrium constant K_P

Using the Mass Action Law (MAL) or Waage-Guldberg’s law, we obtain the equation for K_p as:

(14) $\begin{equation*} K_p=\frac{(y_3\cdot P^{(1-1-2)})}{(y_1\cdot y_2)}=\frac{0.5\cdotX\cdot(1/y_{20} -X)}{(1-X)\cdot (y_{10}/y_{20} -0.5\cdotX)\cdot P^2} \end{equation*}$

with P=1.0 atm

Equilibrium constant K_P can be also determined by the Van’t Hoff equation. Since the reaction enthalpy the integration of the Van’t Hoff equation leads to the following equation:

(15) $\begin{equation*} K_p=exp⁡\left(\frac{-∆G^0}{R\cdot T_0}\right)\cdot exp⁡\left(\frac{(-∆H^0)}{R} \left(\frac{1}{T}-\frac{1}{T_0} \right)\right)\end{equation*}$

with T₀=400 K and the universal gas constant R=1.97 cal/(mol K).

Step 5. Obtain final equation for temperature calculation

By equating the equations (14) and (15), we get the final equation:

(16) $\begin{equation*} \frac{(0.5\cdot X\cdot (1/y_{20} -X_{eq}))}{(1-X_{eq} )\cdot (y_{10}/y_{20} -0.5\cdot X_{eq} )\cdot P^2 }=exp⁡\left(\frac{-∆G^0}{R\cdot T_0}\right)\cdot exp⁡\left(\frac{(-∆H^0)}{R} \left(\frac{1}{T}-\frac{1}{T_0} \right)\right) \end{equation*}$

where the unknown is the temperature T required to obtain the final mole fraction of hydrogen 0.05 (or the equilibrium fraction conversion of hydrogen).

Implementation in SCILAB

First, let’s define parameters occurring in the problem. We will use alphanumerical notation, so subscript ‘i0’ will denote mole fraction of i component at initial time.

function equilgasphasereact
    clc
    disp('this code allows to solve a typical problem of chemical engineering, notably that of chemical thermodynamical equilibrium, molar advancement')
    disp('data')
    dg0=-333 //reaction enthalpy
    dh0=-22580 // reaction enthalpy at 400 K
    T0=400// reference temperature
    R=1.97 //universal gas constant
    y2=0.05// equilibrium mole fraction of H2
    y10=0.4 // initial mole fraction of CO
    y20=0.6// initial mole fraction of H2
    P=1 // total pressure inside reactor
    disp('assumptions: 1) ideal gas law, 2) constant reaction enthalpy dH0')

The fractional conversion of hydrogen is calculated according to Eq (13).

    disp('the equilibrium fractional conversion of hydrogen for y2=0.05')
    Xeq=(1-5/3*y2)/(1-y2)
    disp(Xeq)

Let’s use Eq. (16) to calculate the temperature required to obtain assumed molar fraction of hydrogen. To this end we will write the Eq. (16) as the function g(T) in the SCILAB syntax, as below.

    disp('solve the equation to find the temperature T (K) for which y2=0.05')
    deff('y=g(T)','y=0.5*Xeq*(5/3-Xeq)/((1-Xeq)*(2/3-0.5*Xeq))-exp(333/(1.97*400))*exp(22580/1.97*(1/(T)-1/400))')

The Eq. (16) is then solved through the use of the function « fsolve » that runs the Newton-Raphson algorithm.

    [Teq fval]=fsolve(200,g)
    disp('the temperature Teq is'), disp(Teq)

We will calculate molar fractions of other components according to Eq. (10), Eq. (11) and Eq. (12).

    disp('the final mole fractions of carbon monoxide, hydrogen and methanol are resp')
    y1=(y10/y20-0.5*Xeq)/(1/y20-Xeq); y2=(1-Xeq)/(1/y20-Xeq); y3=0.5*Xeq/(1/y20-Xeq)

In order to verify the correctness of the calculations we should check if the sum of molar fractions is equal one, as in the code below:

    disp([y1 y2 y3])
    disp('verification of computations by the sum of mole fractions'), disp(sum([y1 y2 y3]))
endfunction

The code presented above allows to solve a problem of chemical thermodynamical equilibrium by combining Guldberg-Waage’s law and Van’t Hoff law.

The output of the program is following:

The equilibrium fractional conversion of hydrogen for y2=0.05

0.965

solve the equation to find the temperature T (K) for which y2=0.05

the temperature Teq is

356.06 K=82.91 °C

the final mole fractions of carbon monoxide, hydrogen and methanol are resp y1, y2 and y3

0.2625 0.05 0.6875

verification of computations by the sum of mole fractions
1.

Summary

In this tutorial we have used SCILAB software to solve a problem of chemical thermodynamical equilibrium. Synthesis of methanol, which is an important industrial reaction was used as an example. From mathematical point of view it was necessary to find a root of nonlinear algebraic equation.

We saw that most of equations in chemical engineering are nonlinear due to the kinetics, mass action law, temperature dependence (Arrhenius law) etc., so, this requires the use of software to solve these equations.

Written by Etonde Yannick

I’m a Cameroonian process engineer and the head of the process engineering department of the Saint Jerome Catholic University of Douala (Cameroon). My areas of interest are: process modeling, simulation and optimization, chemical reaction engineering, unit operations, design and scale-up, waste treatment, wastewater treatment and QHSE. I’m a PhD student. I’m working on the modeling, simulation and optimization of the pyrolysis of waste plastics using the Model Predictive Control.