Gibbs free energy minimization method can be used to predict equilibrium compositions in reacting systems involving a large number of species. It is being increasingly preferred to the traditional mass balance method as identification of all the reactions occurring in the mixture along with their extents becomes very difficult due to the large number of species involved as cited by Y Lwin [1]. If the minimum value of Gibbs free energy of a system is predicted correctly, the composition at this value would give the equilibrium composition since Gibbs free energy of any system is minimum at equilibrium.

This minimization is carried out with the help of Solver feature present in the Microsoft Excel spreadsheet package [1]. The equilibrium temperature is estimated by carrying out an enthalpy balance across the reacting system.

The case in consideration is the Claus reaction furnace in the SRU operating on the Claus process [2]. Claus process is the most widely employed industrial process to recover elemental sulphur from gaseous hydrogen sulphide (H₂S). Recovery of sulphur is necessitated because of environmental regulations on emissions for several petroleum products mainly, gasoline, kerosene and diesel [3,4]. Moreover, recovery of sulphur helps in increasing the revenue in view of the wide range of applications of elemental sulphur. Hydro treatment converts the sulphur in petroleum products into H₂S. Acid gas (containing H₂S) is separated from hydrotreated gases by amine treatment. Moreover, H₂S bearing Sour Water coming from different other sources is stripped of its H₂S in Sour Water Stripper (SWS). Both these H₂S containing gases are combined as feed to the SRU.

Claus process involves oxidation of a portion of H₂S in the feed gas in the Claus furnace to produce sulphur dioxide (SO₂). The feed is heated to around 200°C before being fed to the Claus furnace where it reacts with air preheated to the same temperature. Oxidation of 1/3rd of the total H₂S takes place at a pressure of 1-2 bar in the furnace to produce SO₂. The exothermic reaction raises the temperature of the furnace to around 1200°C. The resulting gas mixture is then passed through a series of catalytic Claus reactors where elemental sulphur is formed. The sulphur obtained is condensed and extracted as product.

The major reactions taking place in the Claus furnace are [5].

In addition to the above, many other side reactions take place in the furnace. For example, side reactions involving hydrocarbons and CO₂ present in acid gas can result in the formation of carbonyl sulphide (COS), carbon disulphide (CS₂) and carbon monoxide (CO) in the furnace.

The reacting system comprises of the feed (Acid gas & SWS gas) and combustion air, at an initial temperature of about 200°C and pressure of 1-2 bar. The feed consists primarily of H₂S (from Acid gas & SWS gas), N₂, and O₂ (from air). Besides these, H₂O, CO₂, and NH₃ are also present along with small amounts of light paraffinic hydrocarbons (C1, C2, C3, and C4), Hydrogen and Argon. The reactions yield SO₂, S₂ as the main products along with CO, COS and CS₂.

Gibbs free energy of formation (ΔG^T_f) of any species in the system at any temperature T is given by-

(1)

ΔH^T_f and, ΔS^T_f can be calculated using the following relations:-

(2)

(3)

While the values for ΔH°_f, standard enthalpy of formation for any species at 298.15 K, kJ/mol have been obtained directly from the available literatures, the values for ΔS°_f, standard entropy of formation of any species at 298.15 K, kJ/molK had to be calculated using the following relation [6].

(4)

By convention, ΔH°_f and ΔS°_f values for elements are assumed equal to zero.

Using standard temperature dependent specific heat (C_p) correlations for different species in the reaction system, ΔG^T_f correlations for all species are derived using f equations (1), (2) and (3) [6,7]. These correlations are depicted in Table 1.

The total Gibbs free energy of the system is given by-

(5)

For gas phase reactions,

(6)

Since the standard state is taken as pure ideal gas state at 1 bar, f₁^o =1.

Substituting this we get

(7)

The fugacity coefficient (Φ) of each species in the mixture is calculated using the Soave-Redlich-Kwong (SRK) equation of state [8].

(8)

where Z is also calculated using the SRK equation of state:-

(9)

where,
B_i =b_i P/RT,
B=bP/RT= (Σy_i b_i) P/RT,
A=aαP/ (RT)² =(Σy_i Σy_j aα_ij)P/(RT)²,
aα_ij = (aα_i * aα_j)^½ (1- k_ij),
a_i = 0.42747R² T_c² /P_c
b_i = 0.08664RT_c /P_c,
α_i = [1+m(1-(T/T_c)^½ )]² ,
m= 0.48+1.574ω_i -0.176ω_i²

The following constraints must be satisfied for minimization of total Gibbs free energy:-

• Elemental balance constraint: The number of moles of each element remains conserved throughout the reaction. This can be represented as-
  (10)

  
• Non-negativity constraint: The values of the outlet stream composition (n_i) must always be strictly greater than zero since the calculations involve logarithmic terms.

All the calculations in the preceding stages are done by assuming an outlet stream temperature [7]. Therefore, it is necessary to crosscheck this assumption. This is done by carrying out the following heat balance calculation over the Claus furnace.

(12)

(13)

The heat generated in the reaction is the heat of reaction, which is given by the following relation:-

(14)

Thus, the correct temperature should satisfy the relation given below.

(15)

If the relation given by equation (15) is not satisfied a new value for outlet temperature is assumed and the whole process is repeated till convergence is achieved within acceptable deviation.

Table 1. ΔG^T_f Correlations for species in the reaction mixture

To get the correct heat balance, it is necessary to calculate the heat loss from the walls of the furnace correctly. Only heat loss due to conduction from the shell walls of the Claus furnace needs to be considered in this case [7].

The furnace is made up of three layers, first layer is the thin shell, and other two are insulating firebricks. Heat loss calculation is based on conductivity of the insulating layers, maximum shell wall temperature, wind velocity and ambient temperature (available from plant data).

The heat loss, Q_L is given by the following equation:

(16)

The Gibbs free energy minimization is carried out with the help of Solver function, pre-installed in Microsoft Excel. It applies the Generalized Reduced Gradient (GRG) method to solve non-linear systems for a given set of data [1]. Typical calculation algorithm is shown in Figure 1 and Solver Parameter Dialogue Box in Figure 2.

The input data includes chemical identities of the reaction species, feed composition, furnace inlet temperature and reaction pressure. All other thermodynamic parameters namely, critical temperature, critical pressure, acentric factor, binary interaction parameter, standard heat of formation, standard Gibbs free energy of formation, absolute entropy are supplied by the user on spreadsheet based on literature data [7,9,10].

Besides these, the user has to supply the initial estimates for furnace outlet composition, outlet stream temperature, and compressibility factor, Z. The initial estimate for Z is unity.

Figure 1. Algorithm for Solving the Model

Figure 2. Solver Parameter Dialogue Box

The various equations described in the previous section are entered to calculate the necessary values using above-mentioned input data. A model spreadsheet is shown in Tables 5 (A), 5 (B) and 5 (C). To operate the Solver, the user needs to specify the objective function, the decision variables, and the constraint functions in the 'Solver Parameter' dialogue box.

After these have been defined and the Solver search conditions and parameters have been adjusted, the Solver is run by clicking on the 'Solve' button. When an optimal solution is achieved, Solver displays a message 'Solver found a solution, all constraints and optimality conditions are satisfied'. This indicates that the Solver has found a locally optimal solution. The Solver should be run from several sets of initial estimates to ensure that a globally optimal solution has been found with the minimum possible value of nG. Since Solver is a completely mathematical optimization tool, widely different sets of initial estimates may lead to different sets of solutions. Therefore, to ensure that the globally optimal solutions is obtained, it is advisable to provide initial estimates close to the optimal values and look for the solution with the lowest value of nG.

On this model spreadsheet, the total Gibbs free energy of the system (nG), given by equation (7) occupies the objective function target cell (T90) which needs to be minimized.

The outlet stream composition (moles) and the compressibility factor, Z represented by cell numbers B14 toR14 and E80 respectively, occupy the decision variable change cells.

The equation '$H$7:$M$7=$T$80:$T$85' represents elemental balance constraint on the spreadsheet.

The cells B14 to R14 must be maintained strictly greater than zero to satisfy the non-negativity constraint. This constraint tends to get violated in the GRG method used by the Solver. Thus, to keep the mole numbers strictly greater than zero so that the Solver runs without being interrupted by error messages, MAX function present in the Excel spreadsheet can be used to replace the nonnegativity constraint. For this, the cells B13 to R13 (initial estimates) should be labelled as change cells instead of B14 to R14. Simultaneously, MAX function should be used in the cells B14 to R14 so that they represent the largest value between the corresponding change cells and a small positive quantity. For example, cell B14 should be formulated as 'MAX (0.0001,B13)'.

The value of the function f(Z) given by equation (9) must be equal to zero. It is represented by cell number F80.

The results provided by the Solver for the outlet stream composition and the outlet stream temperature are compared with the plant (referred case study) furnace outlet composition at different SRU plants using corresponding feed compositions and inlet temperature. The comparison is depicted in the following tables (Tables 2- 4).

Solver calculated outlet composition with respect to key components like H₂S, O₂, H₂O, SO₂ and S_x for Plant 1 found to match reasonably well with actual plant outlet composition.

The calculated outlet temperature is for plant 1 found to be 1544 °K, whereas, the outlet temperature according to plant data is 1541 °K.

Estimated temperature for Plant 1 accounts for calculated heat loss of 365293.38 kJ/hr, which is about 2.97% of the heat of reaction.

Solver calculated outlet composition with respect to key components like H₂S, O₂, H₂O, SO₂ and S_x for Plant 2 found to match reasonably well with actual plant outlet composition.

Similarly, the calculated outlet temperature for Plant 2 is found to be 1551 °K, whereas the outlet temperature according to plant data is 1563.5 °K.

Solver calculated outlet composition with respect to key components like H₂S, O₂, H₂O,SO₂ and S_x for Plant 3 found to match reasonably well with actual plant outlet composition.

The calculated outlet temperature for Plant 3 is found to be 1540.5 °K, whereas the outlet temperature according to design data is 1546 °K.

Table 2. Comparison of Solver calculated and Plant Outlet Composition (Plant 1)

Table 3. Comparison of Solver calculated and Plant Outlet Composition (Plant 2)

Table 4. Comparison of Solver calculated and Plant Outlet Composition (Plant 3)

From the results, it is concluded that the Solver Feature in Microsoft Excel can be effectively used to predict the Claus Furnace Outlet Composition based on Gibbs Free Energy Minimization technique. The results give reasonably good prediction with respect to critical components like H₂S, SO₂, S₂, H₂O, CO₂ etc for further calculation of Claus catalytic reaction steps. This in turn provides a good process-engineering basis for the Claus Reaction Furnace and Sulphur Recovery Technology covering the steam generation potential. The deviations in the predictions for non-key components (with low concentrations) like COS, CS₂, NH₃ and CO may be due to convergence criteria of solver correlation used to estimate the Gibbs free energy of these components. However, all above three results indicate that this tool can be used for quick prediction/evaluation of Claus Furnace design/operation with respect to key components like H₂S, O₂, H₂OSO₂, S_x and temperature.

Table 5(a). Model Spreadsheet A1 to T49

Table 5(b). Model Spreadsheet A50 to T105

Table 5( c). Model Spreadsheet A106 to T137

Authors are thankful to Dr. Amalendu Datta for providing his advice during preparation of this article. Authors are also thankful to TPKTI management for giving permission to publish this article.