Dynamic Simulation, Parameter Optimization, and Control of a Reactive Distillation Column for Production of Isopropanol via Propylene Hydration

In this study, a reactive distillation column for production of isopropanol was investigated. Firstly, a dynamic model was developed for the process. The model of the process was then programmed, and the process simulated using a base case obtained from the literature. Results showed that distillate contained more than 58 mol% propylene-free isopropanol. In the next step, optimization of some operating variables was performed to maximize concentration of isopropanol in distillate with condenser temperature as constraint, which was considered to be above the freezing point of water. Several simulations were performed by changing operating parameters, and finally optimum isopropanol content in distillate was obtained above 58 mol%. Results of using classic controllers showed that PID controller had the best performance for both condenser temperature set-point tracking and disturbance rejection.


Introduction
Distillation is one of the important separation processes in chemical industry. Any improvement in this process can have great economic benefits. Reactive distillation (RD) is one of the major steps in separation history to achieve this goal. It is a combination of chemical reaction and physical separation within a single unit operation. Some of its benefits include increasing of reactants conversion and product selectively, heat integration, and reduction of fixed and operating costs 1 . The most important applications of RD include esterification 2,3 , hydration 4-6 , crude oil residue conversion 7 , etc. Reactive distillation was firstly considered by Smith 8 in 1980 for production of methyl tert-butyl ether (MTBE), and its first industrial application was introduced for production of MTBE in 1981 9 . Since then, RD has been used theoretically and experimentally for production of many components.
Isopropanol is an important chemical component with many applications in coatings for metals, painting, preparation of pesticides, pharmaceuticals, production of acetone, etc. 10 Isopropyl alcohol (IPA) industrial production methods include indirect hydration of propylene in the presence of homogeneous acidic catalysts 11 , and direct hydration of propene in the vicinity of a heterogeneous acidic catalyst 12, 13 . In indirect hydration of propene, firstly propene is reacted with 60 % sulfuric acid and converted to isopropyl hydrogensulphate (CH 3 CHCH 3 OSO 2 OH). In the next step, the produced intermediate is hydrated by water to produce IPA and H 2 SO 4 11 . In direct method, propene directly reacts with water in the presence of solid acidic catalysts 12, 13 . Indirect method needs more equipment than the direct method for producing IPA, and due to high corrosivity of sulfuric acid, currently it is replaced by the direct method. Direct method has higher efficiency and lower energy consumption. Catalysts of direct method include solid phosphoric acid 14 , synthetic zeolites 12 , and tungsten-based hetero-poly-acid 15 . Synthetic zeolites have advantages such as higher propylene conversion and lower energy consumption compared to other methods 13 .
In recent years, a new technology for production and separation of IPA using direct hydration of propene in a single catalytic distillation column has been investigated by some authors. Xu et al. 14 investigated production of IPA within a RD column using a zeolite catalyst. They developed their steady-state model in Aspen Plus™ and assumed reactions in chemical equilibrium. Wang and Wong 4 investigated the control of reactive distillation of IPA. The column had 26 ideal plates, and a reaction zone near condenser. They claimed that high conversions of propene were achieved and high-purity IPA exited from the reboiler as product. Niu and Rangaiah 16 investigated retrofitting of RD unit with acidic proton-exchange resin as catalyst, as well as the issue of excess propene in product and proposed two designs. In the first design, they used two RD columns where excess water entered the second column. In the second design, two RD columns received same feeds but with different flow rates. Their bottom products were then mixed, and entered an extractive distillation (ED) column. They used catalyst information only for economic calculations. Their results showed that the first design was more economic than the second one, with 14 % reduction in capital costs compared to a base case. Chua et al. 5 investigated the design and optimization of RD of isopropanol and proposed two designs. In the first design, they used propene in excess to consume water completely to avoid azeotropic conditions. In the second design, they used water in excess to consume propene and avoid its loss. Azeotropic mixture from the bottom of RD column was entered into an ED column. For economic calculations, they considered Amberlyst 36 as catalyst in the reaction zone of RD columns. Their results showed that the second design was more economic than the first one.
In this study, the dynamic model of the reactive distillation column of IPA production by direct hydration of propene on HZSM5 was created. The simulation was then run using a base case to reveal model performance. After validation of the results, optimization of some operating variables was performed by changing their values within predefined ranges. Finally, parameters of classic controllers were tuned using the dynamic response of the IPA composition in order to maintain IPA composition in distillate at its desired and optimum value.
To the best of our knowledge, no work has been done yet considering reaction kinetics in the reaction zone. In addition, all authors used commercial simulation software of Aspen Plus for process simulations, while we coded all the model equations, phase equilibrium calculations, and reaction kinetics in dynamic conditions, from scratch. Only Wang and Wong 4 performed process simulation dynamically, but they also only applied chemical equilibrium for the reaction zone. Another novelty of this study is the optimization of controller parameters, which has not been done by other authors. Fig. 1(a) shows the RD column used in this study for process modeling. Feed streams of pure water and propylene entered the rectifying section of the column at pressure of 20 bar and temperature of 380 K. The reaction zone was located between feed streams, and contained HZSM5 zeolite catalyst. Gas phase reaction was taking place between water and propylene in this zone to produce IPA. The rest of column was used for the separation of components as a normal distillation column. The reaction zone was assumed to be equivalent to three ideal plates.

Process modeling
For process modeling in dynamic mode, stages were numbered from bottom to top of the column. The following assumptions were used for process modeling: a) Pressure is constant at 20 bar within the column with no pressure drop, b) Liquid and gas phases on stages are at thermodynamic equilibrium, c) Gas phase on all stages is ideal, d) All stages are adiabatic except for condenser and reboiler, e) Holdups on stages are constant, f) Reaction only occurs in gas phase in reaction zone, g) Accumulation terms in total mass and energy balance equations are zero, h) Accumulation of components in gas phase is negligible compared to liquid phase, i) Condenser is partial.
The following reaction was considered to take place in the reaction zone: All parameters and variables are introduced in Symbols section. For process modeling, a system with one ideal stage of the column was considered, as shown in Fig. 1(b), and model equations were developed. The model was then used for all stages of the column, including condenser and reboiler.
Total mole balance for stage n assuming steady-state condition: where N+2 refers to the condenser.
Energy balance for stage n:

Mole balance for component i in stage n:
As mentioned previously, holdups were assumed to be constant. They were calculated using column internal diameter and weir height. Calculated holdups for trays were 0.141 m 3 assuming 2 m for column ID and 5 cm for weir height, and excluding 10 % of the tray cross-sectional area for downcomers 18 . Holdups for condenser and reboiler were considered to be three-fold of those for trays, i.e., 0.424 m 3 .
For calculating K-values, gas phase was assumed to be ideal, and activity coefficients of the components in the liquid phase were calculated using UNIFAC method.
Due to low critical temperature of propylene, generally, temperature in stripping section of the column and especially in the reboiler is above the critical temperature of propylene. For these conditions, Henry's law was applied for phase equilibrium calculations of propylene. Henry's law constant was calculated using method of Campanella et al. 19 To determine temperature changes in stages with time, the derivative of equation (8) was taken with respect to time. It was assumed that K-values are functions of temperature and liquid composition only. After simplification, the following equation was obtained: All required derivatives were obtained analytically.

Process simulation and optimization
After the development of dynamic model of the process, model equations were coded in MATLAB ® , and a base case was used for simulation. Initially, some reasonable values were assumed for liquid compositions and temperatures in all stages. Equations (4) and (6) were then solved analytically in each time step assuming steady-state condition in order to calculate liquid and vapor mole flow rates from the stages, as well as heat duties of condenser and reboiler. It was assumed that three parameters, including reflux ratio, distillate rate, and vapor fraction in condenser are known. In the next step, the system of component mole balance equations, (7), was solved numerically using an ordinary differential equation solver in each time step. Vapor mole fractions were calculated assuming thermodynamic equilibrium between phases in stages and using equation (10). As condenser was assumed to be partial, eq. (10) can also be applied to the condenser. In the reaction zone, generation and consumption of components were also included in calculations. Calculations were repeated until reaching final time determined for the simulation.
In this study, values of four factors with high impact on the performance of RD column, including reflux ratio, propene feed flow rate, water to propene feed ratio, and vapor fraction in condenser were examined and optimized by the method of one-factor-at-a-time 20 , and process simulation. Table  1 shows the parameters and their assessed values during optimization.

Process control
Classic controllers were designed and applied for the control of RD column at optimum conditions. Condenser temperature was considered as controlled variable, and reflux ratio as manipulating variable. For controller design, firstly, the response of condenser temperature to a step-change in reflux ratio was obtained. A first-order model with time-de-lay was then proposed for the response, and its parameters were estimated using process reaction curve 21 and curve fitting methods. The response of open-loop system to a step change in its input at t 0 = 0 is called the reaction curve. The step response of the above system can be described as: where A is the step-change amplitude. In the reaction curve method, a tangent line on reaction curve at inflection point is drawn. If intercept of straight line with horizontal axis is shown with t 1 then the parameters of the above system can be obtained using the following equations 21 : where t 2 is the location of the intercept of the tangent line with the horizontal line drawn at ultimate value of the response, y(t→∞) In the curve fitting method, τ d is estimated just like in the reaction curve method. It can then be shown that ln(dy /dt ) is a straight line with respect to t , where -1/τ p and ln(AK p /τ p ) are the slope and y-intercept of the line, respectively. Thus, K p and τ p can be found by linear regression. After system identification, the parameters of controllers were found using Ziegler-Nichols method 22 . Finally, the designed controllers were used for process control and their performances were compared.

Results and discussion
Vapor-liquid equilibrium results Fig. 2 shows T-xy diagram for the binary system of IPA and water at total pressure of 20 bar obtained using phase-equilibria software package written in this work. As may be seen, IPA and water formed an azeotrope at about 65 mol% of IPA, with mixture boiling point of 451.9 K. Therefore, complete separation of these two components in a single RD column could not be achieved and another extractive distillation column was required for this purpose. Fig. 3 shows triangular diagram for the ternary system of propene/IPA/water, which was obtained from Aspen Plus. As may be seen, this system only exhibited two-phase behavior, and the only azeotrope was at about 66 mol% of IPA and 34 mol% of water. Thus, the result validated the result of T-xy diagram calculated by our phase-equilibria package. A base case was considered for simulation and validation of the proposed mathematical model of the RD column based on the data from the literature 4 as well as different runs of the written programs. The values of operating and design parameters for base case are shown in Table 2. Temperature profiles for selected stages are shown in Fig. 4. Because the reaction is exothermic and the column was assumed to be adiabatic, the temperature in stage 25 was higher than the reboiler temperature. Minimum temperature was reported for condenser after 100 min, which was 257.6 K. Temperatures on other stages converged almost to a same value, except for stage 24, which was located in the reaction zone. As mentioned previously, the condenser was assumed to be partial for this process in order to avoid sharp decrease in temperature in the condenser.

F i g . 2 -T-xy diagram for the binary system of IPA (1) and
Water (2) 5 shows profiles of IPA compositions in selected stages of the RD column. As it is evident, IPA composition was near zero in the stripping section, but from reaction zone to the condenser it increased due to production by hydration reaction and separation by condenser. It seemed that IPA content in condenser was lower than in stage 27, which was due to the increase in propene content in cooling region, especially in the condenser (not shown here). However, since propene is a volatile component, it can be easily removed from IPA product outside of the column, and can be returned to the column. Obtained propene-free mole fraction of IPA in condenser was about 0.61, while it was 0.58 in stage 27 after 100 minutes of simulation.
Although good production and separation of IPA was achieved in the distillate, but due to the azeotropic condition of binary mixture of IPA and water, an additional extractive distillation column was required to obtain high purity IPA, which has also been investigated by Niu and Rangaiah 16 and Chua et al. 5 In addition, because condenser temperature is below the freezing point of water, the above result is not feasible, since coolant which is chilled water can't cool down the product to that temperature. To avoid this, condenser temperature has to be maintained above water freezing point. For this reason, and in order to obtain a good quality of output product, operating parameters of the process have to be optimized.

Optimization of IPA reactive distillation column
As mentioned previously, four operating variables were optimized. Column pressure was kept constant at 20 bar, based on the literature 5 . Several simulations were performed with different values of the parameters presented in Table 1. The objective was to maximize IPA propene-free mole fraction in the distillate, X 28 , with the constraint that distillate temperature, T D , had to be above water freezing point. Since water is always present in the product from condenser, and due to the presence of the azeotrope point, cooling the product in condenser to very low temperatures, has no benefit because the azeotrope point prevents high purification of IPA in RD column. Other researchers, such as Niu and Rangaiah 16 , used cooling water as the coolant in condenser as well.
Simulation results for 8 different cases are shown in Table 3. As may be seen, by increasing reflux ratio, T D was increased, and X 28 firstly increased up to the reflux ratio of 24 and then decreased. Thus, the best value of reflux ratio was 24. Increasing reflux ratio to a certain value will increase IPA mole fraction, because reflux stream contains high amounts of unreacted components of water and propene, which can increase the reaction rate. However, since reflux temperature is very low compared to the stage where it enters, it decreases the stage temperature by absorbing heat. This phenomenon extends to the reaction zone causing a decrease in temperature of the zone, and hence a decrease in the reaction rate. The reflux value of 24 is an optimum value at which increasing the amounts of reactants and decreasing temperature in the reaction zone balance each other.
Propene feed flow rate, F 23 was changed at the optimum value of reflux ratio and fixed values of other variables. In this case, increasing F 23 continuously decreased T D as well as X 28 . Thus, it was concluded that 36 kmol h -1 is the optimum value of this parameter. By changing water to propene feed ratio, it was observed that its best value was 1.2. A further decline of this parameter cannot be done because it is desired to have water in excess to have IPA prod- uct in distillate. At steady-state conditions for a conventional distillation column, changing feed flow rate does not affect product compositions. However, in a RD column with a reaction zone, increasing flow rates of reactants causes increased reaction rate, and thus more products are produced, in this case more IPA. In addition, since the reaction is exothermic in this case, increasing feed rate also increases temperature of the reaction zone due to higher reaction rate, and thus affects compositions of all components in RD column because of changing phase equilibrium conditions.
As previously mentioned, a partial condenser was considered to maintain T D at a desired level. This goal was achieved by considering a certain vapor fraction, β 28 , in the condenser. Decreasing its value from the base value of 0.07 resulted in a sharp decrease in T D , while its increase caused T D to in-crease and X 28 to decrease. Thus, the base value of β 28 , i.e., 0.07, was found to be optimal.
Additionally, temperature of vapor outlet from stage 27, which entered into condenser, is shown in 7th row of Table 3, which was about 443 K for the optimum case. This value was high enough to be cooled down to 278.5 K by chilled water at 275.5 K assuming the temperature difference of 3 K between hot and cold stream outlets from the condenser. However, a larger condenser may be needed to perform the cooling task.
Another tip is that the temperature of stage 27, which was near the reaction zone, was the highest one in the column, even higher than the temperature in reboiler, as reported in 8th row of the Table 3. This was because the reaction was highly exothermic and increased the temperature near the reaction zone. Fig. 6 shows the comparison of profiles of IPA propene-free contents at steady-state conditions for the base case and optimum case depicted with respect to the stages of the column. Stage 1 is the reboiler, and stage 28 is the condenser. As may be seen, profiles are very similar for both cases; IPA propene-free mole fraction in condenser for the base case was even slightly higher than in the optimum case. The reason for selecting the optimum case rather than the base case was the higher condenser temperature for the optimum case, which allows using chilled water as a coolant. IPA mole fraction increased sharply from stage 25 due to the reaction taking place in the reaction zone, stages 23 to 25, and also because of decreasing temperature in stages near the condenser.

Process control
As mentioned previously, distillate temperature and reflux ratio were selected as controlled and manipulated variables, respectively. Firstly, reflux ratio was changed from 24 to 26 at time 20 min at opti- mum conditions in order to obtain step response of distillate temperature, shown in Fig. 7. In the next step, parameters of the assumed first order model with time delay, eq. (10), were estimated by the reaction curve 21 and curve fitting methods as described previously. The results are shown in Table  4, and estimated responses with two methods are shown in Fig. 7. As it is evident, the curve fitting method gave a better fitting of the response. The results of system identification were used for the controller design. After calculating ultimate controller gain, K u , and ultimate period, p u , using frequency response method 22 , the appropriate values of controller settings were calculated using Ziegler-Nichols method 22 , as presented in Table 5.
The designed controllers were applied on RD system to control the distillate temperature and, as a result, the IPA content in the distillate. Two experiments, including set-point tracking and disturbance rejection were performed to reveal the performance of controllers. In the first experiment, set-point of distillate temperature was changed suddenly by +20 K from its optimum value after 10 minutes from the start of simulation, and the value of distillate temperature was recorded and depicted for different controllers. Fig. 8 shows the results. As shown, all controllers tracked set-point change rapidly and reached the final value in less than 5 minutes. Both PI and PID controllers had no offset, and PID controller generated smaller overshoot than PI. Even P controller exhibited good performance but with some oscillations and resulted in a very small offset of about 0.05 K. Thus, PID controller had the best performance, although, considering economic issues, P controller may be the best choice.
In the second experiment, the performance of controllers for disturbance rejection was examined. For this purpose, propene feed flow rate was changed stepwise by +10 % from its optimum value at time 10 min. As shown in Fig. 9

Conclusions
In this study, dynamic modeling and simulation of isopropanol production in a reactive distillation column were performed. Direct hydration of propylene on HZSM5 catalyst was considered, and a reaction kinetics from the literature was used for calculation of the reaction rate. Results of simulation of a base case showed that the model was able to simulate RD column with good accuracy, and IPA product was obtained from distillate with a near-toazeotrope composition. Since distillate temperature in the base case was below the freezing point of water, an optimization study was conducted to force the process to obey this constraint while maximizing IPA content in the distillate. Four operating variables, including reflux ratio, propene feed flow rate, water to propene feed ratio, and vapor fraction in condenser were selected, and their values were optimized by performing several simulations. The results of optimization showed that the reflux ratio of 24, propene feed flow rate of 36 kmol h -1 , water to propene feed ratio of 1.2, and condenser vapor fraction of 0.07 were the best values of these variables. IPA propene-free content in the distillate reached the value of 0.583 and distillate temperature was 278.5 K. Control of distillate temperature was then performed at optimum conditions by designing and applying different classic controllers. The controllers designed with Ziegler-Nichols method displayed good performance in terms of rapid response, minimum overshoot, and minimum or no offset. All previous works lack the kinetics for IPA production, but only consider chemical equilibrium. In this work, process simulation was performed by considering reaction kinetics for propene hydration.