A Systematic Review of the Literature on Steady-State Reactive Distillation Modeling and Simulation: Challenges and Opportunities

Abstract

Reactive distillation has already been extensively studied, many applications have been evaluated by simulation and laboratory experiments, and many modifications in the configuration of the column are being proposed, such as the reactive dividing wall column. The steady state plays a key role because it is required to operate a process steadily at the optimal point, and it is necessary to use control strategies to keep the operation at this point. Thus, this chapter presents a systematic literature review, about the modeling and simulation of reactive distillation in steady state presenting some bibliometric results. It also presented some content related to the variations in the configurations of a reactive distillation column, mathematical modeling of the process, and computational simulation. Lately, simulation studies reported in the literature about this subject have mainly used sequential modular simulators like Aspen Plus. But, the solution process of the model’s equations still can be improved. Some opportunities are pointed out, like the use of heuristic algorithms for the solutions of the mathematical equations, better initial guess algorithm for the unknowns, and the use of a more generic scheme for representation of the mathematical model of a reactive distillation stage.

1. Introduction

Many traditional processes related to the chemical industries comprise different stages, in which reaction and separation are among the most important ones [1]. The reaction occurs in different reactor formats (tank or tubular) and operation types (batch or continuous), and the reaction can need some catalyst (homogeneous or heterogeneous). Regarding the separation step, distillation is one of the most widely used unit operations to separate mixtures in chemical industry applications [2, 3]. In recent years, there has been permanently increasing interest in the development of hybrid processes combining reaction and separation, among which can be highlighted reactive distillation processes (configuration in which reaction and separation happen in the same column) [1, 4, 5, 6, 7, 8, 9], known as a reactive distillation column (RDC), and the phenomenon is, thereby, referred to as reactive distillation (RD) process [10, 11, 12, 13, 14, 15, 16]. In this way, several new processing methods have been developed and commercialized [17]. When reactive distillation is compared to conventional processes, one can see that reactive distillation results in a simpler flowchart with fewer recycle streams and a lower number of separation units [18].

Distillation is considered as a mature technology and has been widely used in many processes over the years; however, some improvements are being made in order to reduce energy consumption, since it is well known for its high energy requirements, and improve thermodynamic efficiency [19, 20, 21]. Furthermore, world energy consumption is increasing at a steady rate for many reasons like population growth, industrialization and transportation [22]. Process intensification is a promising tool to increase the sustainability of chemical processes by using novel apparatuses and techniques [3, 23].

Reactive separation conveniently combines the production and removal of one or more products, enabling enhancement of conversion, simplicity, selectivity and yield, which can present various technical, environmental and economic benefits compared to reaction and separation taking place separately [4, 8, 22, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35].

Another advantage of some reactive distillation applications is that the reaction temperature is easy to control in a reactive distillation column. The temperature of each stage is always the bubble point, and it depends on the reactive mixture and system pressure. Therefore, the reaction temperature can be controlled by adjusting the system pressure [36]. The heat of an exothermic reaction is used for the vaporization of liquids, reducing the reboiler duties. The higher temperature value in the reaction section is limited by the boiling point of the mixture of the components involved in the process; thus, the formation of hot spots is severely reduced [10, 37, 38].

The thermodynamic limitations, for example, the presence of azeotropes, make conventional separation systems complex and expensive but can be overcome using reactive distillation. Of course, the reaction temperature must be suitable to the column conditions [24, 39], and the reactants and products relative volatilities must be such that products are removed and reactants retained inside the column [40].

Reversible reactions can be limited to equilibrium, which is difficult to overcome because they have reached a point in which higher product purity is not possible if the reaction conditions are kept [16]. Reactive distillation can be used to overcome this problem by continuous removal of products from the reaction zone [22, 37]. If a liquid-phase reaction must be carried out with a large excess of one reactant, reactive distillation is potentially attractive because it can be carried out closer to stoichiometric feed conditions, avoiding recycling costs [41]. Side reactions can be dampened or even avoided by the constant separation of the products. Also, the formation of azeotropes can be prevented [10, 42].

Reactive distillation has been employed in industry for many decades [24]. However, the potential of reactive distillation has not yet been fully tapped, and there is still ongoing research to improve [43]. Reactive distillation can be used with a wide variety of reactions; among suitable reactive distillation processes, already studied, are acetylation, alkylation, amination, dehydration, etherifications, nitrations, esterifications, hydrolysis, isomerization, transesterifications, polycondensations, and halogenations [1, 7, 44], offering an attractive alternative for the production of many important industrial chemicals [45]. With respect to the exploration of the potential of the reactive distillation process, the steady-state simulation of a reactive distillation column is an important step in the implementation of the technology and the discovering of its potential advantages, and it is essential for design, investigation of new applications, control and optimization [46, 47].

A literature review is fundamental for knowing the cutting edge of a scientific subject and also can be used to find some interesting questions that have not yet been answered, but there is so much scientific content that is important to systematize the review process, because in a review not systematize the research tends to cite and comment the publications that he/she like or known most, and the systematic literature review can help to avoid this making the review process reproducible. Also, a systematic literature review can provide some bibliometric data that can be used to see where (countries, institutions, etc.) the subject is studied most [48, 49]. Among the many systematic literature review methods proposed in the literature, like PRISMA [50], ProKnow-C [51], SIMILAR [52] and NIRP [48], among others, there are some that help the researcher to obtain a rank of the importance of the papers in the final portfolio, among which is Methodi Ordinatio [53, 54]. Knowing the importance of many papers is important to deciding which must receive more and less attention.

In this work, a systematic review of the literature (applying the Methodi Ordinatio) was carried out about the modeling and simulation of steady-state (SS) reactive distillation column aiming to present a collection of the most used methods and tools used in this research area and find some future studies that should be interesting to consider. The steps of the Methodi Ordination are presented in Section 2. Also, in this section, the choices made in the application of the method are presented, as well as some data collected in this application. After getting the list of articles, some aspects related to the subject were collected and presented in Sections 3, 4, and 5. Some challenges and opportunities are presented in six, and the conclusions about the review are outlined in Section 7.

2. Methodology

This section presents the Methodi Ordination, the data obtained in each step of this method, and some bibliometric results.

2.1 Methodi Ordinatio

In this research, a systematic review of the literature was carried out, aiming at publications involving the steady-state modeling and simulation of reactive distillation. The Methodi Ordinatio, proposed by [53, 54], was used in the construction of the literature review to determine which is the most relevant research in the field.

Literature reviews are important in order to organize papers and to find possible lacks and opportunities in the research field. The Methodi Ordinatio is a systematic literature review method based on an index (InOrdinatio) to generate a ranking grouping of the most important papers. The nine steps of the method are presented below:

1. Research topic definition: before you start research, you must define what you want to look for;

2. Preliminary search: some tests can be carried out in order to clarify what you are searching;

3. Definition of keywords, databases and time delimitation: the previous step can help you to define the combination of keywords. Beyond that, you should choose the databases more relevant to the research topic and, in some cases, it is interesting to limit the range of time for the review;

4. Definitive search in the database: Finally, you have the necessary data to do the research in the databases. In this step, and in some of the upcoming ones, it is primordial to use a reference manager software (like EndNote, Mendeley, JabRef, Zotero, etc.);

5. Filtering procedure: mainly if the search is carried out in more than one database, it is normal to have some duplicate papers, and these must be resolved. Some results present no data, which results in the deletion. Also, it is accepted only for research papers. Conference papers, books, book chapters, and review papers are not accepted, even if it is a good research item, because the impact factor indicator is not available for this kind of research item, which makes it impossible to calculate the index (InOrdinatio) which is used to generate the rank of importance of the research items;

6. Searching for papers data: some data are necessary for each paper is necessary for the ordinatio index calculation that is carried out in the next step;

7. InOrdinatio: the order of the article’s importance is determined by applying the data collected in the last step in Eq. (1)

   $$\text{InOrdinatio}=\left(\text{IF}/1000\right)+\alpha \left[10-\left(\text{ResearchYear}-\text{PublishYear}\right)\right]+\mathrm{Ci}$$E1

   where IF is the journal impact factor in which the article was published, $\alpha$ is a value (from 1 to 10) related to the subject importance and newness, ResearchYear is the year of the data collection (2023 for this review), PublishYear is the year of the article publication, and Ci is the number of times the paper has already been cited.

8. Finding texts in full format: it is interesting to perform this action at this point because at this point the number of papers is much smaller;

9. Systematic reading and analysis of articles: some interesting results can be found by analyzing the bibliometric results and mainly by effectively reading the papers.

After all these steps, some lacks and possible opportunities can be found.

2.2 Review steps data

Each step of a systematic review of literature can result in some data that is not necessarily a result, but it is important to understand the procedure. Below it is presented some of them:

1. Research topic definition: the chosen subject was modeling and simulation of reactive distillation in steady state;

2. Preliminary search: this step was fundamental to build the keyword combination;

3. Definition of keywords, databases and time delimitation: Due to the subject, it was not set time delimitation; the databases chosen are ScienceDirect, Scopus and Web of Science and the keywords combination is:

   ("reactive distillation" OR "RD") AND "steady state" AND ("simulation" OR "modelling" OR "modeling")

4. Definitive search in the database: this step resulted in a total of 547 research items:

   • ScienceDirect: 135;

   • Scopus: 293;

   • Web of Science: 119.

5. Filtering procedure: some deletions of research items were made because of:

   • Duplicates: 174;

   • No research item Data: 8;

   • Type of research item (conference paper, review paper, book chapter, etc.): 68;

   • Out of scope (by reading title, keyword and abstract): 144.

   Resulting in 153 papers. Analyzing this step, one can see that title, keywords and abstract must be well written and properly represent the paper content; otherwise, the paper can be incorrectly excluded from the review.

6. Searching for papers data: All data of the final portfolio were found;

7. InOrdinatio: With all data collected in the previous step it was applied Eq. (1) for the ranking composition. The journal impact factor used was JCR. Due to the recentness of the subject is not so relevant in modeling and simulation, a small value of $\alpha$ must be chosen, we choose $\alpha =1$ because there are publications about the subject over more than three decades, therefore it is not a recent research subject. This is subjective choice is a weakness of the Methodi Ordinatio; therefore, there are already some works proposing some improved methods for ranking papers [48, 49]. In this step, nine papers were removed because they resulted in a negative index. The number of citations is from Google Scholar [55, 56]. So, there are 144 articles in the final portfolio.

8. Finding texts in full format: All articles were downloaded.

9. Systematic reading and analysis of articles: Some important features can be found in this step, and some of them are presented in the next section.

2.3 Bibliometric results

After obtaining the final portfolio, some important information can be found by analyzing the article’s bibliometric data. Figure 1 shows the number of articles per country of the institution in which the first author was affiliated when the work was carried out. In this figure, one can see that just four countries (namely, China, Germany, the United States of America and India) are responsible for 70 articles, almost half of the portfolio. Beyond these countries presented in Figure 1, there are three other countries with three articles each, five countries with three articles each, seven countries with two papers each, and seven countries with one article each. In total, the portfolio contains the first authors affiliated with 30 different countries.

Figure 1.

Main countries of steady-state reactive distillation modeling and simulation articles.

Figure 2 shows the journals with a higher number of articles in the final portfolio. Just five journals (namely, Industrial and Engineering Chemistry Research, Chemical Engineering and Processing: Process Intensification, Computers and Chemical Engineering, Chemical Engineering Science and Chemical Engineering Research and Design) are responsible for 79 articles, more than half of the portfolio. Beyond the 18 Journals presented in Figure 2, there are 19 other journals with one article each, that is, the articles from the final portfolio were published in 37 different journals.

Figure 2.

Main journals of steady-state reactive distillation modeling and simulation articles.

As can be seen in Figure 3, the first publication of the portfolio took place in 1991; after this, there was no publication for 5 years, and then the number of publications has varied over the years, with a slight tendency to rise. Among all these years, 2022 has the greatest number of publications.

Figure 3.

Number of articles on steady-state reactive distillation modeling and simulation over the years.

In the final portfolio, there are 351 different authors, among which 249 authors have only one paper, 66 authors have 2 papers from the final portfolio, 20 authors have 3 papers, 9 authors with 4 paper, 3 authors have 5 papers and only 2 authors have 6 and 7 papers from the final portfolio. Figure 4 shows the main authors, in which Švandová, Z. and Markoš, J., the authors who have the most publications in the final portfolio, have seven papers because they published the same seven works; a similar situation occurs with Taylor, R. and Krishna, R., who published the same five articles.

Figure 4.

Main authors of steady-state reactive distillation modeling and simulation articles.

3. Reactive distillation column configurations

Reactive flash is the simplest reactive distillation configuration where separation is carried out in only one vaporization stage [57]. The reactive distillation column is very similar to the conventional distillation column, as illustrated in Figure 5a, in which the separation process occurs through many stages. If the reaction has a rate high enough to be conducted in pressure and temperature of a distillation column conditions without the need for a catalyst, the reaction takes place in all stages of the column, even in the reboiler and condenser, but normally it is considered that the reaction does not occur in the reboiler and condenser to simplify the solution of the mathematical model.

Figure 5.

Variations on reactive distillation.

In the cases that a catalyst is applied in the reactive distillation, that is, a catalyst column, just some stages are reactive, so it is possible to divide the column into three zones like presented in Figure 5b, namely the rectifying zone, reactive zone and stripping zone. If the reaction is heterogeneously catalyzed, the reactive zone is formed by the stages packed with catalyst. In the case of a homogeneously catalyzed reaction, the reactive section is defined by the location of the feed stream and catalyst properties; for example, if the catalyst can be considered to be non-volatile, the reactive zone is formed by feed stage and the stages below that [58].

A side reactor column (SRC), like the one presented in Figure 5c, is a distillation column connected to one or more reactors that can be used if the distillation conditions do not match the conditions for good reaction rates, the column operates at favorable conditions for separation, while the reactor operates at reaction kinetics favorable conditions [59].

In some cases, the separation of mixtures taking place in a conventional distillation sequence can be thermally inefficient. In these cases, one can turn to a thermally coupled distillation sequences (TCDS) solution (that has been used for the separation of multicomponent mixtures achieving energy savings) like the Petlyuk column (total thermal coupled solution) [19], illustrated in Figure 5d, in which it is used a prefractionator column to separate, for example, a mixture of A, B and C, with A the lightest and C the heaviest components in the mixture. An even better solution that can save energy and capital costs is the use of the prefractionator and the main column in a single shell, as shown in Figure 5e. This solution is called a dividing wall column (DWC) because it is a single shell divided by a wall in the middle of some stages of a conventional distillation column. Dividing wall columns enables the separation of four pure fractions. If the dividing wall column is used for a reactive distillation process, it is called a reactive dividing wall column (RDWC).

Reactive dividing wall column (RDWC) is a process in which reactive distillation and dividing wall column are integrated, being that it simultaneously has the advantages of both processes [60]. For heterogeneous catalytic reactions, the reaction zone is normally present on the feed side of the column where the reaction rate is higher because it is where the reactant concentrations are higher, but there are some propositions in the literature for more than one reactive zone [22]. By combining an enzymatic catalyst and a reactive dividing wall column, we have an enzymatic catalyzed reactive dividing wall column (eRDWC).

Heat pump (HP) distillation is a technology that uses the latent heat from the vapor top product as the heat source in the reboiler [61]. There are two mainly conventional ways to use the heat pump in a distillation column (or reactive distillation column): the vapor compression column (VC) and the vapor recompression column (VRC). As shown in Figure 5f and g, in the VC case, the working fluid is evaporated in the condenser, compressed and heated over the reboiler temperature, condensed in the reboiler and cooled down, through a valve, below the condenser temperature, and in the VRC case the working fluid is the vapor leaving the top of the column, and after cooled down through the valve part of the working fluid is refluxed to the column. If a vapor compression or recompression is applied to a reactive distillation column, we have a heat pump for reactive distillation (HPRD).

To improve energy efficiency, the heat integration concept can be applied to a distillation column. Heat-integrated reactive distillation is beneficial for energy saving and economy compared to conventional reactive distillation [62]. The idea of heat integration is the use of hot process streams to raise the temperature of cold process streams. There are several heat integration techniques, one of which is called heat integrated distillation column (HIDiC) technique, Figure 5h. This technique is called heat-integrated reactive distillation column (HIRDiC) if applied to a reactive distillation process [63, 64].

Many other configurations can be found in the literature that are applied in the distillation and reactive distillation, like pressure swing distillation (PSD) or pressure swing reactive distillation (PSRD), which uses two columns operating at two different pressures applied to process separations of pressure-sensitive azeotropes [65].

Using circulating reflux (pump around), normally from the lower stages of a column to stages, makes a stream in a higher temperature to be recirculated in the column, increasing the driving force for heat transfer, which is normal, for example, in a crude oil distillation column that is a multi-product process with many side stream.

Reactive-extractive dividing wall column (REDWC) is a combination of reactive distillation, extractive distillation and dividing wall column, resulting in a highly integrated process that demands less physical space and has a strong interrelation among variables [66].

There are also some combinations that work very well for a specific process that normally need three or more columns, which makes it impractical to be listed here.

Although there are so many configurations in the reactive distillation columns, among which some were presented in this section, the modeling is very similar because it can be made stage-by-stage; therefore, the difference is due to the connections among the stages. Thus, the modeling presented in the next section is focused in the model of one generic stage.

4. Modeling

Although the steady state is not the reality in a chemical process, it is very important to understand the process. The steady-state can show us the limitations of the process, that is, the maximum the process can yield with the current configuration. The steady-state can also be used for optimization and consequently can be important in control because it is desired that the process operates all the time at the optimal point, and the use of control strategies to keep the operation at this point. Because of all these reasons, the mathematical model of a reactive distillation column operating in a steady state is very important to evaluate some features of the process.

A reliable mathematical model is fundamental for understanding and analyzing the behavior of a reactive distillation process [67]. The development of a phenomenological mathematical model of a process is based on applying conservation laws (mass, energy and momentum) and constitutive relations (like chemical equilibrium or chemical reaction kinetics, phase equilibrium, etc.). A mathematical model must represent (approximately) the real behavior of some relevant properties of a process and can be very important to analyze an existing plant or to evaluate the technical viability of a new chemical process plant [68]. Most of the reactive distillation models are originally adapted from conventional distillation calculations based on equilibrium stage model [69]. The reactive distillation column is modeled stage-by-stage with adequate thermodynamic behavior and chemical reactions [14].

For the application of conservation law, for a distillation column is applied in each stage, it is fundamental to know the mass and energy flows. For example, Figure 6 shows a schematic of the flows of a generic stage in a reactive distillation column. This generic stage can be used to represent any stage in the column, even the condenser (stage 1) and reboiler (stage n, last stage). For this reason, it is considered in all stages the liquid side flow from stage j$\left({\mathit{U}}_j\right)$, vapor side flow from stage j$\left(W_j\right)$, heat load from stage j$\left(Q_j\right)$ and flow rate of the feed stream to stage j$\left(F_j\right)$. Furthermore, each stage has an liquid inlet from stage right above $\left(L_{j-1}\right)$, a vapor inlet from the stage right below $\left(V_{j+1}\right)$ and liquid and vapor outlets from stage j (respectively, $L_j$ and $V_j$). If it is considered a homogeneous stage, the liquid and vapor outlets from the stages have the same composition of the stage content, that is, the mole fraction of component i in stage j of liquid da vapor phases ($x_{ij}$ and $y_{ij}$). The feed stream can be in liquid or vapor phase, so $z_{ij}$ is the mole fraction of component i in the feed stream of stage j.

Figure 6.

Schematic of a generic stage in a reactive distillation column.

The equilibrium model of a stage in a reactive distillation column is represented by MESH equations [67, 70, 71, 72, 73, 74]:

• Material balance, global and per component, are represented, respectively, by Eqs. (2) and (3):

$$F_j+R_j+V_{j+1}+L_{j-1}-\left(L_j+{\mathit{U}}_j\right)-\left(V_j+W_j\right)=0$$E2

$$z_{ij}{F}_j+r_{ij}+y_{ij+1}{V}_{j+1}+x_{ij-1}{L}_{j-1}-x_{ij}\left(L_j+{\mathit{U}}_j\right)-y_{ij}\left(V_j+W_j\right)=0$$E3

where $r_{ij}$ is the generation rate of component i on stage j and $R_j$ is the total rate of generation on stage j. The relation between these variables is given by the sum of the generation of all m components in the process, presented by Eq. (4).

$$R_j=\sum _{i=1}^m{r}_{ij}$$E4

• Equilibrium equation. Normally, the liquid-vapor equilibrium, Eq. (5), is considered:

$$y_{ij}=K_{ij}{x}_{ij}$$E5

where $K_{ij}$ is liquid-vapor equilibrium (EQ) relation for the component i on stage j, that for a gamma-phi formulation is given by Eq. (6):

$$K_{ij}=\frac{{\gamma }_{ij}{\varphi }_i^{\mathrm{sat}}{P}_i^{\mathrm{sat}}}{{\widehat{\varphi }}_{ij}P}exp\left({\int }_{P_i^{\mathrm{sat}}}^P\frac{V_i^L}{RT}\mathit{dP}\right)$$E6

where ${\gamma }_{ij}$ is the activity coefficient of component i on stage j, ${\widehat{\varphi }}_{ij}$ is the fugacity coefficient of component i on stage j, ${\varphi }_i^{\mathrm{sat}}$ is the fugacity coefficient of pure component i in saturation state, P is the pressure, $P^{\mathrm{sat}}$ is the saturation pressure, R is the ideal gas constant, T is the temperature, and $V_i^L$ is the component i volume in liquid phase.

• Summation equations. The sum of mole fractions, in the liquid and vapor phases, must be equal to one in all n stages, as presented in Eqs. (7) and (8):

$$\sum _{i=1}^m{x}_{ij}=1$$E7

$$\sum _{i=1}^m{y}_{ij}=1$$E8

• Heat (enthalpy) balance, that is, energy balance presented in Eq. (9):

$$h_j^F{F}_j+H_{j+1}{V}_{j+1}+h_{j-1}{L}_{j-1}-h_j\left(L_j+{\mathit{U}}_j\right)-H_j\left(V_j+W_j\right)-Q_j=0$$E9

where $h_j$ is the enthalpy of liquid stream flow from stage j, $h_j^F$ is the enthalpy of feed stream flow to stage j and $H_j$ is the enthalpy of vapor stream flow from stage j. For this energy balance, the enthalpies must be calculated using the heat of formation as a reference state, so it is not necessary to use an additional term of heat of reaction [1, 75].

Normally, the degrees of freedom in an algebraic system of equations representing the behavior of physical properties in a distillation column are greater than zero. To work around this problem, it is possible to specify some unknown values, but it is habitual to use ratios like those presented in Eqs. (10) and (11):

$${\mathit{U}}_j={\eta }_j^L{L}_j$$E10

$$W_j={\eta }_j^V{V}_j$$E11

where ${\eta }_j^L$ is the ratio between the liquid side stream and the liquid stream outputting stage j and inputting the near stage $j+1$ and ${\eta }_j^V$ is the ratio between the vapor side stream and the vapor stream outputting stage j and inputting the near stage $j-1$.

To evaluate the generation rate of each component in each stage, it is necessary to evaluate the extent of reaction of each reaction in each stage, as shown in Eq. (12):

$$r_{ij}={\delta }_j\sum _{k=1}^p{\nu }_{ki}{\xi }_{kj}$$E12

where ${\xi }_{kj}$ is the extent of reaction k in stage j and ${\nu }_{ki}$ is the stoichiometric coefficient of component i in reaction k.

The calculation of the extent of the reaction can be considered the reaction rate or chemical equilibrium. If the reaction rate is considered, an expression like the one presented in Eq. (13) is used:

$${\xi }_{kj}=V_{hj}{k}_k\left(\prod _{i=1}^m{a}_{ij}^{\left(|{\nu }_{ki}|-{\nu }_{ki}\right)/2}-\frac{1}{K_k}\prod _{i=1}^m{a}_{ij}^{\left(|{\nu }_{ki}|+{\nu }_{ki}\right)/2}\right)$$E13

where $V_{hj}$ is the liquid molar holdup in stage j, $k_k$ is the kinetic constant of reaction k, $K_k$ is the chemical equilibrium constant of reaction k and $a_{ij}$ activity of component i on stage j given by Eq. (14):

$$a_{ij}=x_{ij}{\gamma }_{ij}$$E14

If the reaction is limited by equilibrium, Eq. (15) is used to obtain the values of the mole fraction to reach the chemical equilibrium:

$$K_k=\prod _{i=1}^m{\left(x_{ij}{\gamma }_{ij}\right)}^{{\nu }_{ki}}exp\left(\frac{P_j-P^{°}}{R{T}_j}\sum _{i=1}^m{\nu }_{ki}{V}_i^L\right)$$E15

The chemical equilibrium is reached when the total Gibbs free energy (G) is at its minimum value; in this way, the equilibrium constant is provided by Eq. (16):

$$ln\left(K\right)=-\frac{\mathrm{\Delta }{G}^{°}\left(T\right)}{RT}$$E16

One can see that, for equilibrium limited reaction, the extent of reaction is not present in the model equation, but it appears for applying the definition of the extent of reaction in the calculation of the total number of mols of one component based on the total inlet of this component, as shown in Eq. (17):

$$n_{ij}=n_{ij}^0+\sum _{k=1}^p{\nu }_{ki}{\xi }_{kj}$$E17

where $n_{ij}$ is the quantity of component i in stage j after the reaction reaches the equilibrium based on, $n_{ij}^0$, the total amount of the inlet quantity in the liquid phase (from feed stream or liquid output from stage $j-1$) of component i in stage j. The mole fraction calculated by Eq. (18) must be replaced in Eq. (15)

$$x_{ij}=\frac{n_{ij}}{{\sum }_{i=1}^m{n}_{ij}}$$E18

The Murphree tray efficiency, illustrated in Figure 7 for the McCabe-Thiele method for binary distillation, can be used to better represent the behavior of a mole fraction profile in a reactive distillation column. In this case, the reboiler and condenser are considered in equilibrium, but the tray is considered to have an efficiency compared to the ideal case (vapor-liquid equilibrium). This efficiency can be formulated based on the vapor phase, Eq. (19), or based on the liquid phase, Eq. (19):

Figure 7.

Schematic representation of Murphree tray efficiency in the MaCabe-Thiele diagram.

$$E_{MV}=\frac{y_{ij}-y_{ij+1}}{y_{ij}^{\ast }-y_{ij+1}}$$E19

$$E_{ML}=\frac{x_{ij}-x_{ij-1}}{x_{ij}^{\ast }-x_{ij-1}}$$E20

where $E_{MV}$ is the Murphree tray efficiency for the vapor phase, $E_{ML}$ is the Murphree tray efficiency for the liquid phase, $x_{ij}$ is the liquid mole fraction of component i that actually leaves stage j, $y_{ij}$ is the vapor mole fraction of component i that actually leaves stage j, $x_{ij-1}$ is the liquid mole fraction of component i that actually leave stage $j-1$ (stage right above j) and goes into stage j, $y_{ij+1}$ is the vapor mole fraction of component i that actually leave stage $j+1$ (stage right bellow j) and goes into stage j, $x_{ij}^{\ast }$ is the liquid mole fraction of component i and stage j in equilibrium with the actual vapor composition of stage j$\left(y_{1j},y_{2j},\dots ,y_{mj}\right)$ and $y_{ij}^{\ast }$ is the vapor mole fraction of component i in stage j in equilibrium with the actual liquid composition of stage j$\left(x_{1j},x_{2j},\dots ,x_{mj}\right)$. The relations for $y_{ij}^{\ast }$ and $x_{ij}^{\ast }$ to the equilibrium are presented, respectively, in Eqs. (21) and (22):

$$y_{ij}^{\ast }=K_{ij}\left(x_{1j},x_{2j},\dots ,x_{mj},y_{1j}^{\ast },y_{2j}^{\ast },\dots ,y_{mj}^{\ast },T_j,P_j\right)x_{ij}$$E21

$$y_{ij}=K_{ij}\left(x_{1j}^{\ast },x_{2j}^{\ast },\dots ,x_{mj}^{\ast },y_{1j},y_{2j},\dots ,y_{mj},T_j,P_j\right)x_{ij}^{\ast }$$E22

where $P_j$ is the pressure in stage j.

Another efficiency definition, normally used for avoiding integer optimization (in this case, the stage efficiency can assume continuous variables between 0 and 1), is the bypass efficiency illustrated in Figure 8. This concept assumes that only a fraction (proportional to the bypass efficiency of stage j, ${\epsilon }_j$) of vapor and liquid inlet streams ($L_{j-1}$ and $V_{j+1}$) flow into stage j and reach the equilibrium ($L_j^{\ast }$ and $V_j^{\ast }$), the other fraction of inlet stream bypass and mix with the equilibrium outlet streams [28].

Figure 8.

Illustration of bypass efficiency.

Sometimes, equilibrium (EQ) stage simulations are denominated to be rigorous, but this term is not totally correct because, in actual operation, the phases in a column stage normally operate a little far from the equilibrium [76]. A more realistic (physically consistent) approach available in the literature for modeling the relation between the compositions of vapor and liquid phases is the non-equilibrium (NEQ) in which the finite mass transfer rates across the vapor-liquid interface are accounted for, that is, a rate-based model. However, the rate-based model is much more complicated than the equilibrium model and also more difficult to converge [67]. For this model, it is assumed that the total amount of mass transfer resistance is in the thin films in the border between the vapor-liquid phases, and the mass transfer in these two films is due to molecular diffusion. In the bulk fluid phases, the mixing level is so high that the composition is almost none [1]. The description of interface mass transfer can be based on Maxwell-Stefan theory for the calculation of heat and mass transfer [1, 8, 41, 58, 67, 77, 78]. Near the vapor-liquid interface, there are two boundary layers (one for each phase) in which the composition is variable, as shown in Figure 9, at the vapor-liquid interface can be assumed to be phase equilibrium. In the case this model is considered, energy and molar balances must consider each phase separately [79]. The non-equilibrium mathematical models usually provide more results than the equilibrium models. However, the availability of reliable mass transfer correlations would be a prerequisite for the use of a non-equilibrium stage mode [4].

Figure 9.

Schematic representation of a non-equilibrium stage.

There are also models for heterogeneous catalytic distillation that take into account simultaneous mass transfer and reaction inside the catalyst particle using the Maxwell-Stefan theory [1, 80].

If it is, necessary or desired, many combinations of rate-based or equilibrium reaction and rate-based mass and heat transfer, phase equilibrium or tray efficiency can be used [69].

The transient model also can be used to study the steady-state solution [4, 24, 69, 76, 81, 82]. In this case, the conservation laws become differential equations, and the solution of this model needs a method for solving the stiffness of differential-algebraic equations (DAE) [1, 32, 76, 82, 83, 84, 85]:

$$\frac{d\left(x_{ij}{V}_{hj}\right)}{\mathit{dt}}=z_{ij}{F}_j+r_{ij}+y_{ij+1}{V}_{j+1}+x_{ij-1}{L}_{j-1}-x_{ij}\left(L_j+{\mathit{U}}_j\right)-y_{ij}\left(V_j+W_j\right)$$E23

$$\frac{{\mathit{dV}}_{hj}}{\mathit{dt}}=F_j+R_j+V_{j+1}+L_{j-1}-\left(L_j+{\mathit{U}}_j\right)-\left(V_j+W_j\right)$$E24

$$\frac{d\left(V_{hj}{u}_j\right)}{\mathit{dt}}=h_j^F{F}_j+H_{j+1}{V}_{j+1}+h_{j-1}{L}_{j-1}-h_j\left(L_j+{\mathit{U}}_j\right)-H_j\left(V_j+W_j\right)-Q_j$$E25

where $u_j$ is the internal energy of liquid phase in stage j.

If the liquid molar holdup can be considered as constant, Eqs. (23) and (25) can be simplified as follows:

$$V_{hj}\frac{{\mathit{dx}}_{ij}}{\mathit{dt}}=z_{ij}{F}_j+r_{ij}+y_{ij+1}{V}_{j+1}+x_{ij-1}{L}_{j-1}-x_{ij}\left(L_j+{\mathit{U}}_j\right)-y_{ij}\left(V_j+W_j\right)$$E26

$$V_{hj}\frac{{\mathit{du}}_j}{\mathit{dt}}=h_j^F{F}_j+H_{j+1}{V}_{j+1}+h_{j-1}{L}_{j-1}-h_j\left(L_j+{\mathit{U}}_j\right)-H_j\left(V_j+W_j\right)-Q_j$$E27

There are also in the literature some pseudo-transient mathematical modeling methods that turn the system algebraic eq. (AE) model that describes the process in steady state into a system of differential-algebraic equation, a mathematical model that is statically equivalent to original models [28, 86].

For some cases, molar vapor and liquid flow rates can be considered as constant, in addition to the assumption that there is no heat loss, which eliminates the need for energy balances. Binary distillation is known as the McCabe-Thiele method and can be applied graphically [87]. This principle can be applied to more than two components, and it is known as non-heat effect.

A good mathematical model must represent appropriately the real behavior of the process, and to get this behavior is necessary also a good numerical method for solving the mathematical model. The simulation is the translation of the numerical results obtained by this solution procedure to physical meaning. The next section discusses some methods and strategies used in the simulation of steady-state of reactive distillation.

5. Simulation

The solution process of a mathematical model is known as simulation (nowadays, simulation is synonymous with computational simulation) [68]. Mathematical models used for representing the behavior of some properties relevant to chemical processes are increasingly large due to the inclusion of so many details that aim to reach a realistic description of processes [88].

The complexity of a steady-state reactive distillation model is not due to mass and energy conservation laws, and it comes mainly from constitutive relations of chemical equilibrium or kinetics, mass transfer (in a more simplified way, phase equilibrium), and relations to compute enthalpy [68]. Due to the complexity of the models, it is necessary to apply robust methods or good solution strategies to make it possible to reach solution convergence.

The development of phenomenological mathematical model for a chemical process operating in steady state usually results in a system of non-linear algebraic equations that must be solved iteratively. The model is composed of N algebraic equations, as shown in Eq. (28), and the solution is obtained with respect to a set of N unknown variables, X:

$$F\left(X\right)=\left[\begin{array}{c}{f}_1\left(X\right)\\ f_2\left(X\right)\\ ⋮\\ f_n\left(X\right)\\ \end{array}\right]=\left[\begin{array}{c}{f}_1\left(x_1,x_2,x_3,\dots ,x_N\right)\\ f_2\left(x_1,x_2,x_3,\dots ,x_N\right)\\ ⋮\\ f_n\left(x_1,x_2,x_3,\dots ,x_N\right)\end{array}\right]$$E28

$$X=\left[\begin{array}{c}{x}_1\\ x_2\\ ⋮\\ x_N\end{array}\right]$$E29

The most popular methods for solving this kind of problem are Newton-Raphson and similar ones because these methods present fast convergence when good initial guesses are given, but these methods can fail when no good initial estimates are available, or there are several solutions [75, 89]. The search direction of the Newton-Raphson method is given by the Jacobian matrix, given by Eq. (30):

$$J\left(X\right)=\left[\begin{array}{cccc}\frac{\partial f_1\left(X\right)}{\partial x_1}& \frac{\partial f_1\left(X\right)}{\partial x_2}& \cdots & \frac{\partial f_1\left(X\right)}{\partial x_N}\\ & & & \\ \frac{\partial f_2\left(X\right)}{\partial x_1}& \frac{\partial f_2\left(X\right)}{\partial x_2}& \cdots & \frac{\partial f_2\left(X\right)}{\partial x_N}\\ & & & \\ ⋮& ⋮& \ddots & ⋮\\ & & & \\ \frac{\partial f_N\left(X\right)}{\partial x_1}& \frac{\partial f_N\left(X\right)}{\partial x_2}& \cdots & \frac{\partial f_N\left(X\right)}{\partial x_N}\end{array}\right]$$E30

Using the Newton-Raphson method is a necessary effort to obtain the derivatives analytically, which prevents its use, so it is common to use numerical derivatives. Computer algebra like Maple and Mathematica can make the analytical derivatives obtaining process easier, but these programs tend to be slow when compared to Fortran or C [88].

The iterative process of the Newton-Raphson method to approximate the solution of $k+1$ iteration depends on the solution obtained in the last iteration (k), as shown in Eq. (31):

$$X^{\left(k+1\right)}=X^{\left(k\right)}+\mathrm{\Delta }{X}^{\left(k\right)}$$E31

where $\mathrm{\Delta }{X}^{\left(k\right)}$ is obtained by the solution of the linear system of Eq.(32):

$$J\left(X^{\left(k\right)}\right)\mathrm{\Delta }{X}^{\left(k\right)}=-F\left(X^{\left(k\right)}\right)$$E32

Newton-Raphson method is fast and robust near the solution; however, its performance strongly depends on the choice of good initial estimates [86]. An improvement of the Newton-Raphson method can be obtained by relaxation [1], like the method proposed by [90], Eq. (33), that uses a relaxation parameter s, which is different from the unity when the solution tends to diverge, that is Eq. (34) is not satisfied:

$$X^{\left(k+1\right)}=X^{\left(k\right)}+s_i^{\left(k\right)}\mathrm{\Delta }{X}^{\left(k\right)}$$E33

$$\sqrt{{\sum }_{j=1}^N{f}_j^2\left(X^{\left(k\right)}+s_1^{\left(k\right)}\mathrm{\Delta }{X}^{\left(k\right)}\right)}<\sqrt{{\sum }_{j=1}^N{f}_j^2\left(X^{\left(k\right)}\right)}$$E34

where

$$s_2^{\left(k\right)}=\frac{\sqrt{1+6\eta }-1}{3\eta }$$E35

$$\eta =\frac{{\sum }_{j=1}^N{f}_j^2\left(X^{\left(k\right)}+s_1^{\left(k\right)}\mathrm{\Delta }{X}^{\left(k\right)}\right)}{{\sum }_{j=1}^N{f}_j^2\left(X^{\left(k\right)}\right)}$$E36

Homotopy method can guarantee that the approximate solution is reached if it is chosen an adequate auxiliary homotopy function, G(X), presented in Eq. (37). The solutions to the auxiliary function may be easily guessed/given/known [89]:

$$G\left(X\right)=\left[\begin{array}{c}{g}_1\left(X\right)\\ g_2\left(X\right)\\ ⋮\\ g_N\left(X\right)\end{array}\right]$$E37

Then, we define the homotopy function as presented in Eq. (38):

$$H\left(X,\lambda \right)=\left(1-\lambda \right)G\left(X\right)+\lambda F\left(X\right)$$E38

where $\lambda$ is the homotopy parameter that varies from 0 to 1, being that, for $\lambda =0$$H\left(X,0\right)=G\left(X\right)$ and $\lambda =1$$H\left(X,1\right)=F\left(X\right)$.

The goal is to solve $H\left(X,\lambda \right)=0$ instead of $F\left(X\right)=0$, but in this case, there are $N+1$ unknowns and N equation thus for the solution to be possible, the value of $\lambda$ is fixed, starting in 0, after solving the homotopy function with a fixed value of lambda, its value is sequentially changed until the value $\lambda =1$ is reached. The solution obtained with the last value of $\lambda$ is used as an initial guess to the next step to avoid the situation of divergence.

The homotopy path from zero to one can be followed by taking sequential steps of small values of $\lambda$ and by applying a Newton-Raphson method at each step to track the curve; this is called piecewise linear. However, the parametrized $\lambda$ might lead to poor performance and turning points can be encountered, leading to a failure in the solution procedures. To make the method more efficient, the homotopy function can be turned into an initial value problem (IVP) of a system of ordinary differential equations, and solved by an appropriate numerical method for IVP. The variation in the $\lambda$ value can have the use of a predictor-corrector continuation method applied [91]. Homotopy-continuation methods are mostly applied in problems that present high difficulty solved by other methods, like Newton-Raphson [92].

Some works have adopted a procedure to apply the Newton-Raphson method to solve, all together, the equations obtained from mass and energy balances, phase equilibrium equations, chemical equilibrium or rates of reaction equations, and any additional equation [88]. Other works have applied the Newton-Raphson method with numerical evaluation of the Jacobian matrix to algorithms that divide the models into subsets of equations [59].

A solution strategy, adapted from conventional distillation, for solving the steady-state reactive distillation modeling is called “tearing equations,” shown in Figure 10a, in which, after some rearrangements into the model equations, the solution is divided into many steps, a solution of a subset of unknowns in each step, in a way that it is possible to avoid the solution of non-linear equations sets for all unknowns, mostly the subsets of equations are composed only of linear equations because the other unknowns (not being calculated in the current step) are considered constant [75, 93, 94]. This strategy can make the implementation of the solution algorithm simpler and the solution faster. The subset of equations sequence is arranged in a generic way because the sequence can vary from one method to another and mainly because the strategy of subset grouping can be different, but the subset of unknowns must be contained in the subset of equation, for example, it is possible to use the component mass balance to calculate the mass fraction, use the mass balance in each stage to calculate the liquid stream leaving each stage, the energy balance to calculate the vapor stream leaving each stream, use the chemical equilibrium or reaction rate to calculate the extent of reaction, the summation of molar fractions combined with phase equilibrium to calculate the temperatures.

Figure 10.

Flowchart of the solution algorithms.

An improvement in the “tearing equations” methodology is to add an inner iterative process, Figure 10b, in which just some of the subsets are solved. This can reduce the number of iterations in the outer loop and make the solution even faster, called “inside-out” [74, 95, 96]. The subsets of equations for these methods normally are tridiagonal systems of linear equations, which can be solved by, for example, Thomas algorithm to reduce the computer time needed.

Similar to tearing methods, the simultaneous correction method also separates the equations in many subsets, but the equations are first linearized [96]. Most algorithms for solving reactive distillation steady state models are adaptations of well-known algorithms for solving conventional distillation steady state models, like the Naphtali-Sandholm method; in this method, equations are grouped by stages rather than by components and solved by the Newton-Raphson method [46, 97, 98]. Also, dynamic fuzzy neural networks have been implemented to solve mathematical models of steady-state engineering process that are represented by a system of non-linear equations [99].

The model equations can be implemented using a programming language (like Fortran or C), model builders (like gPROMS, an equation-oriented simulator) or a sequential modular simulator. The use of sequential modular simulators like Aspen Plus, Aspen HYSYS, CHEMCAD, ProSimPlus, etc., is common for analyses of reactive distillation systems [2, 4, 10, 13, 21, 22, 23, 26, 27, 36, 37, 43, 44, 46, 61, 65, 67, 100, 101, 102, 103, 104, 105, 106, 107, 108, 109, 110, 111, 112, 113, 114, 115]. Due to a wide variety of libraries for equipments that they include, thus allowing the users to simulate almost any process plant [116]. Most research about the conceptual design of the reactive distillation process is still carried out on commercial simulators such as Aspen Plus for steady-state design [14].

Aspen (Advanced System for Process Engineering) Plus has been used for many recent process modeling involving physical chemistry, chemical thermodynamics, mass and energy balances, and chemical reaction engineering [117, 118, 119]. This software offers a completely integrated solution to chemical process industries, making it possible to use it in many tasks of process engineering, from design to capital analysis. It has built-in model libraries from many chemical process, including distillation and reactive distillation [44]. It has a strong thermodynamics database and robust numerical solvers [120].

There are in the steady state simulation Aspen Plus robust modules for equilibrium stage and rate-based models. RADFRAC module is based upon a rigorous equilibrium stage model for solving MESH equations, while RATEFRAC is a rate-based model [32, 46, 67, 108, 109, 114, 119, 121, 122, 123, 124].

HYSYS Software is a powerful engineering simulation tool comprised of various components that provide an extremely powerful approach to steady-state modeling [125].

Although the steady-state modeling and simulation of reactive distillation column has been studied for decades, there are still some challenges and opportunities in this research subject, among which some are presented in the next section.

6. Challenges and opportunities

One can see that the reactive distillation process has received many proposals for modification in the hardware structure, and many others can be proposed. Also, there are many study cases reported in the literature, and many others will be studied. However, the greatest challenge that we can see is the solution procedure of the steady-state model.

The availability of many commercial simulators in which the user does not need to know what equations are solved in the simulation is a very interesting tool because the researcher can focus only on the process being studied, but it is important to know what is being solved in a computational simulation to correctly understand the results being reached. Also, the modular sequential simulators can be used to simulate almost any process plant, but there is a disadvantage that arises when a specific process, whose behavior cannot be accurately represented by any of these general models, is present in a process plant. In this case, if the researcher is not capable of obtaining and solving a mathematical model by using programming languages, he/she cannot do any simulation about the process, and there are some characteristics and strategies in the solution procedures that can be improved.

The solution of a steady-state model of a reactive distillation problem involves the solution of a system of highly non-linear algebraic equations [75, 93]. The solution of a system of non-linear equations is very hard and usually needs good initial guesses in the way the method presents convergence [126]. For a steady-state reactive distillation model, a good initialization strategy is required due to the highly non-linear behavior of the thermodynamic model and the kinetic equations [14]. There are very few works that have some propose to obtain the initial estimates of a steady state of a reactive distillation process [75, 93]. Some authors use the solution obtained in other works as initial estimates [88, 127] or use the solution of a more simple model and the results [85].

A methodology for obtaining an initial guess is provided for the solution of steady-state modeling of equilibrium distillation by [94] and steady-state modeling of equilibrium reactive distillation by [75, 93]. These authors used a non-heat effect for flow rates, a linear profile for temperature and a quadratic profile for the extent of the reaction [75, 93]. The quadratic profile for the extent of reaction with the maximum in the feed stage was also used by [74]. The variables already obtained as initial guess are enough to start the numerical solution procedure that uses a tearing equations strategy. For the initial guess of the temperature profile [74, 95], a linear interpolation for the temperature considering the bubble point in the condenser and the dew point in the reboiler dew point for an average composition of all feed streams. The authors in [95] do not present details on the equations, but they also cited that vapor and liquid-phase fractions are initialized by assuming constant molar overflow in the column and also considered ideal equilibrium ratios and non-reactive systems for other initial guesses.

The existence of multiple steady state (MSS) makes the importance of initial estimates even greater. Simulation studies suggest that reactive distillation processes exhibit complex multiple steady-state behavior in many cases because of many physicochemical phenomena interacting deliberately [33, 34, 47, 72, 82, 85, 128, 129, 130, 131]. There is the possibility of using commercial packages to identify the MSS, like ASPEN PLUS [37, 100, 132, 133], HYSYS [18, 73, 78, 134, 135], and SPeedUp™ [136]. Reactive distillation can exhibit input multiplicity and output multiplicity. The input multiplicity occurs when the same output is obtained by a multiple set of inputs, and an output multiplicity occurs when the same input results in multiple sets of outputs. Among those can exist stable and unstable states. To design a reactive distillation process, it is important to discover all steady states. The multiple solutions can be reached by performing a sensitivity analysis in some parameters or varying the initial estimates [37, 137, 138].

[139] proposed a method to find all steady-state solutions of the distillation column. In this work, it is stated that there is no need for initial estimates. The method requires a specific but fairly general block-sparsity pattern, with a linear growth in the computation effort. The algorithm requires that its input system of equations has been permuted in such a way that the Jacobian is in lower block Hessenberg form. The variables are partitioned into subvectors, and the equations are divided into equation subsets in a way that only variables from the subvectors can appear in the subsets of equations, but one variable can appear in more than one subvector. The original problem is reduced to much smaller subproblems. There is some lack of information in that proposal because it is stated that there is no need for initial estimates, but some subset of equations must be non-linear, which makes it necessary to have initial estimates.

Homotopy method can be applied without the supply of good initial estimates, but it turn the computational effort much greater, because the system of equations must be modified and solved many times.

After all, one can see that the best method, considering implementation difficulty, computational effort, and success in obtaining the solution, is the use of methods like Newton-Raphson or Broyden with good initial estimates exploring the sparsity of the model like used in simultaneous correction, Naphtali-Sandholm method, tearing equations strategy or inside-out algorithm, or even solving all equation together.

Heuristic methods, for global optimization, have been receiving some interest in the last years [140]. Advances in computing and information technology allow chemical engineers to solve complex design problems [25] in such a way the genetic algorithm has been used to solve simple problems that are a system of algebraic equations [141], for one just can write the objective function as the sum of the absolute value (or the square) of each algebraic function $\left(F\left(X\right)=0\right)$:

$$F_{\mathrm{obj}}=\sum _i\left|f_i\left(X\right)\right|$$E39

$$F_{\mathrm{obj}}=\sum _i{\left(f_i\left(X\right)\right)}^2$$E40

The use of heuristic optimization methods may require much computational time to solve large systems of algebraic equation, but it can quickly find some good candidates as initial estimates for a fast deterministic numerical solution method like Newton-Raphson or Broyden. And, the nature of heuristic methods can be used to generate different initial estimates to be used in the study of multiple steady states. Also, the method proposed by [75, 93] can be used as one of the possible solution at the start of the heuristic algorithm.

Another important opportunity is about the representation of the generic flows in a column stage because it is a smart choice to make the stage configuration representation as generic as possible in a way that just one model can be used to represent all the stages in a column, being necessary only the specification of the streams connections, this also makes it easier to write the computational code. Many authors used their models generic stages in a way that can be applied to the stage presented in Figure 6. That is a good choice, but its representation and resulting material and energy balances cannot be used, for example, in all elements of columns that have side reactor, side columns, pump around, etc. So, here we are proposing a more generic stage schematic model, presented in Figure 11, that can be used to represent all these. The change is that the side liquid and flow rates can be sent out of the process (product output stream) for side stages (side reactor, heat exchanger and side column stage) or even for another stage in the same column. It turns the model more complex and with many more specifications needed, but makes possible to generate all material and energy balances in the same way, even for the side reactor, pump around, or side columns, which require a considerable amount of attention on the stages numbering and specifications. Even the bypass efficiency can be modeled with this flow scheme, but the column stages must be numbered sequentially from top to bottom, and the other elements (side column, side reactors, pump around, etc.) must be numbered out of this range.

Figure 11.

Schematic of a more generic stage in a distillation column.

For this schematic representation, the global mass balance is represented by Eq. (41):

$$F_j+R_j+V_{j+1}+\sum _{\begin{array}{c}l=1\\ l\ne j\end{array}}^n{\omega }_{lj}{W}_l+L_{j-1}+\sum _{\begin{array}{c}l=1\\ l\ne j\end{array}}^n{\lambda }_{lj}{\mathit{U}}_l-\left(L_j+{\mathit{U}}_j\right)-\left(V_j+W_j\right)=0$$E41

$$F_j^L+F_j^V+R_j++\sum _{\begin{array}{c}l=1\\ l\ne j\end{array}}^n{\tau }_{\mathit{lj}}{V}_l+\sum _{\begin{array}{c}l=1\\ l\ne j\end{array}}^n{\sigma }_{lj}{V}_l-L_j-V_j-P_j^L-P_j^V=0$$E42

where ${\omega }_{lj}$ is the fraction of side vapor flow rate that leaves stage l and goes into stage j, ${\lambda }_{lj}$ is the fraction of side liquid flow rate that leaves stage l and goes into stage j, ${\omega }_{jj}$ is the fraction of side vapor flow rate that leaves the process (product output stream) from stage j and ${\lambda }_{jj}$ is the fraction of side liquid flow rate that leaves the process (product output stream) from stage j.

In another generic stage modeling we are proposing, as shown in Figure 12, it is considered separated liquid and vapor feed streams ($F_j^L$ and $F_j^V$), as same as product streams ($P_j^L$ and $P_j^V$) in a way that it is not considered side stream, like in the previous model, but the liquid and vapor outputting stage j can go into any other stage. Thus, the stages can be numbered freely, even column side elements. For this schematic representation, the global mass balance is represented by Eq. (42).where ${\tau }_{lj}$ is the fraction of vapor flow rate that leaves stage l and goes into stage j, ${\sigma }_{lj}$ is the fraction of liquid flow rate that leaves stage l and goes into stage j, $F_j^L$ is the feed flow rate in the liquid phase in stage j, $F_j^V$ is the feed flow rate in the vapor phase in stage j, $P_j^L$ is the product flow rate in the liquid phase from stage j and $P_j^V$ is the product flow rate in the vapor phase from stage j.

Figure 12.

Schematic of another generic stage in a distillation column.

7. Conclusions

This work was carried out as a systematic literature review about modeling and simulation of steady-state reactive distillation processes. Using the Methodi Ordinatio, a portfolio was obtained composed of 144 papers because it eliminated the articles that resulted in a negative Ordinatio index. With the bibliometric article’s data, it was possible to find the countries, journals, and authors with the most publications in the portfolio, and it was also possible to see the temporal evolution of the number of articles. Reactive distillation is an evolution of process in which the reaction and separation take place in different units, but there are many modifications in the configuration of the reactive distillation columns aiming to improve the product quality or reduce the capital or operational costs, being that the reactive dividing wall column can be highlighted and there are many modifications in a specific case that are being evaluated. A steady-state simulation is a fundamental tool for evaluating a change in the process without effectively doing the modification in practice, lowering the cost in analysis, and doing the modification only if it is favorable. Lately, these simulations are mainly conducted by commercial simulators, like Aspen Plus, which are excellent tools that greatly reduce the investigation time for the process modification. There are some strategies that can be improved in the modeling and in the solution of the models, mainly about the initial estimates, a very important component in the solution process that has few proposes in the literature; it also can be possible to propose some more general flow schemes in and out of a stage in a way that only one general balance can be applied in the main column stages and side elements like side reactor, pump around and side columns.

Conflict of interest

The authors declare no conflict of interest.

Abbreviations

AE	algebraic equation
CD	catalytic distillation
DAE	differential-algebraic equation
DWC	dividing wall column
EQ	equilibrium
eRDWC	enzymatic reactive dividing wall column
HIDiC	heat-integrated distillation column
HIRDiC	heat-integrated reactive distillation column
HP	heat pump
HPRD	heat pump reactive distillation
IVP	initial value problem
MSS	multiple steady states
NEQ	Non-equilibrium
PSD	pressure swing distillation
PSRD	pressure swing reactive distillation
RD	reactive distillation
RDC	reactive dividing column
RDWC	reactive dividing wall column
REDWC	reactive-extractive dividing wall column
SRC	side reaction column
SS	steady state
TCDS	thermally coupled distillation sequences
VC	vapor compression column
RC	vapor recompression column
Symbols
$a_{ij}$	activity of component i on stage j
$E_{ML}$	Murphree tray efficiency for liquid phase
$E_{MV}$	Murphree tray efficiency for vapor phase
F	system of non-linear algebraic equation
$F_j$	feed flow rate to stage j
$F_j^L$	feed flow rate to stage j in the liquid phase
$F_j^V$	feed flow rate to stage j in the vapor phase
G	Gibbs-free energy
$G^{°}$	Gibbs-free energy in standard state
$h_j$	enthalpy of liquid stream flow from stage j
$h_j^F$	enthalpy of feed stream flow to stage j
$H_j$	enthalpy of vapor stream flow from stage j
J	Jacobian matrix
$K_{ij}$	liquid-vapor equilibrium relation for the component i on stage j
$k_k$	kinetic constant of reaction k
$K_k$	chemical equilibrium constant of reaction k
$L_j$	liquid flow rate from stage j
$L_j^{\ast }$	liquid flow rate from stage j in equilibrium
n	total number of stages
$n_{ij}$	quantity of component i in stage j after the reaction reached the equilibrium (for the extent of reaction calculation)
$n_{ij}^0$	total amount of the inlet quantity in the liquid phase (from feed stream or liquid output from stage $j-1$) of component i in stage j (for the extent of reaction calculation)
m	total number of components
$Q_j$	heat load from stage j
P	pressure
$P_j$	pressure in stage j
$P_j^L$	product flow rate from stage j in the liquid phase
$P_j^V$	product flow rate from stage j in the vapor phase
$P^{°}$	pressure in standard state
$P^{\mathrm{sat}}$	saturation pressure
$r_{ij}$	generation rate of component i on stage j
R	ideal gas constant
$R_j$	total rate of generation on stage j
T	temperature
$T_j$	temperature of stage j
$u_j$	internal energy of liquid phase in stage j
${\mathit{U}}_j$	liquid side flow rate from stage j
$V_j$	vapor flow rate from stage j
$V_j^{\ast }$	vapor flow rate from stage j in equilibrium
$V_{hj}$	liquid molar holdup in stage j
$V_i^L$	component i volume in liquid phase
$W_j$	vapor side flow rate from stage j
$x_{ij}$	liquid mole fraction of component i in stage j
$x_{ij}^{\ast }$	equilibrium liquid mole fraction of component i in stage j
X	set of unknowns
$y_{ij}$	vapor mole fraction of component i in stage j
$y_{ij}^{\ast }$	equilibrium vapor mole fraction of component i in stage j
$z_{ij}$	feed mole fraction of component i in stage j
${\alpha }_i$	relative volatility of component i
${\gamma }_{ij}$	activity coefficient of component i on stage j
${\delta }_j$	binary variable used to identify reactive stages, equal to 1 if the reaction occurs in stage j, 0 otherwise
${\epsilon }_j$	bypass efficiency of stage j
$\eta$	reduction in the Euclidean norm of the set of non-linear functions between two steps of a Newton-Raphson relaxation method
${\eta }_j^L$	ratio between the liquid side stream and the liquid stream outputting stage J and inputting the near stage $j+1$
${\eta }_j^V$	ratio between the vapor side stream and the vapor stream outputting stage J and inputting the near stage $j-1$
${\lambda }_{lj}$	fraction of side liquid flow rate that leaves stage l and goes into stage j${\nu }_{ki}$ stoichiometric coefficient of component i in reaction k
${\xi }_{kj}$	extent of reaction k in stage j
${\sigma }_{lj}$	fraction of vapor flow rate that leaves stage l and goes into stage j
${\tau }_{lj}$	fraction of liquid flow rate that leaves stage l and goes into stage j
${\widehat{\varphi }}_{ij}$	fugacity coefficient of component i on stage j
${\varphi }_i^{\mathrm{sat}}$	fugacity coefficient of pure component i in saturation state
${\omega }_{lj}$	fraction of side vapor flow rate that leaves stage l and goes into stage j