RETROMEZCLA: A DYNAMIC STIRRED TANK REACTOR SIMULATOR Santiago Builes, Federico Calle, Carlos A. Henao, Jorge A. Velásquez Department of Chemical Engineering Universidad Pontificia Bolivariana Medellín, Colombia

Correspondence: Dr. Jorge A. Velásquez Dept. of Chemical Engineering Universidad Pontificia Bolivariana Circular 1ª #70 - 01 Medellín, Antioquia Colombia [email protected] Tel. (574)-4159020 Ext. 9598

Retromezcla, page 1

ABSTRACT Retromezcla is a computer program based on MS Excel designed to simulate the dynamic operation of stirred tank reactors. It allows to simulate batch, semi-batch and continuous reactors. Retromezcla calculates the molar fractions of the substances, the temperature and volume inside the reactor as a function of time. Retromezcla can simulate multiple reactions using first order, reversible and/or Yang and Hougen mechanisms. Retromezcla has a database to store the properties and interaction coefficients, the database can be edited and new substances can be added. Retromezcla uses Patel-Teja EOS and Wong-Sandler mixing rule to predict the thermodynamic properties of the substances crossing the boundaries of the reactor. It was evaluated using the heterogeneous esterification of acetic acid with methanol, and it provided accurate results during the whole simulation range. Keywords: Stirred tank reactors, dynamic simulation, chemical reactions, thermodynamic properties, Excel.

Retromezcla, page 2

INTRODUCTION Control and optimization of processes require a model to describe the system. Models are based on conservation laws of mass and energy, which give a set of equations that describe the behavior of the system. Simulators employ those models to give security limits, optimum security values inside the reactor, or to give a linear approximation describing the system near its operation limits [1]. They can give optimum values of the process variables inside the reactor, the feed, the exit and the cooling systems. Optimum values are usually the ones that maximize the conversion of a reactant or the selectivity of a reaction [2]. Simulation in steady state has been useful for the design and the evaluation of process equipment. Sometimes, its application is extended to time dependent operating conditions, such as feedforward systems. This extension is not possible on batch or semi-batch processes because temporal variation of operating conditions requires dynamic simulation. New generations of dynamic simulators for PCs have put in the hands of chemical engineers a tool easy to configure, fast and reliable, although usually their licenses have high prices [3]. It is possible to find a great amount of simulations done with stirred tank reactors. ModelMaker and other similar simulators are especially useful to make parametric sensitivity analysis. Such analysis gives ranges of values for which is impossible to guarantee a stable operation of the reactor [2]. Rashid and Bogle used SPEEDUP simulator to simulate a reaction with an irreversible first order mechanism, created by Perkins and Sargent (1982). They concluded that, although it is certain the greater the temperature the greater the conversion, it is recommendable to operate the reactor in conditions below the optimal in order to ease the control tasks [4]. Retromezcla, page 3

Kowar and Pagone created REACTR, which simulates a cylindrical stirred reactor with a jacket for cooling or heating. The reactor can operate in batch, continuous and semibatch; and in adiabatic, isothermal or nonisothermal conditions. One of the greatest limitations of this simulator has to do with the kinetic of the reactions. It only allows simulations of irreversible reactions of isomerization or first or second order additions with unitary coefficients [5]. The previous simulators have reduced options for two main reasons. First, they are developed for particular cases. On the other hand, their models have assumptions that make them unrealistic, such as not to consider density as a function of temperature. Retromezcla is different from previous developments for several reasons: on first instance, its model has few suppositions. It is possible to specify a great amount of reactions with simple or complex mechanisms, in addition to special arrangements of reactions. Second, the program has built-in the thermodynamic package and the algorithm for the solution of equations. These tools, added to the user-friendly interface, avoid the use of several different programs. In addition, operating as a complement of MS Excel® allows the user to manipulate the results of the simulation and easily elaborate tables and graphs for the results.

RETROMEZCLA’S MODEL DEVELOPMENT It is necessary to develop an algorithm to solve the mass and energy balances in the reactor on every step during the whole simulation range. The computer program must contain the algorithm and numerical methods to solve the resultant

Retromezcla, page 4

set of equations. Models to calculate thermodynamic properties of multicomponent mixtures and equations to calculate the reaction kinetics are also required. The mass balance for the reactor is shown on Eq. (1). NR

Fi ,in − Fi , out + Vrxn ∑ α ij ⋅ rj = j =1

dNi dt

i = 1,2..., NC

(1)

The transient energy balance is shown on Eq. (2):

Fin ⋅ ein − Fout ⋅ eout + Q + W =

dEsys

(2)

dt

Eq. (3) corresponds to the energy balance for the reactor at hand.

Fin ⋅ hin − Fout ⋅ hout + Q =

dU dt

(3)

Eq. (3) is derived from Eq. (2), according to the relationship between internal energy and enthalpy and making the following assumptions: • The shaft work of the stirrer has an insignificant value compared to other forms of energy. • The electric, magnetic, gravitational, motion and surface tension effects can be ignored (simple compressible system). • The kinetic and potential energies have insignificant values compared to other forms of energy.

Thermodynamic properties The control volume of the simulated reactor contains solely liquid phase. The program does not consider the gaseous phase that could be associated to the liquid. Thus, the selected thermodynamic package must be robust for the

Retromezcla, page 5

calculation of properties in the liquid phase. The equation of state (EOS) of Patel and Teja (1982) was selected. The chosen mixing rule was that of Wong and Sandler (1992). The molar heat capacities were determined based on the model developed by Bureš (1986). Patel Teja EOS accurately predicts the volumetric properties of liquids. For nonpolar substances its results are comparable with those from the equations of state by Peng and Robinson (1976) and the one by Soave, Redlich and Kwong (1972). For polar substances, the values are adjusted to produce more accurate results than the aforementioned equations [6]. Wong and Sandler developed a mixing rule based on the Helmholtz excess free energy. This rule has received special attention for being density independent and for allowing its extension to high-pressure systems from parameters determined in low-pressure systems by means of extrapolation [7]. Yang et al (1997) extended this mixing rule to the 3-parameters Patel and Teja equation of state. Yang et al established that their model is applicable with a good level of reliability to the prediction of the vapor-liquid equilibrium of numerous systems, including highly asymmetric systems. In addition, it can be used to the determination of some volumetric properties of liquids with reasonable accuracy [8]. The parabolic and hyperbolic equations used to evaluate the molar heat capacities as functions of the temperature do not allow a reliable extrapolation. This difficulty can be overcome by using the equation developed by Bureš. This equation uses a three parameter non-polynomial function to evaluate the heat capacities [9].

Numerical methods Retromezcla, page 6

The dynamic model of a stirred tank reactor generates a system of algebraicdifferential equations of high complexity that does not allow an analytical resolution. It is necessary to resort to numerical methods to obtain a solution. The selected method for the solution of the resulting system of differential equations is Runge-Kutta’s fourth order method (RK4). It is a stable method, and it does not require the calculation of the derivative of the functions. RK4 is a technique that offers an excellent relationship between precision and number of calculations. Brent’s method was selected for the solution of the reactor’s outlet temperature on every step. This method combines the linear interpolation and the inverse quadratic interpolation with bisection, according to the conditions of the problem, to determine the root of an equation in a given interval [10]. This method was modified in order to increase its speed. The modification is simple: the algebraicdifferential system of equations formed by the additional equations and mass and energy balances requires for its solution a set of initial conditions. Within this set, there is an initial value of the temperature inside the reactor. This value is provided to Brent’s method intending that it always use such value as the first approach to the root. As a result, the solution interval is quickly reduced. Therefore, the simulation concludes in a smaller time than with the original method.

Algorithm An algorithm that organizes the mentioned set of equations and the specifications given by the user in order to complete the simulation on every step of the simulation range was developed. The algorithm has the following steps:

Retromezcla, page 7

a) Initialize the values that will remain constant throughout the process. It is necessary to mention that the assumption that these values remain constant is inherent to this algorithm and not properly to the exposed model (these values are: Pin, Tin, Fm in, xi,in, with i=1,2,…,NC; Pout, Fm,out and Q). b) Specify the initial values for the conditions that are variable through time (these values are: xi,out, with i =1,2,…,NC, Vrxn and Tout ). c) From values in step a), calculate the dependent magnitudes of these values (these values are: Fin, Fi,in with i =1,2,…,NC, and hin), by means of the Eqs. (4) to (6).

⎛ NC ⎞ Fin = Fm,in ⎜ ∑ xi ,in ⋅ Mwi ⎟ ⎝ i =1 ⎠

(4)

Fi ,in = Fin ⋅ x i ,in

(5)

hin = f 1 (Tin , Pin , {xi ,in })

(6)

d) Based on the values from a) and b), calculate the initial values of the other variables by using Eqs. (7) to (14) (Fout(0), Fi,out(0), υout(0), hout(0), Ci,out(0), rj(0), Ni,(0), U(0)). These results are needed to solve the reactor balances.

⎛ NC ⎞ Fout = Fm,out ⎜ ∑ xi ,out ⋅ Mwi ⎟ ⎝ i =1 ⎠

(7)

Fi ,out ( 0 ) = Fout ⋅ x i ,out ( 0 )

(8)

υ out ( 0) = f 2 (Tout ( 0) , Pout , {xi ,out (0) })

(9)

hout (0) = f 3 (Tout ( 0) , Pout , {xi ,out (0) })

(10)

C i ,out ( 0 ) = ⋅

xi ,out ( 0)

υ out ( 0 )

(11)

Retromezcla, page 8

r j ( 0 ) = f 4 (Tout ( 0 ) , Pout , {Ci ,out ( 0) })

(12)

N i ( 0) = Vrxn(0) ⋅ Ci ,out (0)

(13)

U ( 0 ) = ∑ N i ( 0 ) ⋅ (hout ( 0 ) − Pout ⋅ υ out ( 0 ) )

(14)

NC

i =1

e) Find the values of Ni and U by solving the system of differential equations corresponding to the mass and energy balances. NR dNi = Fi ,in − Fi ,out + Vrxn ∑ α ij ⋅ r j dt j =1

(1)

dU = Fent ⋅ hent − Fsal ⋅ hsal + Q dt

(3)

f) Find the values of the simulation variables xi,out(1), Tout(1) and Vrxn(1) using the values obtained in step e) and assuming perfect mixing inside the reactor for the current solution step. First, calculate the compositions of the components by means of the Eq. (15). xi ,out (1) =

N i (1)

(15)

NC

∑N i =1

i (1)

The inner reactor temperature for the step is determined by using Brent’s method in equation (16). The temperature is an implicit variable in such equation. NC

[

U (1) = ∑ N i (1) ⋅ hout (1) (Tout ) − Pout ⋅ υ out (1) (Tout )

]

(16)

i =1

Once the inner reactor temperature has been determined, the molar volume of the liquid mixture is calculated. Its value is found by using the expression (17). NC

Vrxn (1) = ∑ N i (1) ⋅ υ out (1)

(17)

i =1

Retromezcla, page 9

g) The algorithm is repeated from step d), for each one of the values of time until the end of the defined interval, according to the step size fixed on a). It is possible to obtain negative values for some of the calculated variables during the execution of the algorithm. Mass, absolute temperatures, and absolute pressures are variables that cannot be negative, thus it is not coherent to obtain a value less than zero for any variable. Whenever this happens, the program must stop and the user has to review the specifications or the kinetic laws in order to detect the error. Sometimes, negative values obtained when calculating the volume of reaction are associated with the total consumption of the limiting reagent. The precedent algorithm was programmed on VBA (Visual BASIC for Applications) as a MS Excel® add-in giving origin to a computer program called Retromezcla. This simulator allows the user to simulate continuous reactors (CSTR), reactors with one inlet or one exit (Semibatch) or discontinuous reactors (Batch) under diverse conditions. In addition, the user can simulate single or multiple reactions with different kinetic laws.

RETROMEZCLA SIMULATOR This

simulator

has

three

main

parts:

“Substances”,

“Reactions”

and

“Specifications”, each of them composed by a few windows. On every window a help frame with tips and information is available. The “Substances” part is shown in Figure 1. In this part of the program, the user has to create a list of substances, importing them from the database, containing all the reagents, products and/or inerts of the simulation. It is possible to add new substances to the simulator’s database by using the button “Add new substance”. It Retromezcla, page 10

is also possible to modify either the properties of the existing substances or the binary coefficients of a selected pair of them by using the buttons “Edit Substance” or “Binary coefficients”, respectively. In the “Reactions” part, the user has to specify the data related to the stoichiometry of reactions taking place in the reactor as well as their kinetic laws. It is possible to specify multiple reversible or irreversible reactions with elemental or complex kinetics such as Yang & Hougen’s. An image of the “Reactions” part in its stoichiometry window is shown in Figure 2. In the “Specifications” part, the user has to assign the values of the input variables mentioned on steps a) and b) of the algorithm. According to the values specified for the inlet and the outlet, the reactor can be a batch, a semibatch or a continuous reactor.

Figure

3

shows

the

window

“Operating

conditions”

from

the

“Specifications” part of Retromezcla. Apart from the three sections mentioned above, Retromezcla has an additional section called “Simulation” in which it is possible to decide whether the reactor is isothermal or not. Furthermore, this section allows the user to choose the step size of the RK4 method, the size of the solution interval and their units. In addition, it shows a summary of the data entered by the user. Figure 5 shows the “Simulation” part of Retromezcla. Once the simulation is completed, Retromezcla opens a new MS Excel sheet in which two tables and three graphs appear with the simulation results. The graphs done by the simulator are molar fractions vs. time, temperature vs. time and volume vs. time.

Retromezcla, page 11

RESULTS AND DISCUSSION In order to show the performance of Retromezcla, an isothermal batch reactor was simulated. The simulated reaction was the reversible esterification of acetic acid by methanol to produce methyl acetate and water, as shown below. SAC −13 CH 3COOH + CH 3 OH ←⎯ ⎯ ⎯→ CH 3 COOCH 3 + H 2 O

(18)

This reversible reaction is catalyzed by SAC-13, a solid catalyst that consists of an acid resin supported by a porous matrix of silica. The action of this catalyst changes the conventional second order kinetics of the reaction into the heterogeneous-catalytic shown below [11]. − rA =

kCC (CM C A − CW C E / K e ) 1 + K M CM + K AC A + KW CW

(19)

Taking into account the data given by Liu et al Eq. (19) becomes Eq. (20) with -rA in kmol/m3h [11]. − rA =

0.218 ⋅ CM C A − 0.03516⋅ CW C E 1 + 0.16 ⋅ C M + 0.13 ⋅ C A + 3.11⋅ CW

(20)

The simulated reactor was an isothermal batch reactor operating at 60°C and 1 atm, during 700 minutes with an initial molar ratio of 2:1 between methanol and acetic acid. Liu et al maintained constant the reaction volume by the addition of THF (Tetrahydrofuran) [11]. To compare the results obtained with Retromezcla to those of Liu et al, acetic acid conversion was computed as an arithmetic average between conversion based in molar concentration and the one based in the number of moles. This was done to balance the effect of the reaction volume change over the calculations. Taking into account the whole simulation range, 4 points were selected form the ones reported Retromezcla, page 12

by Liu et al. These points were compared to those obtained by Retromezcla. The results are shown in Table 1 and Figure 5. Figure 5 shows that the curve predicted by Retromezcla fits the experimental data over the whole simulation range. This is supported by the results presented in Table 1. From the data in this table is possible to see that the maximum percentage difference obtained was 2.24%. Considering that the average percentage error is 1.42%, it is proper to conclude that Retromezcla accurately simulates reactions with complex mechanisms. It is required to consider the absolute difference between experimental and simulated, and not the percentage error, given that conversion is always between 0 and 1. The last column of Table 1 shows small differences among experimental and predicted values, which are in mean 3.34×10-3 points of conversion. This verifies the statement that Retromezcla gives accurate results when compared to actual data.

CONCLUSION Retromezcla is a versatile stirred tank reactor simulator because it allows the user to simulate different types of reactors in which multiple reactions may occur at a wide variety of operating conditions. Its database and the number of substances in it can be modified. Retromezcla works as a MS Excel add-in, which makes it very accessible for more industries than other commercial simulators and allowing it to export easily its results to create tables and graphs. The comparison made with experimental values reported for the esterification of acetic acid with methanol

Retromezcla, page 13

showed that Retromezcla is very accurate, even when kinetic laws are complicated, such as the reversible catalytic-heterogeneous ones.

NOTATION

Ci,out: Outlet molar concentration of component i. ein: Molar energy of the fluid entering the reactor. eout: Molar energy of the fluid exiting the reactor. Esys: Total energy of the fluid inside the reactor. Fin: Inlet molar flow. Fi,in: Inlet molar flow of the component i. Fi,out: Outlet molar flow of the component i. Fm,in: Inlet mass flow. Fm,out: Outlet mass flow. Fout: Outlet molar flow. hin: Inlet molar enthalpy. hout: Outlet molar enthalpy. Mwi: Molecular weight of component i. N: Total number of moles inside the reactor. Ni: number of moles of component inside the reactor e i en el interior del reactor. Ni: Number of moles of component i in side the reactor. NC: Number of components inside the reactor. NR: Number of reactions that occur inside the reactor. Q: Net heat flow that enters the reactor. rj: Global reaction rate of reaction j. Retromezcla, page 14

Tin: Inlet temperature. Tout: Outlet and inner temperature. U: Total internal energy of the fluid inside the reactor. Vrxn: Liquid volume inside the reactor. W: Net work that enters the reactor. xi,in: Inlet molar fraction of the component i. xi,out: Outlet molar fraction of the component i. GREEK LETTERS

αij: Stoichiometric coefficient of component i in the reaction j. υsal: Outlet molar volume.

REFERENCES [1] Gatzke E. P. & Doyle F.J. III. (1999). Multiple model approach for CSTR control. 14th IFAC World Congress, 7, 343. [2] Díaz de los ríos, Manuel. (1998). Simulación dinámica de reactores tipo tanque agitado. Ing. Quim., 30, 97. [3] Feliu, Josep A. et al. (2003). Match your process constraints using dynamic simulation. Chem. Eng. Prog, 99, 42. [4] Rashid, M. & Bogle, I.D.L. (1989) Dynamic operabiltiy analysis and simulation of a CSTR with exothermic reaction. Comput. Chem. Eng., 13, 327. [5] Kowar, Thomas R. & Pagone, Franco M. (2001). REACTR: An industrial chemical reactor dynamic simulation computer program. Org. Process Res. Dev., 5, 393.

Retromezcla, page 15

[6] Patel, Navin C. & Teja, Amyn S. (1982) A new cubic equation of state for fluids and fluid mixtures. Chem. Eng. Sci., 37, 463. [7] Wong, David S. H. & Sandler, Stanley I. (1992) A theoretically correct mixing rule for cubic equations of state. AlChE J., 38, 671. [8] Yang, T. et al. (1997). Extension of the Wong-Sandler mixing rule to the threeparameter Patel-Teja equation of state: application up to the near critical region. Chem. Eng. J., 67, 27. [9] Bureš, M. (1986). A nonlinear equation describing the molar heat capacities of gases as a function of temperature. International Chemical Engineering, 26, 160. [10] Brent, R. P. (1971) An algorithm with guaranteed convergence for finding a zero of a function. The Computer Journal, 14, 422. [11] Liu, Yijun, Lotero, Edgar & Goodwin Jr., James G. (2006). A comparison of the esterification

of

acetic

acid

with

methanol

using

heterogeneous

versus

homogeneous acid catalysis. J. Catal., 242, 278.

Retromezcla, page 16