Inside–Out Method for Simulating a Reactive Distillation Process

Reactive distillation is a technical procedure that promotes material strengthening and its simulation plays an important role in the design, research, and optimization of reactive distillation. The solution to the equilibrium mathematical model of the reactive distillation process involves the calculation of a set of nonlinear equations. In view of the mutual influence between reaction and distillation, the nonlinear enhancement of the mathematical model and the iterative calculation process are prone to fluctuations. In this study, an improved Inside–Out method was proposed to solve the reaction distillation process. The improved Inside–Out methods mainly involved—(1) the derivation of a new calculation method for the K value of the approximate thermodynamic model from the molar fraction summation equation and simplifying the calculation process of the K value, as a result; and (2) proposal for an initial value estimation method suitable for the reactive distillation process. The algorithm was divided into two loop iterations—the outer loop updated the relevant parameters and the inside loop solved the equations, by taking the isopropyl acetate reactive distillation column as an example for verifying the improved algorithm. The simulation results presented a great agreement with the reference, and only the relative deviation of the reboiler heat duty reached 2.57%. The results showed that the calculation results were accurate and reliable, and the convergence process was more stable.


Introduction
Reactive distillation combines the chemical reaction and distillation separation process in a single unit operation [1]. As reactive distillation has the advantages of improving reaction conversion rate and selectivity and reducing costs, the research on design and optimization of the reactive distillation processes has increased rapidly in recent years [2]. Now reactive distillation technology has been successfully applied to the areas of esterification, etherification, hydrogenation, hydrolysis, isomerization, and saponification [3].
The production of methyl tert-butyl ether (MTBE), methyl tert-pentyl ether (TAME), and ethyl tert-butyl ether (ETBE), etc., adopts the reactive distillation technology, and the reaction conversion rate was significantly higher than the fixed-bed reactors [4]. Using the reactive distillation technology to produce methyl acetate can reduce investment and energy costs to one-fifth, and the reaction conversion rate is close to 100% [5]. Zhang et al. proposed a novel reaction and extractive distillation (RED) process and used it for the synthesis of isopropyl acetate (IPAc). The results showed that the purity of IPAc reached 99.5% [6].

Mathematical Model
The schematic diagram of the reaction distillation column is shown in Figure 1. The whole column is composed of N stages. The condenser is regarded as the first stage and the reboiler is regarded as the last stage. There can be a liquid phase product, a vapor phase product, and a free water product at the top of the tower, and a liquid phase product at the bottom of the column. Chemical reactions can take place anywhere in the column. The general model of the jth theoretical stage is shown in Figure 2, where F j is the feed flow rate of stage j, L j is the liquid flow rate outputting stage j, and inputting stage j + 1, V j is the vapor flow rate outputting stage j and inputting stage j − 1, U j is the liquid side flow rate outputting stage j, W j is the vapor side flow rate outputting stage j, Q j is the heat duty from stage j, R r,j is the reaction extent of the rth reaction on stage j. The schematic representation from Figure 2 to represent all stages are shown in Figure 3. stage is shown in Figure 2, where Fj is the feed flow rate of stage j, Lj is the liquid flow rate outputting stage j, and inputting stage j + 1, Vj is the vapor flow rate outputting stage j and inputting stage j − 1, Uj is the liquid side flow rate outputting stage j, Wj is the vapor side flow rate outputting stage j, Qj is the heat duty from stage j, Rr,j is the reaction extent of the rth reaction on stage j. The schematic representation from Figure 2 to represent all stages are shown in Figure 3.
The stages of the column are assumed to be in equilibrium conditions, the vapour and liquid phases leaving the stage are assumed to be in thermodynamic equilibrium, the stage pressure, temperature, flow, and composition are assumed to be constant at each stage. Five sets of equations are used to describe the equilibrium state of the stage-the material balances, the phase equilibrium, the molar fraction summation, the enthalpy balance, and the chemical equilibrium equations.
A set of highly nonlinear equations was obtained by modeling the steady state process of the reactive distillation column. The non-linearities were in the phase equilibrium, enthalpy balances, and chemical equilibrium equations. The solution of model equations is very difficult and the convergence depends strongly on the goodness of the initial guesses. After choosing appropriate iterative variables and providing the initial value of variables, the iterative calculation process can start. When the iterative calculation converges, the stage temperature, flow, and composition profiles can be obtained.   Uj is the liquid side flow rate outputting stage j, Wj is the vapor side flow rate outputting stage j, Qj is the heat duty from stage j, Rr,j is the reaction extent of the rth reaction on stage j. The schematic representation from Figure 2 to represent all stages are shown in Figure 3.
The stages of the column are assumed to be in equilibrium conditions, the vapour and liquid phases leaving the stage are assumed to be in thermodynamic equilibrium, the stage pressure, temperature, flow, and composition are assumed to be constant at each stage. Five sets of equations are used to describe the equilibrium state of the stage-the material balances, the phase equilibrium, the molar fraction summation, the enthalpy balance, and the chemical equilibrium equations.
A set of highly nonlinear equations was obtained by modeling the steady state process of the reactive distillation column. The non-linearities were in the phase equilibrium, enthalpy balances, and chemical equilibrium equations. The solution of model equations is very difficult and the convergence depends strongly on the goodness of the initial guesses. After choosing appropriate iterative variables and providing the initial value of variables, the iterative calculation process can start. When the iterative calculation converges, the stage temperature, flow, and composition profiles can be obtained.   The stages of the column are assumed to be in equilibrium conditions, the vapour and liquid phases leaving the stage are assumed to be in thermodynamic equilibrium, the stage pressure, temperature, flow, and composition are assumed to be constant at each stage. Five sets of equations are used to describe the equilibrium state of the stage-the material balances, the phase equilibrium, the molar fraction summation, the enthalpy balance, and the chemical equilibrium equations.

Inside-Out Method
The Inside-Out method divides the mathematical model into an inside and outer loop for iterative calculation, and uses two sets of thermodynamic models, a strict thermodynamic model and an approximate thermodynamic model (K value model and enthalpy value model). The simple A set of highly nonlinear equations was obtained by modeling the steady state process of the reactive distillation column. The non-linearities were in the phase equilibrium, enthalpy balances, and chemical equilibrium equations. The solution of model equations is very difficult and the convergence depends strongly on the goodness of the initial guesses. After choosing appropriate iterative variables and providing the initial value of variables, the iterative calculation process can start. When the iterative calculation converges, the stage temperature, flow, and composition profiles can be obtained.

Inside-Out Method
The Inside-Out method divides the mathematical model into an inside and outer loop for iterative calculation, and uses two sets of thermodynamic models, a strict thermodynamic model and an approximate thermodynamic model (K value model and enthalpy value model). The simple approximate thermodynamic model was used for frequent inside loop calculations to reduce the time it takes to calculate the thermodynamic properties. The strict thermodynamic model was used for the outer loop calculations, and the parameters of the approximate thermodynamic model were corrected at the outer loop.
Defining the inside loop variables-relative volatility α i,j , stripping factors S b,j of the reference components, liquid phase stripping factors R L,j , and the vapor phase stripping factors R V,j .
The inside loop uses the vapor and liquid phase flow rates of each component instead of the composition, for the calculations. Their relationship with the vapor and liquid phase composition and flow is as follows: According to the inside loop variables defined by the mathematical model of the Inside-Out method, the material balance Equation (M) and the phase equilibrium Equation (E) can be rewritten into the following relationship.
Material Balance: Phase balance:

Strict Thermodynamic Model
The strict model is mainly used for the outer loop to correct the parameters in the approximate thermodynamic model, including the calculation of the phase equilibrium constant K value, and the vapor-liquid enthalpy difference, as follows: Define a new variable for the K value of the reference component, K b,j , and the reference component b is an imaginary component. K b,j was calculated by the following formula, where ω i,j is a weighting factor.
The relationship between K b,j and the stage temperature T j is related using Equation (17), T* is the reference temperature.
The value of the parameter B j can be obtained by the differential calculation of 1/T by the K b , and the temperatures T j−1 and T j+1 of the two adjacent stages are selected as T 1 and T 2 .
2 Enthalpy model The calculation of the enthalpy difference is the main time-consuming part of the entire enthalpy calculation, so the calculation of the enthalpy difference was simplified in the approximate enthalpy model, and a simple linear function was used to fit the calculation of the enthalpy difference.

Improved K-Value Model
The K-value model was improved by simplifying the calculation of the weighting factor. This was derived from the molar fraction summation equation. The derivation process is shown below.
The relationship between the K value of the reference component K b,j and the K value of the common component is: The weighting factors need to satisfy: The stage temperature needs to satisfy the following relationship: The introduction of new variables u i,j was defined as: The new relationship is: Since K b,j was calculated by K i,j through a weighting function, the trend with temperature changes was the same, so u j,i is a temperature-independent parameter, so at the stage temperature: The derivative of Φ T j and Φ T j at the stage temperature was equal to: The weighting function used the new calculation method:

Initial Value Estimation Method
The system of equations obtained from the modeling of the steady-state reactive distillation process was non-linear and was very difficult to solve. Therefore, it was important to provide good initial estimates, otherwise it would not be possible to reach a solution.
Before providing good initial estimates, all feed streams were mixed to perform chemical reaction equilibrium calculations to obtain a new set of compositions and flow rate, which were used to calculate the initial value of the variable.
Temperature: Calculate the dew and bubble point temperature by flash calculation, which was used as the initial temperature at the top and at the bottom of the tower. The initial temperature of the middle stages was obtained by linear interpolation.
Composition: The composition calculated by isothermal flash calculation was used as the initial value of the vapor-liquid composition of each stage.
Flow: According to the constant molar flow rate assumption and the total material balance calculation, the initial value of the vapor-liquid flow rate of each stage was obtained Reaction extent: The reaction extent calculated by the chemical reaction equilibrium was taken as the maximum value of the reaction extent in the reaction section. According to the feeding condition of the reaction section, the position with the maximum reaction extent was selected. If there was only one feed in the reaction section, the position of the feed stage was selected. If there were multiple feeds in the reaction section, the position of the feed stage near the reboiler was chosen. If there was no feed in the reaction section, the position of the reaction stage closest to the feed was selected. The initial value of the reaction extent of the other stages were calculated according to the maximum value and a certain proportion of attenuation. The calculation is shown below. Between the first reaction stage and the stage with the maximum reaction extent: Between the stage with the maximum reaction extent and the last reaction stage: where R max is the maximum value of the reaction extent in the reaction section, j m is the stage position where the maximum reaction extent occurs; j 1 is the first stage position in the reaction section; and j n is the last stage position in the reaction section.

Calculation Steps and Block Diagram
The reactive distillation column adopts an improved Inside-Out method for calculation. The main working idea is that the the outer loop used a strict thermodynamic model to calculate the phase equilibrium constant and the vapor-liquid phase enthalpy difference, and the results were used to correct the approximate thermodynamic model parameters and update the reaction extents of each stage. The inside loop uses an approximate thermodynamic model to solve the MESHR equation to obtain the stage temperature, flow rate, and composition. After the inside loop calculation converges or reaches the number of iterations, the model returns to the outer loop to continue the calculation. When the inside and outer loop converges at the same time, the calculation ends.
The calculation steps are shown below:

Given initial value
(1) According to the feed, pressure, and operation specifications, the initial values of the temperature, flow rate, composition, and the reaction extent are provided (Section 3.2).

Outer loop iteration
(2) Calculate the phase equilibrium constant K i,j with a strict thermodynamic model, and then calculate parameters K b,j , A, and B in the K-value model.
(3) Calculate the vapor-liquid phase enthalpy difference using a strict thermodynamic model, and fit the parameters c j , d j , e j , and f j in the approximate enthalpy model.
(4) Calculate the stripping factors S b,j , relative volatility α i,j , liquid phase stripping factors R L,j , and vapor phase stripping factors R V,j .
(5) Calculate the reaction extent according to the reaction conditions. The chemical reaction can be a specified conversion rate or a kinetic reaction, and the kinetic reaction rate is calculated by a power law expression.
Power law expression: T 0 is not specified: T 0 is specified: Processes 2020, 8, 604 9 of 17 (6) Outer convergence judgment conditions: 3. Inside loop iteration (7) According to the material balance Equation (13), a tridiagonal matrix was constructed, and then the liquid flow rate l i,j of each component was obtained by solving the matrix. Using the phase equilibrium Equation (14), the vapor phase flow rate v i,j of each component was obtained.
(8) The vapor-liquid flow rates V j and L j can be calculated by Equation (9), the vapor-liquid phase composition x i,j and y i,j can be calculated by Equation (8).
According to Equation (40), a new set of stage temperature can be calculated using the new K b,j values.
Now, the correction values of the vapor-liquid flow rate, composition, and temperature are obtained. They satisfy the material balance and phase equilibrium equation, but does not satisfy the enthalpy equilibrium equation. In the following, the iterative variable of the inside loop should be modified according to the deviation of the enthalpy balance Equation (10). Select lnS b,j as the inside loop iteration variable. If there are side products, the side stripping factors lnR L,j and lnR V,j should be added as the inside loop iteration variables. Taking the enthalpy balance equation as the objective function, the enthalpy balance equation of the first and last stage should be deleted and replaced with two operation equations of the reactive distillation column.
(11) Calculate the vapor-liquid phase enthalpy according to the simplified enthalpy model. (12) Calculate the partial derivative of the objective function to the iterative variable and construct the Jacobian matrix. The Schubert method was used to calculate the modified value of the inside iteration variable [22,23]. Damping factors can be used if necessary.
(13) Using the new inside iteration variable value, repeat steps (7) to (9) to obtain the new stage flow, composition, and temperature, and calculate the error of the objective function. If the error was less than the convergence accuracy of the inside loop, return to the outer loop to continue the calculation, otherwise repeat the calculation from Step (11) to Step (13). The calculation block diagram of the improved Inside-Out method is shown in Figure 4.

Example 1: Isopropyl Acetate
Taking the process of synthesizing isopropyl acetate (IPAc) with acetic acid (HAc) and isopropanol (IPOH) as an example, the esterification reaction is a reversible exothermic reaction, and

Example 1: Isopropyl Acetate
Taking the process of synthesizing isopropyl acetate (IPAc) with acetic acid (HAc) and isopropanol (IPOH) as an example, the esterification reaction is a reversible exothermic reaction, and the conversion rate is limited by the chemical equilibrium. Zhang et al. proposed reactive and extractive distillation technology for the synthesis of IPAc, and the extraction agent dimethyl sulfoxide (DMSO) was added in the tower to obtain high-purity IPAc [6]. Zhang et al. conducted a steady-state simulation of the reactive distillation column, and performed dynamic optimization control, based on the results. In this work, the improved Inside-Out method was used to solve the esterification process, and the steady-state simulation results were compared with that of Zhang et al. [6]. The chemical reaction equation of this esterification reaction was as follows: Kong et al. determined the kinetic model of the esterification reaction through kinetic experiments [24]. The kinetic expressions and parameters are as follows: where r is the reaction rate, mol·L −1 ·min −1 ; C i is the molar volume concentration of the i component, mol·L −1 ; k + , and kare the forward and reverse reaction rate constants; R is the gas constant, with avalue of 8.314 J·mol −1 ·K −1 ; and T is the temperature, in K. The improved Inside-Out method was used to simulate the isopropyl acetate process. The property method was non-random two liquid (NRTL). The operating conditions of the isopropyl acetate reactive distillation column is shown in Table 1. The schematic diagram of the isopropyl acetate reactive distillation column and the parameters of the feed stream are shown in Figure 5. Figure 6 is a comparison between the initial value of the reaction extent and the simulation result. It can be seen that the initial value and the simulation results were relatively close, which increased the possibility of convergence. The error behavior at each iteration, calculated by Equation (38), for the program using the improved Inside-Out method is shown in Figure 7. The outer loop calculation reached the convergence through 10 iterations, and the error decreased steadily, without oscillation. The schematic diagram of the isopropyl acetate reactive distillation column and the parameters of the feed stream are shown in Figure 5.  Figure 6 is a comparison between the initial value of the reaction extent and the simulation result. It can be seen that the initial value and the simulation results were relatively close, which increased the possibility of convergence. The error behavior at each iteration, calculated by Equation (38), for the program using the improved Inside-Out method is shown in Figure 7. The outer loop calculation reached the convergence through 10 iterations, and the error decreased steadily, without oscillation.   Table 2 shows the comparison between the simulation results of the isopropyl acetate reactive distillation column and the literature. It can be seen that the temperatures at the top and bottom of tower obtained in this work were similar to that in the literature, and the relative deviation was within 1%. The relative deviation of the condenser heat duty was 1.82% and the relative deviation of the reboiler heat duty was 2.47%.  Figure 6 is a comparison between the initial value of the reaction extent and the simulation result. It can be seen that the initial value and the simulation results were relatively close, which increased the possibility of convergence. The error behavior at each iteration, calculated by Equation (38), for the program using the improved Inside-Out method is shown in Figure 7. The outer loop calculation reached the convergence through 10 iterations, and the error decreased steadily, without oscillation.   Table 2 shows the comparison between the simulation results of the isopropyl acetate reactive distillation column and the literature. It can be seen that the temperatures at the top and bottom of tower obtained in this work were similar to that in the literature, and the relative deviation was within 1%. The relative deviation of the condenser heat duty was 1.82% and the relative deviation of the reboiler heat duty was 2.47%.   Table 2 shows the comparison between the simulation results of the isopropyl acetate reactive distillation column and the literature. It can be seen that the temperatures at the top and bottom of tower obtained in this work were similar to that in the literature, and the relative deviation was within 1%. The relative deviation of the condenser heat duty was 1.82% and the relative deviation of the reboiler heat duty was 2.47%.  Table 3 shows the comparison between the simulation results of the molar composition of the top and bottom products and the literature. It can be seen that the molar composition of the top and bottom product obtained in this work presented a great agreement with the literature, except for the smaller part, and the simulation results were accurate and reliable.

Example 2: Propadiene Hydrogenation
Taking the depropanization column in the propylene recovery system of the ethylene production process as an example. The Inside-Out method was used to simulate the depropanizer column, and the results were compared with the measured values. The top of the depropanizer column was equipped with a catalyst. In the reaction section, propadiene was hydrogenated into propylene and propane. The conversion rate of the first reaction was 0.485, the conversion rate of the second reaction was 0.015. The chemical reaction equation were as follows: The improved Inside-Out method was used to simulate the depropanization column. The property method was Soave-Redlich-Kwong (SRK). The depropanizer column had 3 feeds, and the feed stream parameters are shown in Table 4. The operating conditions of the depropanization column are shown in Table 5.  Table 6 shows the comparison between the simulation results of the depropanization column and the measured value. It can be seen that the temperatures in the top and bottom products, as calculated by the Inside-Out method were close to the measured value, and the relative error was about 1%. The liquid product flow rate at the top and bottom of the tower coincided with the measured value, only the relative deviation of the vapor product flow at the top of the tower was 2.74%. Table 7 shows the comparison between the mass composition of the top and bottom liquid products and the measured value. The simulation results of the product mass compositions were basically consistent with the measured value, except for the smaller composition values.

Conclusions
Based on the process principle of reactive distillation, a mathematical model of reactive distillation process was established. The improved Inside-Out method was provided in this work for the solution of a reactive distillation process. In view of the nonlinear enhancement of the reactive distillation mathematical model and the difficulty of convergence, the calculation of the K value of the approximate thermodynamic model was improved to simplify the calculation process. The initial value estimation method suitable for the calculation of reactive distillation was proposed, which increased the possibility of convergence.
The algorithm was verified by using an isopropyl acetate reactive distillation column and a depropanization column as examples. In Example 1, the reaction extent calculated by the initial value estimation method was close to the simulation results, it facilitated the convergence process of the solution algorithm. The simulation results obtained in this work were compared with Zhang et al. [6], presenting great agreement with the reference, only the relative deviation of the reboiler heat duty reached 2.57%. In Example 2, the simulation results of the depropanization column presented a good agreement with the measured value, except for the smaller composition value. The results showed that the improved Inside-Out method calculation results were accurate and reliable.