*Article* **Minimum Entropy Generation Rate and Maximum Yield Optimization of Sulfuric Acid Decomposition Process Using NSGA-II**

## **Ming Sun 1, Shaojun Xia 1, Lingen Chen 2,3,\*, Chao Wang <sup>1</sup> and Chenqi Tang <sup>1</sup>**


Received: 24 July 2020; Accepted: 19 September 2020; Published: 23 September 2020

**Abstract:** Based on the theory of finite-time thermodynamics (FTT), the effects of three design parameters, that is, inlet temperature, inlet pressure, and inlet total mole flow rate, of a tubular plug-flow sulfuric acid decomposition reactor on the total entropy generation rate (EGR) and SO2 yield are analyzed firstly. One can find that when the three design parameters are taken as optimization variables, the minimum total EGR and the maximum SO2 yield of the reference reactor restrict each other, i.e., the two different performance objectives cannot achieve the corresponding extremum values at the same time. Then, the second-generation non-dominated solution sequencing genetic algorithm (NSGA-II) is further used to pursue the minimum total EGR and the maximum SO2 yield of the reference reactor by taking the three parameters as optimization design variables. After the multi-objective optimization, the reference reactor can be Pareto improved, and the total EGR can be reduced by 9% and the SO2 yield can be increased by 14% compared to those of the reference reactor. The obtained results could provide certain theoretical guidance for the optimal design of actual sulfuric acid decomposition reactors.

**Keywords:** finite-time thermodynamics; sulfuric acid decomposition; tubular plug-flow reactor; entropy generation rate; SO2 yield; multi-objective optimization

## **1. Introduction**

At present, the Hybrid-Sulphur (H-S) thermochemical cycle and the Sulphur-Iodine (S-I) thermochemical cycle are considered to be the two most promising recycling methods in the preparation of hydrogen from water by thermochemical cycles [1], and the schematic diagram of S-I thermochemical cycle is shown in Figure 1. Both the H-S and the S-I cycles contain the sulfuric acid decomposition process. Therefore, it is important and necessary to improve the performance of the sulfuric acid decomposition process.

The S-I thermochemical cycle consists of three main chemical reactions: (1) the endothermic decomposition of hydrogen iodide in gas phase; (2) the spontaneous absorption of sulfur dioxide in liquid phase; (3) the sulfuric acid decomposition reaction. The corresponding reaction equations are given as follows:

$$\mathrm{H\_2SO\_4} \stackrel{800\mathrm{K}}{\rightarrow} \mathrm{SO\_3} + \mathrm{H\_2O} \tag{l}$$

$$\text{SO}\_3 \stackrel{1100\text{K}}{\rightarrow} \text{SO}\_2 + \frac{1}{2}\text{O}\_2 \tag{\text{II}}$$

**Figure 1.** The schematic diagram of S-I thermochemical cycle.

Reaction type (I) is the spontaneous decomposition of H2SO4 into SO3 and H2O at 400–500 ◦C. Reaction type (II) is the reaction of SO3 over 750 ◦C to produce SO2 and O2 under the action of a catalyst. In this process, a great deal of heat is consumed, which is also the main energy consumption process in the S-I and H-S thermochemical cycles.

In the aspect of thermodynamic analysis and optimization of sulfuric acid decomposition, Van der ham et al. [1] assumed that the reaction mixture satisfies the ideal gas equation of state, established the physical model of sulfuric acid decomposition reaction, and analyzed the minimization of entropy generation rate (EGR) of a sulfuric acid decomposition reactor by using the optimal control theory. Kuchi et al. [2] carried out a numerical simulation of a high-temperature shell and tube heat exchanger and decomposer, investigated the fluid flow, heat transfer, and chemical reaction processes in the decomposer by using the porous media method, and established a two-dimensional axisymmetric tubular plug-flow reactor model. Ponyavin et al. [3] studied the sulfuric acid decomposer process in a high-temperature ceramic heat exchanger and established a three-dimensional calculation model of the reactor. Van der ham et al. [4] further compared two methods to improve the efficiency of sulfuric acid decomposition reactor and proposed two design schemes to improve the efficiency of the reactor. On the basis of Ref. [1], Wang et al. [5,6] optimized the decomposition of sulfuric acid in the tubular plug-flow reactor with the goal of maximum yield [5], further analyzed the influences of the design parameters of the reactor on the SO2 yield and specific EGRs [6], and obtained the optimal parameters corresponding to the minimum specific EGRs.

Many scholars have optimized other types of thermochemical reaction processes by using the theory and method of finite-time thermodynamics (FTT) [7–22]. For example, Wang et al. [23] investigated the isotherm chemical reaction *A*⇔*B*⇔*C* and obtained the best concentration configuration of the reaction. Johannessen and Kjelstrup [24] studied the EGR minimization of sulfur dioxide oxidation process. The second-generation non-dominated solution sequencing genetic algorithm (NSGA-II) has been widely used in multi-objective optimization of various engineering problems [25–30].

On the basis of Refs. [1,5,6], this paper will further analyze the effects of reactant inlet temperature, pressure, and total molar flow rate on total EGR and SO2 yield, and perform the multi-objective optimization of the process by using the NSGA-II algorithm by applying FTT.

#### **2. Modeling of the Sulfuric Acid Decomposition Process**

A reference reactor used in the performance analysis and optimization as well as the kinetics and thermodynamics models will be introduced in this section.

## *2.1. Reference Reactor*

The model of a tubular plug-flow reactor for sulfuric acid decomposition is shown in Figure 2. It is assumed that the temperature (*T*w) of the outer wall of tubular plug-flow reactor does not change with time and its distribution is linear along the axial direction of the reactor. The distribution follows *T*<sup>w</sup> = 975 + 148*z*/*L* (K). The reaction mixture in the reactor is regarded as an ideal gas and only flows along the axial direction of the reactor. The radial concentration gradient and temperature gradient of

the reaction mixture in the reactor are ignored without both radial diffusion and back-mixing. The total molar flow rate and velocity of the reaction mixture at the cross-section of the reactor are as follows:

$$F\_{\text{tot}} = \sum\_{i} F\_{i} \tag{1}$$

$$\mathbf{v} = \frac{F\_{tot}}{A\_c} \frac{R}{P \times 10^5} \tag{2}$$

where *Fi* is the molar flow rate of reaction component *i*, i.e., H2SO4, SO3, H2O, SO2 and O2; *Ac* is the radial cross section area of the reactor, and *R* is the universal gas constant.

**Figure 2.** Schematic of tubular plug-flow reactor.

The data of catalyst selection, reactor structure, and thermodynamic parameters of the reaction mixture are determined according to Ref. [1], as listed in Table 1.


**Table 1.** Parameters of the reference reactor.

#### *2.2. Models of Kinetics and Thermodynamics*

The fluid flow, heat transfer, and chemical reaction of the reaction mixture in a tubular plug-flow reactor follow momentum, energy, and mass conservation equations, respectively, which are given by:

$$\frac{d\mathbf{dP}}{d\mathbf{z}} = -\left[\frac{150\eta}{D\_{\mathbf{P}}^2} \frac{\left(1-\varepsilon\right)^2}{\varepsilon^3} + \frac{1.75\rho\_{\rm in}\mathbf{v}\_{\rm in}}{D\_{\mathbf{P}}} \frac{1-\varepsilon}{\varepsilon^3}\right] \mathbf{v} \tag{3}$$

$$\frac{\mathrm{d}T}{\mathrm{d}z} = \frac{\pi D I\_{\mathrm{q}} + A\_{c} \rho\_{p} \sum\_{j} \left[ r\_{\mathrm{m},j} (-\Delta\_{r} H\_{j}) \right]}{\sum\_{i} \left( F\_{i} \mathbb{C}\_{p,i} \right)} \tag{4}$$

$$\frac{\mathrm{d}F\_{\mathrm{H}\_{2}\mathrm{SO}\_{4}}}{\mathrm{d}z} = -A\_{\mathrm{c}}\rho\_{\mathrm{p}}r\_{\mathrm{m},1} \tag{5}$$

$$\frac{d\mathbf{F}\_{\rm FL\_2O}}{d\mathbf{z}} = A\_{\rm c} \rho\_{\rm P} r\_{\rm m,1} \tag{6}$$

$$\frac{d F\_{\rm SO\_3}}{d z} = A\_{\rm c} \rho\_{\rm P} (r\_{\rm m,1} - r\_{\rm m,2}) \tag{7}$$

$$\frac{d F\_{\rm SO\_2}}{dz} = A\_{\rm c} \rho\_{\rm p} r\_{\rm m,2} \tag{8}$$

$$\frac{\mathrm{d}F\_{\mathrm{O}\_{2}}}{\mathrm{d}z} = \frac{1}{2} A\_{\mathrm{c}} \rho\_{\mathrm{p}} r\_{\mathrm{m},2} \tag{9}$$

where ρin and vin are the density and flow velocity of the reaction mixture on the entrance section, respectively; subscript *j* = 1, 2 represents the reaction types (I) and (II); *rm*,*<sup>j</sup>* is the reaction rate of mass per unit catalyst, and they are *rm*,1 = *r*1/ρ<sup>p</sup> and *rm*,2 = *r*2; *Cp*,*<sup>i</sup>* and Δ*rHj* are the component molar constant-pressure heat capacity and the reaction enthalpy of the reaction type *j*, and their expressions are given in the Appendix A.

The heat transfer from the heat source outside the tube to the reaction mixture inside the tube follows Newtonian heat transfer law:

$$J\_{\mathbf{q}} = \mathcal{U}(T\_{\mathbf{w}} - T) \tag{10}$$

For different reaction conditions and mechanisms, the driving force in the kinetic equation could be written as different mathematical forms, and the corresponding coefficients in the kinetic equation should be determined by experiments and also be different for different choices of the driving force. According to Ref. [1], the condition that the chemical reaction occurred at the vicinity of the equilibrium is assumed to be satisfied, and all components are assumed to have stoichiometric reaction order, so the chemical reaction rates of reaction types (I) and (II) are as follows:

$$r\_1 = k\_1 \left( P\_{\text{H}\_2\text{SO}\_4} - \frac{P\_{\text{H}\_2\text{O}}P\_{\text{SO}\_3}}{K\_1} \right) \tag{11}$$

$$r\_2 = k\_2 \left( P\_{\rm SO\_3} - \frac{P\_{\rm SO\_2} \sqrt{P\_{\rm O\_2}}}{K\_2} \right) \tag{12}$$

where *k*<sup>1</sup> and *k*<sup>2</sup> are the reaction rate constants of reaction types (I) and (II), according to Ref. [1], *<sup>k</sup>*<sup>1</sup> <sup>=</sup> <sup>10</sup><sup>−</sup>3mol(*SO*3)/(Pa·m3·s), *<sup>k</sup>*<sup>2</sup> <sup>=</sup> 4.7×10−<sup>3</sup> exp( <sup>−</sup>99·10<sup>3</sup> *RT* ) mol(*SO*3)/(Pa·kg·s); *P* represents the partial pressure of the corresponding component; *Kj* <sup>=</sup> exp( <sup>Δ</sup>*rG*◦ *T*,*j* <sup>−</sup>*RT* ) is the equilibrium constant of the chemical reaction type *<sup>j</sup>*; <sup>Δ</sup>*rG*◦ *<sup>T</sup>*,*<sup>j</sup>* is the standard Gibbs free enthalpy of the reaction type *j*, and the expression is given in the Appendix A. The driving force in the kinetic Equation (12) is written as *r*<sup>2</sup> = *k*<sup>2</sup> - *P*SO3 − *P*SO2 \$ *P*O2/*K*<sup>2</sup> , and effects of the different forms of the driving force on the optimization results will be considered in another paper in the future.

The SO2 yield of the tubular plug-flow reactor is as follows:

$$
\Delta F\_{\text{SO}\_2} = F\_{\text{SO}\_2, \text{out}} - F\_{\text{SO}\_2, \text{in}} \tag{13}
$$

The local EGR of the tubular plug-flow reactor is as follows:

$$\begin{aligned} \sigma\_{\text{tot}} &= \sigma\_{\text{ht}} + \sigma\_{\text{f}} + \sigma\_{\text{cr}} \\ &= \pi D I\_{\text{f}} \Big( \frac{1}{T} - \frac{1}{T\_{\text{w}}} \Big) + A\_{\text{c}} \text{v} \Big[ -\frac{1}{T} \Big( \frac{\text{d}P}{\text{d}z} \Big) \Big] + A\_{\text{c}} \rho\_{\text{b}} \sum\_{j} r\_{m,j} \Big( -\frac{\Delta\_{r} G\_{j}}{T} \Big) \end{aligned} \tag{14}$$

where subscripts ht, f, and cr represent the local EGRs of heat transfer, fluid flow, and chemical reaction, respectively.

The total EGR is obtained by integrating the local EGR, i.e.,

$$
\Sigma\_{\rm tot} = \int\_0^L \sigma\_{\rm tot} dz \tag{15}
$$

## **3. Parameter Analyses of Sulfuric Acid Decomposition Reactor**

By changing the inlet parameters of the reference reactor, including the inlet temperature *T*in, pressure *P*in and the total molar flow rate *F*tot,in, the total EGR and the SO2 yield of the reference reactor are analyzed, and the influences of the initial inlet conditions on the two performance objectives can be obtained. The variation ranges of the initial inlet parameters are: 750 K ≤ *T*in ≤ 900 K, 4 MPa ≤ *P*in ≤ 9.5 MPa, and 0.0027 mol/s ≤ *F*tot,in ≤ 0.1 mol/s.

Figure 3 shows the effects of the temperature *T*in of the reaction mixture on the total EGR and the SO2 yield. It can be seen that the total EGR decreases nonlinearly with the increase of the temperature *T*in, and the decreasing trend is fast firstly and then slow; when the temperature *T*in increases from 750 ◦C to 900 ◦C, the total EGR decreases from 0.331 W/K to 0.189 W/K, i.e., decreases by 43%. The main reason is that with the temperature *T*in of the reaction mixture increases, the heat transfer temperature difference between the reaction mixture and the external heat source decreases, which reduces the local EGR of heat transfer and the total EGR. The SO2 yield increases very slowly with the increase of the temperature *T*in, and when the temperature *T*in increases from 750 ◦C to 900 ◦C, the SO2 yield increases by only 0.4%. It can be seen that the total EGR can be reduced by increasing the temperature *T*in of the reaction mixture, i.e., the irreversibility of the sulfuric acid decomposition process could be reduced by increasing the *T*in of the reaction mixture. However, it is not significant to increase the SO2 yield by increasing the temperature *T*in of the reaction mixture.

Figure 4 shows the effects of the pressure *P*in of the reaction mixture on the total EGR and the SO2 yield. It can be seen that the curve of the total EGR is concave and parabolic-like with the increase of the pressure *P*in, and the minimum value is 0.224 W/K when the pressure *P*in is about 0.85 MPa. The SO2 yield decreases linearly with the increase of the pressure *P*in. When the pressure *P*in increases from 0.4 MPa to 1 MPa, the SO2 yield decreases from 0.0118 mol/s to 0.0105 mol/s, i.e., decreases by 11.02%.

Figure 5 shows the effects of the molar flow rate *F*tot,in of the reaction mixture on the total EGR and the SO2 yield. It can be seen that the total EGR and the SO2 yield increase with the increase of the molar flow rate *F*tot,in, and the minimum total EGR and the maximum SO2 yield are mutually restricted. When the molar flow rate *F*tot,in increases from 0.027 mol/s to 0.10 mol/s, the total EGR and the SO2 yield increases by 4.8 times and 1.8 times, respectively.

**Figure 3.** The effects of *T*in on the total EGR and the SO2 yield.

**Figure 4.** Effects of *P*in on the total EGR and the SO2 yield.

**Figure 5.** Effects of *F*tot,in on the total EGR and the SO2 yield.

## **4. Multi-Objective Optimization and Result Analyses**

From the analyses in Section 3, when the three inlet parameters are chosen as optimization variables, and the minimum total EGR and the maximum SO2 yield are taken as optimization objectives, respectively, there is no optimal solution to achieve the extremum values of the total EGR and SO2 yield at the same time. Therefore, how to select the appropriate initial inlet conditions to achieve the relative optimal total EGR and SO2 yield is very important. The NSGA-II algorithm is one of the excellent algorithms to solve multi-objective optimization problems, and can give a series of non-inferior solutions (solutions that cannot be optimized for arbitrary objectives without making other objectives worse) of multi-objective problems. The corresponding improvement process is called Pareto improvement, the corresponding set of non-inferior solutions is called the Pareto-optimal solution set, and the corresponding objective function solution is called the Pareto-optimal front.

Figure 6 shows the flow chart of the NSGA-II algorithm. In this section, all of the *T*in, *P*in and *F*tot,in are taken as the optimization variables to minimize the total EGR and maximize the SO2 yield. The optimization intervals of the variables are consistent with the previous single-variable analysis.

**Figure 6.** Basic flow chart of NSGA-II algorithm.

Figure 7 is Pareto optimal frontier of a reference reactor based on the objective of minimizing total EGR and maximizing SO2 yield, where points A and B represent the solution of the maximum SO2 yield and the minimum total EGR, respectively. At point A, the weighting coefficient of SO2 yield in multi-objective optimization is 1, and the weighting coefficient of total EGR is 0, it is also the solution of maximizing the SO2 yield. Similarly, point B is the solution of minimizing the total EGR. From Figure 7, it can be seen that the minimum total EGR and the maximum SO2 yield are mutually constrained, and they cannot achieve the extremum values at the same time. Only the relative optimal solutions of the two objectives under different weighting coefficients can be found, that is, the non-inferior solution. One can select the appropriate optimal solution from the Pareto-optimal solution set according to different needs to meet the different demands of decision-making purposes. Commonly used multi-objective decision-making methods are Shannon, LINMAP, and TOPSIS, but in the actual decision-making process, decision-making is usually based on actual engineering experience and personal preferences of decision-makers, there is no universal way to make decisions.

**Figure 7.** Pareto optimal frontiers of reference reactor.

In this paper, in order to facilitate the comparison with the reference reactor, a suitable multi-objective decision point (point C) is selected for comparison. Because the solution of the minimum specific EGRs is the solution of the total EGR and the yield under a certain ratio, the decision point of the minimum specific EGR must be on the Pareto-optimal front, which can be used as an important basis to verify the accuracy of the NSGA-II algorithm results.

Figure 8 is the bar chart of the target value of the reference reactor under optimization and non-optimization. Table 2 lists the results of each optimization target condition. It can be seen that compared with the reference reactor, the SO2 yield of the reactor with the maximum yield increases by 118%, but the total EGR increases by 222%; the total EGR of the minimum EGR reactor decreases by 40%, and the corresponding SO2 yield also decreased by 22%; the total EGR and the SO2 yield of the reactor with the minimum specific EGR decrease by 38% and 16%, respectively. From Figure 7, it can be easily concluded that the reference reactor is not located at the Pareto optimal frontier, so the reference reactor can be optimized by Pareto improvement. A non-inferior solution (point C) is obtained by the multi-objective optimization method, in which the total EGR of the reactor decreases by 9% and the SO2 yield of the reactor increases by 14% compared to the reference reactor. Also, from Figure 7, it can be seen that a series of non-inferior solutions located at the upper left of the decision point (point E) of the reactor have good properties of reducing the total EGR and increasing the SO2 yield.

Figures 9–11 show the distribution of the *T*in, *P*in and *F*tot,in in Pareto-optimal fronts, and the black and white dots in the figures represent the total EGR and the SO2 yield, respectively, which exist in pairs. As seen from Figures 9 and 10, the *T*in and *P*in of the reaction mixture in Pareto-optimal fronts are mainly distributed in high-temperature (892–896 K) and high-pressure (9.0–9.2 bar) area, so increasing the *T*in and *P*in of the reaction mixture is an important means for Pareto improvement. Figure 11 shows that the *F*tot,in of the reaction mixture in Pareto-optimal fronts distributes uniformly in

its optimal range, which indicates that adjusting the *F*tot,in of the reaction mixture in Pareto-optimal fronts is an important means to reconcile the contradiction between the minimum total EGR and the maximum SO2 yield.

**Figure 8.** Comparison of total EGR and the yield of optimized objectives.



**Figure 9.** Distribution of inlet temperature in Pareto-optimal fronts.

**Figure 10.** Distribution of inlet pressure in Pareto-optimal fronts.

**Figure 11.** Distribution of total inlet molar flow rate in Pareto-optimal fronts.

## **5. Conclusions**

In this paper, the effects of reaction mixture inlet parameters on the total EGR and SO2 yield of the tubular plug-flow sulfuric acid decomposition reactor are analyzed, and the multi-objective optimization for the two performance objectives are carried out by using FTT. The results show that:


**Author Contributions:** Conceptualization, M.S., S.X. and L.C.; Funding acquisition, S.X.; Methodology, M.S., S.X. and L.C.; Software, M.S., C.W. and C.T.; Supervision, L.C.; Validation, M.S. and C.T.; Writing—original draft, M.S. and S.X.; Writing—review & editing, L.C. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by National Natural Science Foundation of China (Grant Nos. 51976235) and 51606218) and the Hubei Province Natural Science Foundation of China (Grant No, 2018CFB708).

**Acknowledgments:** The authors wish to thank the Academic Editor and reviewers for their careful, unbiased and constructive suggestions, which led to this revised manuscript.

**Conflicts of Interest:** The authors declare no conflict of interest.

## **Nomenclature**


## **Appendix A**

According to the Refs. [31], the component molar constant-pressure heat capacity, molar enthalpy and molar Gibbs energy can be calculated by the following formula:

$$\mathcal{L}\_{p,l}^{\circ} = A\_l + B\_l \frac{T}{1000} + C\_l \left(\frac{T}{1000}\right)^2 + D\_l \left(\frac{T}{1000}\right)^3 + E\_l \left(\frac{1000}{T}\right)^2\tag{A1}$$

$$H\_{T,l}^{\circ} = A\_l \frac{T}{1000} + \frac{B\_l}{2} \left(\frac{T}{1000}\right)^2 + \frac{C\_l}{3} \left(\frac{T}{1000}\right)^3 + \frac{D\_l}{4} \left(\frac{T}{1000}\right)^4 - E\_l \left(\frac{1000}{T}\right) + F\_l \tag{A2}$$

$$S\_{T,l}^{\circ} = A\_l \ln \left( \frac{T}{1000} \right) + B\_l \frac{T}{1000} + \frac{C\_l}{2} \left( \frac{T}{1000} \right)^2 + \frac{D\_l}{3} \left( \frac{T}{1000} \right)^3 - \frac{E\_l}{2} \left( \frac{1000}{T} \right)^2 + G\_l \tag{A3}$$

$$
\Delta\_{\rm r}H\_T^\circ = \sum\_l \upsilon\_l H\_{T,l}^\circ \tag{A4}
$$

$$
\Delta\_r S\_T^\circ = \sum\_l \upsilon\_l S\_{T,l}^\circ \tag{A5}
$$

$$
\Delta\_{\varGamma} \mathcal{G}\_{\varGamma}^{\diamond} = \Delta\_{\varGamma} H\_{\varGamma}^{\diamond} - T \Delta\_{\varGamma} \mathcal{S}\_{\varGamma}^{\diamond} \tag{A6}
$$

where *Ai*∼*Gi* are the thermodynamic coefficients of the formula, which are listed in Table A1; *vi* is the stoichiometric number of reaction component *i*.

**Table A1.** Thermodynamic coefficients.


#### **References**


© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).
