A reactive distillation process for deep hydrodesulfurization of diesel: Multiplicity and operation aspects

Computers and Chemical Engineering 34 (2010) 196–209

Contents lists available at ScienceDirect

Computers and Chemical Engineering journal homepage: www.elsevier.com/locate/compchemeng

A reactive distillation process for deep hydrodesulfurization of diesel: Multiplicity and operation aspects J. Carlos Cárdenas-Guerra a , Teresa López-Arenas b , Ricardo Lobo-Oehmichen a , Eduardo S. Pérez-Cisneros a,∗ a b

Departamento de Ingeniería de Procesos e Hidráulica, Universidad Autónoma Metropolitana-Iztapalapa, Av. San Rafael Atlixco No. 186, C.P. 09430, México, D.F., Mexico Departamento de Procesos y Tecnología, Universidad Autónoma Metropolitana-Cuajimalpa, Artificios No. 40, C.P. 01120, México, D.F., Mexico

a r t i c l e

i n f o

Article history: Received 25 November 2008 Received in revised form 28 July 2009 Accepted 29 July 2009 Available online 7 August 2009 Keywords: Diesel deep HDS Multiple steady states Reactive distillation

a b s t r a c t A systematic study of the operating conditions and parameter sensibility under which multiple steady states may occur in a reactive distillation column (RDC) for diesel deep hydrodesulfurization is presented. The multiplicity analysis is performed through bifurcation diagrams for two feed case scenarios. The main variables that affect the steady state behavior are the reflux ratio, the hydrocarbon feed flowrate, the reboiler heat duty, and the catalyst load (liquid holdup). For the first case only a single steady state was found, while in the second case input and output multiplicities were determined. It is shown that the introduction in the feed of the recalcitrant 4,6-DMDBT, plays a crucial role in the RDC design and operation. A set of values of the bifurcation parameters for an appropriate operation of the RDC are recommended to avoid the multiplicity region and to reach the required output targets. © 2009 Elsevier Ltd. All rights reserved.

1. Introduction Conventional hydrotreating is a commercially proven refining process that passes a mixture of heated feedstock and hydrogen (H2 ) through a catalytic reactor to remove sulfur and other undesirable impurities. A review of the technologies for producing ultra-low sulfur diesel (ULSD) fuel reveals that current technologies can be modified to produce diesel with less than 10 parts per million (ppm) of sulfur. Nevertheless, only a small number of refineries currently produce diesel with sulfur in the 10 ppm range on a limited basis. The existence of the required technology does not ensure that all refineries will have that technology in place and in time to meet the new ULSD standards because these plants are characterized by a wide range of size, complexity, and quality of crude oil inputs. It is generally believed that a two-stage deep hydrodesulfurization (HDS) process will be required by most, if not all refiners, to achieve a diesel product with less than 10 ppm of sulfur. A design consistent with this technological approach has been proposed by Haldor Topsøe (DOE/EIA Report, 2001), which includes a first stage that reduces the sulfur content to around 250 ppm or lower, and a second stage that completes the reduction to less than 10 ppm. Knudsen, Cooper, and Topsøe (1999) suggested that to deep desulfurize cracked stocks, a higher operating pressure of the reactor is necessary. Pressure requirements would depend on the quality of

∗ Corresponding author. Tel.: +52 55 5804 4642; fax: +52 55 5804 4900. E-mail address: [email protected] (E.S. Pérez-Cisneros). 0098-1354/$ – see front matter © 2009 Elsevier Ltd. All rights reserved. doi:10.1016/j.compchemeng.2009.07.014

the crude oil and the setup of the individual refinery; thus, the level of pressure required for deep desulfurization is a key uncertainty in assessing the cost and availability of the technology. Based on the philosophy of previous works and designs (Sakanishi, Ando, Abe, & Mochida, 1991; Sakanishi, Ando, & Mochida, 1992; Takatsuka, Inoue, & Wada, 1997), it could be considered a two-stage deep desulfurization process where the first stage is a currently conventional HDS process that reduces the sulfur content up to 500 ppm, and a hypothetical second stage being a reactive distillation column (RDC), which should remove the sulfur from diesel to the desired specification. This RDC would operate at lower pressures and, hence, it will reduce H2 consumption, and energy requirements. Reactive distillation is considered as a highly promising process because it combines the requirements of in situ separation with reaction. This integration brings several advantages, i.e., energy and capital savings, increased reactant conversion, enhanced product selectivity, and improved heat integration. Nevertheless, the presence of chemical reactions in a distillation column leads to difficulties in the modeling (Taylor & Krishna, 2000). The increasing interest in reactive distillation has been accompanied by the development of various simulation algorithms in order to study the operation and control of this process (Al-Arfaj & Luyben, 2002; Hung et al., 2006; Kienle & Marquardt, 2003; Kumar & Kaistha, 2008; Monroy-Loperena, Pérez-Cisneros, & Álvarez-Ramírez, 2000; Wang, Wong, & Lee, 2003). Up to now, only few papers have addressed the application of reactive distillation to the diesel deep HDS. Taylor and Krishna

J. Carlos Cárdenas-Guerra et al. / Computers and Chemical Engineering 34 (2010) 196–209

(2000) discussed the possibility of applying the reactive distillation concepts to crude oil fractions HDS, while Krishna (2002) showed how this technology could be used. An analysis of the operating conditions to obtain ULSD in a conventional HDS process (Knudsen et al., 1999; van Hasselt et al., 1999) suggests that reactive distillation could be an interesting technological alternative for diesel deep HDS. In a reactive distillation process, the countercurrent flow is the natural operation mode and the internal flowrates requirement can be obtained through the catalyst packing arrangement, regulating the reflux and/or the boilup ratio, and placing properly the sulfured hydrocarbon feed. Viveros-García, OchoaTapia, Lobo-Oehmichen, de los Reyes-Heredia, and Pérez-Cisneros (2005) showed a comparison between deep HDS in a conventional countercurrent trickle-bed reactor and the operational and design alternatives offered by a reactive distillation process. They noted that all the operation requirements for a deep HDS in the reactor could be fulfilled by a reactive distillation operation with an appropriate process design. Design, operability, and control of a RDC become more difficult due to the interactions between chemical reaction and separation. Such complex interactions, primarily between kinetic expressions and thermodynamic models, lead to a highly nonlinear behavior of the RDC indicating the possible existence of multiple steady states (MSS). Different authors have identified MSS (input and output multiplicities) in conventional, azeotropic, and reactive distillation systems (Baur, Taylor, & Krishna, 2003; Ciric & Miao, 1994; Gani & Jørgensen, 1994; Jacobsen & Skogestad, 1991; Kannan, Joshi, Reddy, & Shah, 2005; Kienle & Marquardt, 2003; Mohl et al., 1997, 1999; Mohl, Kienle, Sundmacher, & Gilles, 2001; Müller & Marquardt, 1997; Singh, Singh, Kumar, & Naistha, 2005; Wang, Chang, Wang, & Li, 2008; Yang, Wu, Zhao, Wang, & Lu, 2006). Some of them concluded that MSS can arise in binary distillation if the flowrates are given on a mass basis instead of a molar basis. Gani and Jørgensen (1994) reported MSS in heterogeneous azeotropic distillation columns through the study and simulation of the separation of ethanol–water–benzene ternary system. They pointed out that the movement of the internal composition front could be the cause of the possible connection between multiplicity and operability. They realized that the system reached a different steady state as fronts corresponding to temperature and composition profiles move up or down depending on small positive or negative pulse disturbances in the reboiler heat duty. Müller and Marquardt (1997) verified the experimental existence of MSS for first time in a heterogeneous azeotropic distillation by rigorous simulations and bifurcation diagrams. Jacobsen and Skogestad (1994) studied the stability of conventional (non-reactive) distillation columns and showed that columns displaying MSS had at least one solution that corresponded to an unstable point and that instability was more likely with large internal flowrates. The RADFRAC module (from ASPEN PLUS process simulator) has been used by Jacobs and Krishna (1993) to study MSS in a RDC for MTBE synthesis. Their results showed high and low conversion steady states when methanol is fed to stages 10 or 11 using the Jacobs–Krishna column configuration. Physical explanation for the occurrence of MSS in the MTBE process was provided by Hauan, Hertzberg, and Lien (1995, 1997). Mohl et al. (1999) employed a pilot-scale column to produce MTBE and TAME. MSS were found experimentally when the column was used to produce TAME, but not in the MTBE process. Baur et al. (2003) presented a bifurcation analysis for the synthesis of TAME in a RDC with two different methods for describing the reaction kinetics: pseudo-homogeneous models and heterogeneous models. Both reaction models showed the possibility of MSS. Therefore, it is evident that the analysis of existence of a single steady state (SSS) or MSS should give insights into the reactive distillation process, help to avoid unsafe operating conditions, and further facilitate subsequent studies such as control, monitoring, data reconcilia-

197

tion, parameter estimation, and optimization of existing reactive distillation process. In academia and industry the nonlinear analysis of chemical engineering models has great significance (Krasnyk, Ginkel, Mangold, & Kienle, 2007). In particular, the steady state and dynamic simulations of chemical processes can be done through two computational approaches: oriented-equation models and modular models. Flowsheet simulators for chemical process, based on modular models (such as ASPEN PLUS, HYSYS, PRO/II, and CHEMCAD), have the main advantage of including: property databases, thermodynamic models, process unit models, numerical methods, etc. However, a drawback of them is that they do not provide algorithms for nonlinear analysis explicitly (i.e., bifurcation algorithms and stability analysis). While working with equation-based models, it is possible to use continuation algorithms included in packages such as AUTO, CONT, HOMPACK, PITCON, BIFPACK, DIVA (see Mangold, Kienle, Gilles, & Mohl, 2000 and references in) or ProMoT (Krasnyk et al., 2007). Specifically DIVA has been used to study MSS in a RDC of MTBE and TAME (Katariya, Kamath, Moudgalya, & Mahajani, 2008; Mangold et al., 2000). Some reported benefits of DIVA is that it includes a model library with models of standard operation units and new models are written in symbolic form, also it comes with an inbuilt package for continuation and stability analysis for DAEs systems. On the other hand, notwithstanding the modular simulators do not allow detailed nonlinear analysis, the aim of several works has been to develop procedures for the construction of bifurcation diagrams in modular simulators (Kannan et al., 2005; Vadapalli & Seader, 2001; Yang et al., 2006). Vadapalli and Seader (2001) proposed two continuation algorithms that can be added into ASPEN PLUS to calculate steady state branches in regions of multiplicity. The natural arclength algorithm for an adiabatic continuous stirred tank reactor problem and pseudo-arclength algorithm for a homogeneous azeotropic distillation problem were successful in detecting both turning points and all three steady state branches. However, the limitations were the impossibility of determining whether the steady states were stable or unstable, and of detecting exact bifurcation points because the Jacobian matrix generated by ASPEN PLUS is not accessible by the user. Kannan et al. (2005) identified MSS in homogeneous azeotropic distillation using the HYSYS process simulator in steady state mode, while the stability of the different branches was determined by dynamic simulations also in HYSYS. Nevertheless, most of the reported works have been for conventional (non-reactive) distillation columns or chemical reactors, but for reactive distillation columns (RDCs) few works have considered the construction of bifurcation diagrams using modular simulators. Yang et al. (2006) found input and output multiplicities in a RDC for synthesis of ethylene glycol, defining a procedure in ASPEN PLUS with the method of sensitivity analysis, but without establishing the stability of the branches. Thus, in this paper we present a systematic study of the operating conditions and parameter sensibility under which MSS may occur, and to assess their influence on a RDC for the deep HDS of diesel. A procedure for multiplicity analysis in ASPEN PLUS based on bifurcation diagrams is presented. The multiplicity analysis is performed for two feed case scenarios, showing its implication on the design and operation of the RDC. Specifically, we study the effect of the operating conditions and parameter sensibility over the main variables to monitor or control hereinafter: the recalcitrant organosulfur compounds conversion and the product purity.

2. The HDS reactive system There are two possible reaction pathways for sulfur removal from the organo-sulfur compounds as illustrated for dibenzothio-

198

J. Carlos Cárdenas-Guerra et al. / Computers and Chemical Engineering 34 (2010) 196–209

Fig. 1. Reaction pathways for dibenzothiophene hydrodesulfurization.

phene (DBT) and 4,6-dimethyldibenzothiophene (4,6-DMDBT) in Figs. 1 and 2, respectively. The first pathway is the sulfur atom direct extraction (hydrogenolysis) from the sulfured molecule. The second pathway is the hydrogenation of one aromatic ring followed by the sulfur atom extraction. Van Parijs and Froment (1986) suggested a reaction network in which thiophene (Th) is hydrodesulfurized to give hydrogen sulfide (H2 S), 1-butene, cis-

and trans-2-butene, followed by the secondary hydrogenation of the butenes into butane. For benzothiophene (BT), Van Parijs, Hosten, and Froment (1986) proposed a parallel pathway similar to Th, indicating also that such reactions occur at different sites. In fact, to obtain an ULSD, the most recalcitrant sulfur compounds must be eliminated from the hydrocarbon feed, this is, DBT and 4,6-DMDBT.

Fig. 2. Reaction pathways for 4,6-dimethyldibenzothiophene hydrodesulfurization.

J. Carlos Cárdenas-Guerra et al. / Computers and Chemical Engineering 34 (2010) 196–209

The liquid-phase HDS reactions for DBT and 4,6-DMDBT have been thoroughly studied and revisited by several authors (Broderick & Gates, 1981; Froment, Depauw, & Vanrysselberghe, 1994; Houalla, Nag, Sapre, Broderick, & Gates, 1978; Vanrysselberghe & Froment, 1996; Vanrysselberghe, Gall, & Froment, 1998). The reaction network for DBT hydrodesulfurization proposed by Houalla et al. (1978) is shown in Fig. 1, and indicates two primary pathways: hydrogenation of one aromatic ring gives an equilibrium mixture of 1,2,3,4-tetrahydrodibenzothiophene and 1,2,3,4,10,11-hexahydrodibenzothiophene; HDS of these two intermediates occurs rapidly to give cyclohexylbenzene (CHB). This combination of reactions is designated as the hydrogenation pathway. The second reaction route, direct HDS of DBT to give biphenyl (BPH) is referred as the hydrogenolysis pathway. BPH reacts further with H2 to give CHB, but this reaction is typically two orders of magnitude slower than the rate of DBT hydrogenolysis (Broderick & Gates, 1981). For the 4,6-DMDBT hydrodesulfurization, Vanrysselberghe et al. (1998) have proposed the reaction network shown in Fig. 2. In the hydrogenation pathway, the 4,6-DMDBT is hydrogenated into partially hydrogenated dimethyldibenzothiophenes prior to sulfur removal. These intermediates are highly reactive and they are directly desulfurized into 3-methylcyclohexyltoluene (3-MCHT) and H2 S. The hydrogenolysis reaction route gives 3,3 dimethylbiphenyl (3,3 -DMBPH) and H2 S. 3,3 -DMBPH is then hydrogenated into 3-MCHT. Experiments by Vanrysselberghe et al. (1998) have demonstrated that 4,6-DMDBT is mainly converted through the hydrogenation reaction pathway, while DBT is mostly converted by the hydrogenolysis reaction pathway. In general, when alkyl substituents are attached to the carbon atoms adjacent to the sulfur atom, the rate for direct sulfur extraction is diminished whereas the sulfur removal rate via the hydrogenation pathway is relatively unaffected.

3. Conceptual design of the reactive distillation column Viveros-García et al. (2005) developed a conceptual design of a RDC for deep HDS of diesel through a thermodynamic analysis considering the following aspects: (i) the volatility of the organosulfur compounds; (ii) the different reactivities of the organo-sulfur compounds; and (iii) the computation of non-reactive and reactive residue curve maps for DBT elimination. An equilibrium stage model considering homogeneous liquid chemical reactions was used. The validation of such model has previously been discussed (Granados-Aguilar, Viveros-García, & Pérez-Cisneros, 2008; PérezCisneros, Granados-Aguilar, Huitzil-Meléndez, & Viveros-García, 2002; Viveros-García et al., 2005). Recently, Granados-Aguilar et al. (2008) shown that the equilibrium model predicts consistently experimental values (Broderick & Gates, 1981) when only the hydrogenolysis reaction pathway of DBT is considered. The nonidealities of the vapor and liquid phases are calculated through the Peng–Robinson equation of state (Peng & Robinson, 1976) with the appropriate binary interaction parameters. The missing critical properties and acentric factors for DBT and 4,6-DMDBT as well as its respective hydrogenolysis and hydrogenation reaction products, were estimated using the Joback and Reid method (Joback & Reid, 1987). The complete set of physical properties for vapor liquid equilibrium calculations is given in Table 1. The RDC configuration obtained consisted of 14 stages with two reactive zones and three non-reactive zones (Fig. 3). Target conversion of 99% for the DBT and 4,6-DMDBT and complete elimination of Th and BT were assumed as part of the design specifications. Table 2 shows the column configuration details used for the simulations. Granados-Aguilar et al. (2008) also found that an operating pres-

199

Table 1 Physical properties required for VLE calculations from ASPEN PLUS. Compound

Tb (◦ C)

Tc (◦ C)

Pc (atm)

ω

Th BT DBT 4,6-DMDBT H2 BD ET BPH 3,3 -DMBPH H2 S n-C11 n-C12 n-C13 n-C14 n-C16

84.1 219.9 331.4 282.7 −252.7 −4.4 136.2 255 231.5 −60.3 195.9 216.3 235.4 253.5 286.8

306.2 480.8 623.8 530.7 −239.9 151.8 344 499.8 446.5 100.3 365.8 384.8 401.8 419.8 449.8

56.1 40.8 38 39.6 12.9 42.6 35.6 33.3 31.3 88.4 19.2 17.9 16.5 15.4 13.8

0.1969 0.2954 0.3982 0.5322 −0.2159 0.1950 0.3034 0.4028 0.5060 0.0941 0.5303 0.5763 0.6173 0.6430 0.7174

Tb is the boiling point, Tc is the critical temperature, Pc is the critical pressure, and ω is the acentric factor.

sure of 30 atm and a H2 to hydrocarbon (HC) feed ratio of 3 were the optimal values in a simple reactive batch distillation process to reach deep HDS of diesel. 4. Multiplicity analysis RDCs mathematical models are highly nonlinear, and MSS solutions have been reported by many researchers, i.e., Baur et al. (2003), Chen, Huss, Doherty, and Malone (2002), Katariya et al. (2008), Kienle and Marquardt (2003), Mangold et al. (2000), Mohl et al. (1999), Wang et al. (2008), Yang et al. (2006). However, none of these works has addressed the HDS process; most of them studied the MTBE and TAME cases. In the present work, the MSS solutions are analyzed through bifurcation diagrams, which are built using a continuation method (Guckenheimer & Holmes, 1983), raking the steady state of the model when a bifurcation parameter value is increased or decreased. Two multiplicity types can be found: input and output multiplicities, but a combined input–output multiplicity may also be present (see Fig. 4). Output multiplicity occurs when

Fig. 3. Reactive distillation column configuration for ultra-low sulfur diesel production.

200

J. Carlos Cárdenas-Guerra et al. / Computers and Chemical Engineering 34 (2010) 196–209

Table 2 Design specifications of the reactive distillation column. Specification

Value

Specification

Number of total stages Stages of reactive zone I Stages of reactive zone II HC feed stage H2 feed stage HC feed flowrate (kmol/h)

14 5–7 10–12 9 12 100

H2 /HC feed ratio Distillate flowrate (kmol/h) Reflux ratio Column pressure (atm) Partial condenser temperature (C) Holdup (kg of catalyst)

one set of input variables (manipulated variables) results in two or more unique and independent sets of output variables (measured variables). In chemical reactors and conventional (non-reactive) distillation columns there are usually three steady states associated to ignition (high conversion), extinction (low conversion) and medium unstable conversion, as shown in Fig. 4(a). On the other hand, input multiplicity occurs when two or more unique sets of input variables produce the same output condition (Fig. 4(b)). This type of multiplicity has important implications for close loop control since it is related to the so-called zero dynamics of the system, which is associated with unusual, unexpected, or inverse responses of the outputs after a step change, has been applied to the inputs. As RDCs share some common features with chemical reactors and conventional (non-reactive) distillation columns, the RDC behavior may exhibit input–output multiplicity with three steady states as shown in Fig. 4(c), or output multiplicity, with a large number of different steady states induced by ignition and extinction of every single reactive column tray (Kienle & Marquardt, 2003). For typical bifurcation diagrams with three MSS (as shown in Fig. 4), the stable steady states are denoted by a continuous line and the unstable ones by a dashed line. The latter corresponds to the steady states associated to medium conversion in an output multiplicity. 4.1. Basic procedure for tracking multiple steady states with ASPEN PLUS In general, the construction of bifurcation diagrams using modular simulators is not straightforward; because most applied

Fig. 4. Steady state multiplicities: (a) output multiplicity; (b) input multiplicity; (c) input–output multiplicity; and (d) basic procedure for multiplicity analysis.

Value 3 340 0.5 30 225 10,000

numerical algorithms for nonlinear equation systems, as incorporated into steady state process simulators seek, at best, one solution. So that additional steady state solutions may be found by defining case studies, and performing a sensitivity analysis on one or more parameters, or varying the initial guesses of the previous solution. It should be understood that, in this case, by sensitivity analysis, we mean the study of the effect of changes in input variables on process (model) outputs. Specifically, it means to perform in ASPEN PLUS the following four steps: (i) specification of measured (sampled) variable(s); (ii) specification of manipulated (varied) variable(s); (iii) specification of the range(s) for manipulated (varied) variable(s); and (iv) specification of the quantities to calculate and tabulate. Some works have been focused to develop systematic algorithms to do this through ASPEN PLUS (Vadapalli & Seader, 2001; Yang et al., 2006) and HYSYS (Kannan et al., 2005). In this paper, the goal is not to develop a new generalized procedure but to establish a basic methodology for tracking MSS of the RDC for the HDS process in ASPEN PLUS. First, the RDC was simulated using the RADFRAC module. For this purpose, ASPEN PLUS required the specification of components, property model, feed conditions, operating pressure, column configuration (number of stages, feed stage location, reactive zone designation, holdups, type of condenser and reboiler), two operating specifications (in our case, reflux ratio and distillate flowrate), and the reaction type. In the present study the reaction type was defined as a kinetically controlled system and the reaction rate expressions (given in Appendix A) were incorporated by a usersupplied Fortran subroutine. Afterwards, similar to the procedures proposed by Vadapalli and Seader (2001) and Yang et al. (2006), a sensitivity analysis was carried out to assess the influence of the operating parameters on the RDC performance. However, the methodology proposed here differs mainly in the tracking of zones close to the turning points and of the middle branch. The sensitivity option in ASPEN PLUS was selected so that current results were used as initial estimates for the next steady state solution. In this case, only one operating parameter and the initial estimates were variable. Having in mind a typical bifurcation curve as shown in Fig. 4(d), the procedure was divided in three steps as follows: 4.1.1. Step 1: Upper branch Starting from the operating conditions of the reference steady state (point A in Fig. 4(d)), which is located out of the MSS region, one operating parameter was selected to be varied. Then, using the sensitivity analysis, the simulations were trivially converged to extreme (very high) values of the parameter at first (point B). Once again, starting from the initial point A, the parameter was varied continuously with a small step towards the other extreme (lower values) trying to reach the first turning point (TP1). However, as the region of multiplicity passes, the upper branch no longer exits. The program fails to converge and suddenly jumps to the lower branch (point C). To track TP1 as close as possible, small steps (equivalent to a large number of solution points) should be used. In fact, 150 solution points were evaluated from point A to B and other 150 points from point A to TP1. Moreover, to refine the value of the TP1, the simulation was reinitialized at each step with current solution as

J. Carlos Cárdenas-Guerra et al. / Computers and Chemical Engineering 34 (2010) 196–209

initial estimates, and then the step size was reduced by half until the step size was small enough and until divergence occurred (meaning that the TP1 was closely evaluated since its Jacobian becomes infinity). 4.1.2. Step 2: Lower branch The described procedure in step 1 was repeated but starting from point C. First, the operating parameter was decreased with a small step towards extreme (very low) values, reaching point D. Then, the parameter was increased towards the second turning point (TP2). Again, 150 points were evaluated for each one of these two regions (from point C to D and from point C to TP2). Finally, the value of the TP2 was refined. 4.1.3. Step 3: Middle branch The presence of upper and lower branches, or equivalent to the presence of two turning points, assesses the existence of an intermediate branch (commonly unstable). For tracking this branch initial guesses were required to converge to one point (point E). For this purpose, a bisection algorithm was applied. Taking all results (i.e., vapor and liquid compositions, vapor and liquid flowrates, and temperatures at each stage of the RDC) from the generated solution for TP1 and TP2, new initial guesses were computed as the average of them. The simulation was reinitialized with the operating parameter set at the average of values of TP1 and TP2. If the simulation converged to a new steady state (point E), the procedure of steps 1 and 2 was followed for tracking the complete middle branch (in this case, the evaluation of 75 points from point E to TP1 and other 75 points from point E to TP2). Otherwise, if the simulation converged to the higher or lower branch, a new value (slightly higher or lower) of the operating parameter was set around the average of values of TP1 and TP2. This iterative algorithm was applied until convergence at the point E is achieved. Generally, convergence was obtained after few iterations depending on the bifurcation parameter. It is worth of mention that the total number of evaluated points for the bifurcation diagrams (shown in Figs. 7, 8, 11, 12, 13 and 14) was 750 points, which guarantees that the increased/decreased step for the sensitivity analysis was small enough so that no bifurcation branches were skip or avoid. Nevertheless, in these figures not all of the evaluating points are shown in order to have a better visualization of the branches. This basic simulation procedure worked favorably for the RADFRAC module of ASPEN PLUS in order to build the bifurcation diagrams when reflux ratio, reboiler heat duty, and holdup were used as bifurcation parameters. But, when hydrocarbon feed flowrate is used as bifurcation parameter, a convergence problem using the sensitivity analysis tool of ASPEN PLUS was found. So, in this case the methodology was followed as abovementioned but the increment and decrement of the operating parameter was manually done. Afterwards, all the solution

201

curves or bifurcation diagrams were analyzed and the multiplicity existence in the operating region was determined. The methodology described also has some limitations as the ones reported by some authors (Dalal & Malik, 2003; Kannan et al., 2005; Vadapalli & Seader, 2001; Yang et al., 2006) when (modular) steady state simulators are used. For instance, the Jacobian matrix generated in the calculations is not accessible by the user and the stability of the solution cannot be verified straightforward. Only equation-oriented simulators can explicitly show this information. Otherwise the stability of the branches could be verified through modular dynamic simulators. It has been widely reported that the S-shaped bifurcation curve of three steady state solution branches can have only two possible stability characteristics (Vadapalli & Seader, 2001): the three solution branches are either stable–unstable–stable or unstable–stable–unstable. Nonetheless, most of the reported case studies in chemical processes correspond to the former stability case, including conventional (non-reactive) and reactive distillation columns. 4.2. Case studies With the purpose of finding SSS or MSS regions for the RDC, two feed composition cases considering different HDS reaction pathways are studied. The cases are as follows: Case 1: HC feed composition with Th, BT, and DBT as the organosulfur compounds present in the HC feed (given in Table 3) and following only the hydrogenolysis reaction pathway (see Fig. 1). Case 2: HC feed composition considering Th, BT, DBT, and 4,6DMDBT as the organo-sulfur compounds present in the HC feed (given in Table 3) and following both, the hydrogenolysis and hydrogenation reaction pathways (see Figs. 1 and 2). Two reference steady states are defined to achieve a design target conversion of 99.99% for DBT and 99.3% for 4,6-DMDBT considering the feed composition and the reaction route described in Table 3. Afterwards, bifurcation diagrams are constructed to find the SSS and MSS regions by increasing or decreasing several operating or design parameters such as the reflux ratio, hydrocarbon feed flowrate, reboiler heat duty, and catalyst load (liquid holdup). In particular, to show the effect of liquid holdup, it is assumed that the catalyst in the reactive zones is completely wet. The HDS reaction rate expressions, reaction rate coefficients, and adsorption equilibrium constants used for the determination of the steady state solutions are given in Appendix A. 4.3. Reference steady states Figs. 5 and 6 show the liquid composition profiles for H2 S and the organo-sulfur compounds (Th, BT, DBT, and 4,6-DMDBT), the

Table 3 Hydrocarbon feed composition. Reaction pathway Component

Hydrogenolysis Feed composition for case 1 (mole fraction)

Hydrogenolysis and hydrogenation Feed composition for case 2 (mole fraction)

Th BT DBT 4,6-DMDBT n-C11 n-C12 n-C13 n-C14 n-C16

0.0087 0.0087 0.1000 0.0000 0.4966 0.3166 0.0089 0.0015 0.0589

0.0080 0.0080 0.1000 0.0200 0.4890 0.3160 0.0080 0.0010 0.0500

Total sulfur content (wt%)

2.25

2.60

n-C11 : n-Undecane; n-C12 : n-dodecane; n-C13 : n-tridecane; n-C14 : n-tetradecane; n-C16 : n-hexadecane.

202

J. Carlos Cárdenas-Guerra et al. / Computers and Chemical Engineering 34 (2010) 196–209

Fig. 5. Reactive distillation column profiles for the reference steady state (case 1): (a) liquid composition; (b) generation rates; and (c) temperature.

Fig. 6. Reactive distillation column profiles for the reference steady state (case 2): (a) liquid composition; (b) generation rates; and (c) temperature.

Table 4 HDS reactive distillation column simulation results for the reference steady states. Parameter

SSS-1

Distillate flowrate (kmol/h)

Vapor distillate

Liquid distillate

Composition (mole fraction)

Distillate flowrate (kmol/h) Composition (mole fraction)

yH2 = 0.8933 yH2 S = 0.0374 yTh = 1.7631E−33 yBT = 7.2130E−07 yDBT = 3.2859E−09 yBD = 0.0026 yBPH = 0.0009

xH2 = 0.0408 xH2 S = 0.0097 xTh = 2.7986E−33 xBT = 9.7036E−06 xDBT = 3.0022E−07 xBD = 0.0014 xBPH = 0.0278

Bottoms flowrate (kmol/h)

Bottoms

Composition (mole fraction)

Reflux ratio Condenser heat duty (kW) Bottoms temperature (◦ C) Reboiler heat duty (kW) DBT conversion design target (%) 4,6-DMDBT conversion design target (%) Total sulfur content in ULSD (wt%) BD: 1,3-Butadiene; ET: ethylbenzene.

xH2 = 0.0663 xH2 S = 0.0006 xTh = 3.2149E−20 xBT = 9.2188E−06 xDBT = 1.4438E−05 xBD = 4.0844E−05 xBPH = 0.1832

SSS-2 303.3 yn-C11 = 0.0484 yn-C12 = 0.0150 yn-C13 = 0.0001 yn-C14 = 4.3185E − 06 yn-C16 = 5.4725E − 06 yET = 0.0018

36.6 xn-C11 = 0.6218 xn-C12 = 0.2863 xn-C13 = 0.0037 xn-C14 = 0.0001 xn-C16 = 0.0005 xET = 0.0074

47.3 xn-C11 = 0.2567 xn-C12 = 0.3501 xn-C13 = 0.0150 xn-C14 = 0.0030 xn-C16 = 0.1239 xET = 0.0007

0.5 4630.6 414.6 5662.1 99.99 – 0.0005

yH2 = 0.8881 yH2 S = 0.0434 yTh = 0.0000 yBT = 6.2932E−07 yDBT = 7.4351E−10 y4,6-DMDBT = 3.3341E−07 yET = 0.0016 y3,3 -DMBPH = 0.0008

301.3 yn-C11 = 0.0471 yn-C12 = 0.0151 yn-C13 = 0.0001 yn-C14 = 2.9451E − 06 yn-C16 = 4.8680E − 06 yBD = 0.0024 yBPH = 0.0009 38.6

xH2 = 0.0401 xH2 S = 0.0112 xTh = 0.0000 xBT = 8.5223E−06 xDBT = 6.8807E−08 x4,6-DMDBT = 1.2790E−05 xET = 0.0068 x3,3 -DMBPH = 0.0168

xn-C11 = 0.6043 xn-C12 = 0.2877 xn-C13 = 0.0034 xn-C14 = 0.0001 xn-C16 = 0.0004 xBD = 0.0013 xBPH = 0.0274 45.6

xH2 = 0.0634 xH2 S = 0.0010 xTh = 0.0000 xBT = 8.4708E−06 xDBT = 0.0002 x4,6-DMDBT = 0.0002 xET = 0.0006 x3,3 -DMBPH = 0.0237

xn-C11 = 0.2476 xn-C12 = 0.3483 xn-C13 = 0.0138 xn-C14 = 0.0020 xn-C16 = 0.1091 xBD = 3.5268E−05 xBPH = 0.1894 0.5 4698.1 414.7 5562.4 99.9 99.3 0.0099

J. Carlos Cárdenas-Guerra et al. / Computers and Chemical Engineering 34 (2010) 196–209

203

consumption profiles of the organo-sulfur compounds, and the temperature profile along the RDC for the two feed cases, respectively. It may be noted in Figs. 5(a) and 6(a) that the maximum component compositions are located at the HC feed stage (stage 9). The high concentration of DBT at this point is due to the amount of this compound in the feed stream and its high boiling temperature. It is also observed that as the organo-sulfur compounds go into the reactive zones (stages 5–7 and 10–12), the concentrations of the three or four (according to the feed case) species fall sharply to zero. The liquid composition profile suggests that the H2 S is being concentrated mainly at the rectifying section. Therefore, the catalyst activity inhibition due to the presence of H2 S may be reduced at the stripping section leading to a higher conversion of the heaviest organo-sulfur compounds at this zone. Figs. 5(b) and 6(b) show that the sulfur species are mainly consumed at stages 7 and 10, which are the incoming stages to the reactive zones. It should be noticed that between the reactive zones there are two non-reactive stages (8 and 9) performing the separation of the paraffin and organo-sulfur compounds. The temperature profiles for both feed cases (Figs. 5(c) and 6(c)) show that reactive zone I operates close to isothermal conditions, while at the reactive zone II the temperature increases around 15 ◦ C. The isothermal behavior of the reactive zone I can be explained considering that the mixture incoming to this zone is a regular paraffin solvent mixture with little amount of the heaviest organo-sulfur compounds (see Figs. 5(a) and 6(a)). Table 4 shows the simulation results for both reference single steady states (SSS-1 for case 1 and SSS-2 for case 2). For case 2, it can be noted that the H2 S concentration increases and the H2 concentration decreases at the exit streams. It should pointed out that in both feed cases there is a significant reduction of the sulfur content from 2.25 to 0.0005 wt% (weight percentage) for case 1 and from 2.60 to 0.0099 wt% for case 2. In general, it may be said that the RDC design proposed allows the elimination of the organo-sulfur compounds for both feed cases and such sulfur reduction is achieved under similar operating conditions (heat duties of condenser and reboiler, bottoms temperature, distillate and bottoms flowrates). 4.4. Bifurcation diagrams for case 1 Bifurcation diagrams were built considering the reflux ratio (RR), the hydrocarbon feed flowrate (HCFF), the reboiler heat duty (QR), and the liquid holdup (HLP) or amount of catalyst loaded as bifurcation parameters. These parameters were selected since they are key variables (i.e., manipulated or disturbance variables) to operate the process. The conversion of the most recalcitrant organo-sulfur compounds (DBT and 4,6-DMDBT) are the output variables determining the steady state multiplicity existence. Fig. 7(a–d) shows the influence of such variables for case 1, where the nonlinearities of 4,6-DMDBT kinetic expressions are not incorporated in the computation. The vertical tie line on these figures represents the reference single steady state (SSS-1) for case 1. It can be observed in Fig. 7(a–d) that all bifurcation diagrams exhibited no evidence of multiplicity (i.e., only SSS). Fig. 7(a) illustrates the conversion of DBT with the variation of the RR and it shows that with a RR of 0.15 or greater, a complete elimination of DBT is achieved. However, an excessive reflux leads to an increase of the “cold” energy required to condensate the stream returning to the column, especially when the highly volatile H2 and H2 S are present. Fig. 7(b) depicts the conversion of DBT with the variation of the HCFF and, from 350 kmol/h a considerable fall of the DBT conversion occurs. This is an expected behavior since the amount of DBT fed into the RDC is augmented as the HCFF increases, while the amount of H2 fed remains constant, this is, the H2 /HC feed ratio is reduced from its optimal value of 3. Fig. 7(c) shows the conversion of DBT with the variation of the QR. It can bee

Fig. 7. Bifurcation diagrams for DBT conversion (case 1) with the variation of: (a) reflux ratio; (b) hydrocarbon feed flowrate; (c) reboiler heat duty; and (d) holdup.

seen at Fig. 7(c) that the amount of energy used in the SSS-1 may be reduced to a value of around 4098 kW and still the complete elimination of DBT is achieve. Notice also that, a heat reduction from 2500 kW would lead to a sharp decrease of the DBT conversion. Fig. 7(d) shows the conversion of DBT with the variation of the catalyst loaded at the reactive zones. It can be noted that the DBT conversion design target (99.99%) is only achieved for amounts greater than 5000 kg of catalyst loaded. 4.5. Bifurcation diagrams for case 2 4.5.1. Influence of the reflux ratio The influence of the reflux ratio (RR) as bifurcation parameter considering the feed composition for case 2 is shown in Figs. 8–10. In this case, the conversions of DBT and 4,6-DMDBT are the output variables determining the steady state multiplicity existence. It can be noted in Fig. 8 that the MSS region is found between the RR values of 0.22 and 0.345. The stable steady states with high conversion of DBT and 4,6-DMDBT are located in RR values greater than 0.345, outside of the MSS region. It is interesting to observe that,

204

J. Carlos Cárdenas-Guerra et al. / Computers and Chemical Engineering 34 (2010) 196–209

Fig. 8. Bifurcation diagram for DBT and 4,6-DMDBT conversion (case 2) with the variation of the reflux ratio.

while the DBT conversion is highly affected by the variations of the reflux ratio from 0.1 to 0.345, the 4,6-DMDBT conversion is practically constant (∼99%). In fact, at this reflux ratio range, it is found that there are three steady states for DBT conversion (low, medium, and high) and an output multiplicity region is determined for DBT conversion. However, if a zoom of the high 4,6-DMDBT conversion is performed, an input–output multiplicity region is found. In this region, there exist three possible input solutions (for the reflux

Fig. 10. Reactive distillation column profiles for case 2 with a reflux ratio of 0.3: (a) H2 liquid composition; (b) n-C16 liquid composition; (c) n-C12 liquid composition; and (d) n-C11 liquid composition.

Fig. 9. Reactive distillation column profiles for case 2 with a reflux ratio of 0.3: (a) DBT liquid composition; (b) H2 S liquid composition; (c) 4,6-DMDBT liquid composition; and (d) temperature.

ratio) to get the same output design target (i.e., 99.3% conversion of 4,6-DMDBT) approximately at RR values of 0.257, 0.317, and 0.36. Moreover, in the same region, for 4,6-DMDBT conversions between 99% and lightly less of 99.3%, there exist four possible input solutions which implies to keep strict operating conditions and a robust controller. Such type of multiplicity is relevant since it shows the strong influence of the nonlinearity of the kinetic model for the 4,6-DMDBT elimination. Fig. 9 shows the composition profiles for DBT, H2 S, and 4,6DMDBT as well the temperature profile along the RDC for the three steady states evaluated at RR of 0.3. It can be noted that while the DBT shows three liquid composition profiles (see Fig. 9(a)) corresponding to high, medium, and low conversion, the 4,6-DMDBT liquid composition profile and the temperature profile are not significantly affected (see Fig. 9(c) and (d)). As seen in Fig. 9(b) there are three H2 S liquid composition profiles corresponding to a low, medium, and high conversion of the organo-sulfur compounds. It should be noticed that for a low conversion of DBT, this is a higher liquid composition profile of DBT (see Fig. 9(a)), corresponds to a lower H2 S liquid composition profile in Fig. 9(b). Fig. 10 shows the liquid composition profiles for H2 and some paraffin (n-C11 , n-C12 and n-C16 ) along the RDC for the three different steady states evaluated at RR = 0.3. It can be observed in Fig. 10(a) that the H2 liquid composition has a maximum value at the exit of the reactive zone I (stage 5) and it decreases sharply at the top of the column. On the other side, at the bottom of the RDC, a minimum value for the H2 liquid composition is obtained for the low conversion steady state. It should be pointed out that the low conversion steady state

J. Carlos Cárdenas-Guerra et al. / Computers and Chemical Engineering 34 (2010) 196–209

Fig. 11. Bifurcation diagram for DBT and 4,6-DMDBT conversion (case 2) with the variation of the hydrocarbon feed flowrate.

in Fig. 10(a) corresponds to the lowest temperature profile for the RDC showed in Fig. 9(d). That is, Fig. 9(d) shows that there exist a temperature difference, between the high conversion and low conversion steady states profiles of around 5–15◦ , and even, with this small differences, the H2 solubility and the DBT conversion are strongly affected, as shown in Figs. 9(a) and 10(a), respectively. The above may be explained considering that, a lower temperature promotes an increment in the H2 solubility (this is, a higher H2 liquid composition along the RDC), but it reduces the reaction rate for DBT elimination. Fig. 10(b) shows that the liquid composition profile for n-C16 is almost not affected by these temperature differences while for n-C12 (Fig. 10(c)) and n-C11 (Fig. 10(d)) exhibit a similar composition trend that the hydrogen liquid composition profile. 4.5.2. Influence of the hydrocarbon feed flowrate The influence of the hydrocarbon feed flowrate (HCFF) as bifurcation parameter considering the feed composition for case 2 is shown in Fig. 11. As before, the conversions of DBT and 4,6-DMDBT are the output variables determining the steady state multiplicity existence. The MSS region is found between the HCFF values of 173.3–326 kmol/h (Fig. 11). The stable steady states with high conversion of DBT and 4,6-DMDBT are located in HCFF values smaller that 173.3 kmol/h, outside of the MSS region. It should be pointed out that, while the HCFF is increased, the H2 fed to the RDC is the same. Therefore, a lower conversion may be expected since the H2 /HC feed ratio is diminished as the HCFF is augmented. A close view of the high 4,6-DMDBT conversion exhibits an input–output multiplicity region, which shows that the 4,6-DMDBT conversion design target (99.3%) is not achieved for HCFF values greater than 120 kmol/h. Convergence problems were found when the HCFF is reduced to values where the drying of some stages is promoted (HCFF < 100 kmol/h). In such case the increment and decrement of this bifurcation parameter was performed manually.

205

Fig. 12. Bifurcation diagram for DBT conversion (case 2) with the variation of the hydrocarbon feed flowrate and the reflux ratio.

4.5.3. Combined influence of the hydrocarbon feed flowrate and the reflux ratio The combined influence of the hydrocarbon feed flowrate (HCFF) and the reflux ratio (RR) as bifurcation parameters considering the feed composition for case 2 is presented in Fig. 12 and Table 5. For this situation only the conversion of DBT is shown to be the output variable determining the steady state multiplicity existence. Nonetheless, the kinetic expressions for 4,6-DMDBT elimination are also considered in the computation of this bifurcation diagram. A multiplicity region can be identified for a combination of low reflux ratios (RR = 0.3 or RR = 0.5) and for low hydrocarbon feed flowrates (64 < HCFF < 326 kmol/h). On the other hand, for high reflux ratios (RR = 0.7 or RR = 0.8) and HCFF > 326 kmol/h, a SSS region is observed. The bifurcation curve for RR = 0.7, exhibits a sharp decrement in the DBT conversion between 478 and 490 kmol/h, indicating that near to this value (RR < 0.7) there exists an higher order singularity (Kienle & Marquardt, 2003), which is the limit to pass from only SSS regions to MSS regions. In addition, when low reflux ratios are considered, RR = 0.3, high nonlinear regions are found. Fig. 12 also shows the possible existence of an isola region consisting of medium and high DBT conversions. However, an isola which is a separate-closed homotopy curve disconnected from the main homotopy curve could not be found due to physical constraints, i.e., dried stages occur for small values of HCFF. Table 5 summarizes the type of multiplicity, the multiplicity regions, the heat duties for the condenser and reboiler (QC and QR), and a set of recommended operating HCFF values for the combined influence of the RR and HCFF. It may be observed that as the RR is augmented, the required heat duties increase due to the increment of the internal flowrates in the RDC, but also the liquid flowrate at the bottoms tends to decrease. A set of operating HCFF values are also recommended in order to achieve the conversion design targets, to keep reasonable heat duties and to reach accept-

Table 5 Combined influence analysis for case 2. Bifurcation parameters

Multiplicity type

Output multiplicity region

Heat duty

HCFF (kmol/h)

RR

DBT

4,6-DMDBT

Min. value

Max. value

QC (kW)

QR (kW)

100 100 100 100

0.3 0.5 0.7 0.8

O O X X

I–O I–O X X

59 173.3 – –

199 326 – –

3067 4698.1 6404.1 7287.9

3931.3 5562.4 7269.4 8153.3

X: No multiplicity; O: output multiplicity; I–O: input–output multiplicity.

Operability range (recommended values)

Conversion design targets are not achieved 100 ≤ HCFF ≤ 105 94 ≤ HCFF ≤ 99 94 ≤ HCFF ≤ 96

206

J. Carlos Cárdenas-Guerra et al. / Computers and Chemical Engineering 34 (2010) 196–209

Fig. 13. Bifurcation diagram for DBT and 4,6-DMDBT conversion (case 2) with the variation of the reboiler heat duty.

Fig. 14. Bifurcation diagram for DBT and 4,6-DMDBT conversion (case 2) with the variation of the holdup.

able bottoms flowrates (i.e., to maintain their values around SSS-2 given in Table 4). It should be noted that at the lowest and highest values of RR the conversion design targets are not achieved. 4.5.4. Influence of the reboiler heat duty The influence of the reboiler heat duty (QR) as bifurcation parameter, considering the feed composition for case 2 is shown in Fig. 13. In this circumstance, the conversions of DBT and 4,6-DMDBT are the output variables determining the steady state multiplicity existence. There appears in Fig. 13 that the MSS region is found between the QR values of 3298.1 and 4352 kW. The stable steady state with high conversion of DBT and 4,6-DMDBT is located when the reboiler duties are greater than 4352 kW, outside of the MSS region. It is interesting to observe that, while the DBT conversion is highly affected by the variations of the reboiler heat duty from 3500 to 4500 kW, the 4,6-DMDBT conversion is practically constant (∼99%). If a zoom of the high 4,6-DMDBT conversion is performed, an input–output multiplicity region is found, just as with RR and HCCF. Thus, 4,6-DMDBT conversion design target (99.3%) may be achieved at four different QR values, approximately at 3709.7, 4108.3, 4414.4, and 6554.7 kW.

Fig. 14 shows the influence of the HLP at the reactive zones on the performance of the RDC. An output multiplicity for DBT conversion and an input–output multiplicity for 4,6-DMDBT conversion are found at HLP values between 4300 and 6650 kg (Fig. 14). The presence of input–output multiplicity for 4,6-DMDBT conversion means that one value for 4,6-DMDBT conversion can be obtained at two or three different values of the HLP. On the other hand, regarding the three steady states for DBT conversion, two correspond to low and high conversion and are stable, while the medium conversion one is unstable. Along the branch containing unstable steady state solutions, open loop operation is not possible and the control of unstable states becomes more difficult than controlling stable states. This means that the stable steady states of the upper branch (i.e., with high conversion) outside the MSS region are a better selection as operating points (set points). Thus, the holdup should be kept at values greater than 6650 kg, but considering that the design targets must be met then HLP values greater than 8000 kg should be selected. It is important to consider that the liquid holdup or catalyst loaded at the reactive zones is rather a design parameter, so a change in such parameter would mean an operation shut down of the RDC.

4.5.5. Influence of the holdup Although the catalyst loaded at the reactive zones is a design constraint, the assumption of a completely wetted catalyst leads to consider the holdup (HLP) or catalyst load as a bifurcation parameter. The amount of catalyst loaded (up to 8000 kg) into the reactive zones is in agreement with the extremely low reaction rate expressions for 4,6-DMDBT elimination (Vanrysselberghe et al., 1998). Also, the holdup values are in the range of the catalyst loads used in the HDS conventional catalytic reactors (Furimsky, 1998).

4.5.6. Summary of the steady state multiplicity Table 6 shows the summary of MSS found for the two feed composition cases. It is clear that for case 1, when the 4,6-DMDBT is not considered, there is not multiplicity with the variation of the bifurcation parameters. However, for case 2, output and input–output multiplicities were found, and the minimum and maximum values for the bifurcation parameters defining the multiplicity region are given in Table 6. Finally, a RDC with the operating recommended values of the bifurcation parameters for case 2 are given in order to

Table 6 Multiplicity analysis summary. Bifurcation parameter

Multiplicity type

Case 1

Reflux ratio Hydrocarbon feed flowrate (kmol/h) Reboiler heat duty (kW) Holdup (kg of catalyst)

Output multiplicity region for case 2

Case 2

DBT

4,6-DMDBT

DBT

4,6-DMDBT

X X

– –

O O

I–O I–O

X X

– –

O O

I–O I–O

X: no multiplicity; O: output multiplicity; I–O: input–output multiplicity.

Min. value

Max. value

0.22 173.3

0.345 326

3298.1 4300

4352 6650

Operability range (recommended values)

0.4 ≤ RR ≤ 0.69 100 ≤ HCFF ≤ 105 4736 ≤ QR ≤ 7192 HLP ≥ 8000

J. Carlos Cárdenas-Guerra et al. / Computers and Chemical Engineering 34 (2010) 196–209

achieve the conversion design targets. It should be pointed out that the parameters with less possibility to be perturbed are the RR and the HCFF. In fact, recommended operation values for RR are around the reference on (RR = 0.5). The values for HCFF are also around the design target of 100 kmol/h, confirming that the optimal value of the H2 /HC feed ratio must be kept at 3. The reboiler heat duty used by the reference single steady state (SSS-2) still can be improved and an energy saving of around 15% may be obtained. While the amount of catalyst loaded at the reactive zones must be greater than 8000 kg in order to achieve an ULSD production.

207

was found that, the DBT conversion was highly affected by the variations of the reboiler heat duty, while the 4,6-DMDBT conversion was practically constant (∼99%). However, the 4,6-DMDBT conversion design target (99.3%) was achieved at four different reboiler duties, indicating the existence of input-output multiplicities. It was also determined that the amount of catalyst loaded at the reactive zones must be greater than 8000 kg in order to achieve an ULSD production. Also, a set of values of the bifurcation parameters are recommended to achieve an appropriate RDC operation. It is clear, that the introduction of 4,6-DMDBT in the HC feed leads to deal with highly nonlinear reactions kinetic expressions and, with this, the possibility of finding MSS is increased. Furthermore, it could be established that the absence of H2 in the reactive zones implies operating conditions similar to those of conventional (non-reactive) distillation columns with their respective MSS. On the other hand, when an excess of H2 is present in the reaction zones, the behavior of the RDC is very likely to a HDS conventional catalytic reactor. Therefore, it may be concluded that the multiplicities found in this work are highly related to the specific phenomena (reaction-separation) involved. As an overall conclusion, our findings suggest that, under optimal design and operating conditions, reactive distillation could be considered as a viable technological alternative to produce ULSD.

5. Conclusions A systematic study of the operating conditions and parameter sensibility under which MSS may occur in a RDC for diesel deep HDS has been performed. Through the analysis of the bifurcation diagrams it may be concluded that the operation of the RDC is highly sensitive to perturbations in the reflux ratio, the hydrocarbon feed flowrate, and the reboiler heat duty. The presence of 4,6-DMDBT in the feed lead to bifurcation diagrams with multiplicity existence, while only single steady state were found in its absence. For case 1, only a SSS was found. It was shown that small changes in the reflux ratio and reboiler heat duty lead to the complete elimination of DBT. In this sense, the RDC behavior is similar to that of a conventional (non-reactive) distillation column. For case 2, a multiplicity region was found between the RR values of 0.22 and 0.345. The stable steady states with design target conversions of DBT and 4,6-DMDBT were located at RR values greater than 0.345, outside of the MSS region, and an input–output multiplicity for 4,6-DMDBT at the MSS region was found. Thus, the reflux ratio is a key parameter for an appropriate RDC design and operation. Additionally, an increment in the hydrocarbon feed flowrate leads to lower conversions of the organo-sulfur compounds with the reduction of the H2 /HC feed ratio. In addition, it

Acknowledgements We thank CONACyT for the financial support through the project U45160-Y. J. Carlos Cárdenas-Guerra thanks CONACyT for the scholarship provided. Appendix A. HDS reaction rate expressions The kinetic models for hydrogenolysis ( active sites) and hydrogenation ( active sites) reaction pathways of Th, BT, DBT, and

Table A1 Reaction rate expressions and kinetic parameters for thiophene and benzothiophene hydrodesulfurization (cases 1 and 2). Component

Thiophene (Van Parijs & Froment, 1986) rTh, =

Reaction rate expression

kTh, KTh, K 2

H2 ,

1/2

(1+KH , P 2 H

2

2

+KTh, PTh +KH S, (PH S /PH )) 2 2 2



Reaction rate coefficients 6 J/kmol

kTh, = 52, 222 exp − 125.2×10 RT

 44.6×106

Adsorption equilibrium constants

Kinetic parameters

Benzothiophene (Van Parijs et al., 1986)

PTh PH

KTh, = 5.53 × 10−9 exp



 J/kmol

RT

−6

3

rBT, =

[kmol/kgcat s]

kBT, KBT, KH , PBT PH 2 2 [1+(KH , PH )1/2 +KH S, (PH S /PH )+KBT, PBT ] 2 2 2 2 2



6 J/kmol

kBT, = 0.94 exp − 73.6×10RT

[kmol/kgcat s]

KBT, = 1.90476 × 10−4

[Pa−1 ]



[kmol/kgcat s]

3

[kmol/kgcat s]

[Pa−1 ]

−6

−1

KH2 , = 3.53318 × 10 [Pa−1 ] KH2 S, = 2.08240 × 10−3 [Pa−1 ]

KH2 , = 5.22729 × 10 [Pa ] KH2 S, = 9.12 × 101 [dimensionless] Where Pi is the partial pressure of component i in Pa, T is in K, and R is 8314 J/kmol K.

Table A2 Reaction rate expressions and kinetic parameters for dibenzothiophene hydrodesulfurization (cases 1 and 2). Reaction pathway Hydrogenolysis

Hydrogenolysis and hydrogenation

Component

Dibenzothiophene (Froment et al., 1994)

Reaction rate expression

rDBT, =

Dibenzothiophene (Vanrysselberghe et al., 1998)

k1 KDBT, KH , CDBT CH 2 2



(1+KDBT, CDBT +

3

[kmol/kgcat s]

rDBT, = (1+KDBT, CDBT +



Reaction rate coefficients 6 J/kmol

k1 = 7.84 × 108 exp − 158×10RT Kinetic parameters



[kmol/kgcat s]

kDBT, KH, KDBT, CDBT CH



KH , CH +KH S, (CH S /CH )) 2 2 2 2 2

2 3

KH, CH +KBPH, CBPH +KH S, CH S +K4,6−DMDBT, C4,6−DMDBT ) 2 2 2



Adsorption equilibrium constants KDBT, = 5.31 KH2 , = 4.02 KH2 S, = 1.72

KH, = 3.36312 × 10−11 exp

[m3 /kmol]

KDBT, = 7.56868 × 101

[m3 /kmol]

 113.232×106

−8

KH2 S, = 1.47118 × 10 3

Where Ci is the molar concentration of component i in kmol/m , T is in K, and R is 8314 J/kmol K. a Data taken from Vanrysselberghe and Froment (1996).

J/kmol





RT

 48.214×106 J/kmol  RT  105.670×10  6 J/kmol

[kmol/kgcat s]a [m3 /kmol]a

[m3 /kmol]a

KBPH, = 3.84984 × 10−4 exp

[dimensionless]

6 J/kmol

kDBT, = 6.78711 × 106 exp − 122.770×10 RT

exp

RT

[m3 /kmol]a [m3 /kmol]a

[kmol/kgcat s]

208

J. Carlos Cárdenas-Guerra et al. / Computers and Chemical Engineering 34 (2010) 196–209

(m /kmol) exp

KBPH, = 4.96685 × 10

KDBT, = 2.50395 × 10−7 exp

K4,6-DMDBT, = 1.58733 × 10−8 exp

RT

 76.840×106 RT  6 KH, = 1.40255 × 10−15 exp 142.693×10 RT  37.899×106 −4

 90.485×106

[kmol/kgcat s]





2

k4,6-DMDBT, = 1.0228 × 1024 exp − 299.042×10 RT

6 J/kmol



3 2

KH, CH +KBPH, CBPH +K4,6-DMDBT, C4,6-DMDBT )

k4,6-DMDBT, KH, K4,6-DMDBT, C4,6-DMDBT CH

(1+KDBT, CDBT +

r4,6-DMDBT, =

References

Where Ci is the molar concentration of component i in kmol/m , T is in K, and R is 8314 J/kmol K. a Data taken from Vanrysselberghe and Froment (1996).

3

[m3 /kmol] K4,6-DMDBT, = 1.80397 × 101

Adsorption equilibrium constants

 Kinetic parameters

Reaction rate coefficients



(1+KDBT, CDBT +

r4,6-DMDBT, =

k4,6-DMDBT, = 1.79044 × 104 exp − 106.223×10 RT

2

6 J/kmol



[kmol/kgcat s]

[kmol/kgcat s] 2

KH, CH +KBPH, CBPH +K4,6−DMDBT, C4,6−DMDBT )

4,6-dimethyldibenzothiophene (Vanrysselberghe et al., 1998) r4,6-DMDBT = r4,6-DMDBT, + r4,6-DMDBT, [kmol/kgcat s] Component Reaction rate expression

k4,6-DMDBT, KH, K4,6-DMDBT, C4,6-DMDBT CH

Hidrogenolysis and hydrogenation Reaction pathway

Table A3 Reaction rate expressions and kinetic parameters for 4,6-dimethyldibenzothiophene hydrodesulfurization (case 2).

 J/kmol [m3 /kmol] RT  J/kmol 3 [m /kmol]a  J/kmol [m3 /kmol]a  J/kmol 3 a

[kmol/kgcat s]

4,6-DMDBT reported by Van Parijs and Froment (1986), Van Parijs et al. (1986), Froment et al. (1994), Vanrysselberghe and Froment (1996), and Vanrysselberghe et al. (1998) were used. See Tables A1–A3.

Al-Arfaj, M. A., & Luyben, W. L. (2002). Comparative control study of ideal and methyl acetate reactive distillation. Chemical Engineering Science, 57, 5039. Baur, R., Taylor, R., & Krishna, R. (2003). Bifurcation analysis for TAME synthesis in a reactive distillation column: Comparison of pseudo-homogeneous and heterogeneous reaction kinetics models. Chemical Engineering and Processing, 42, 211. Broderick, D. H., & Gates, B. C. (1981). Hydrogenolysis and hydrogenation of dibenzothiophene catalyzed by sulfided CoO-MoO3 /␥-Al2 O3 : The reaction kinetics. AIChE Journal, 27, 663. Chen, F., Huss, R. S., Doherty, M. F., & Malone, M. F. (2002). Multiple steady states in reactive distillation: Kinetic effects. Computer and Chemical Engineering, 26, 81. Ciric, A. R., & Miao, P. (1994). Steady state multiplicities in an ethylene glycol reactive distillation column. Industrial and Engineering Chemistry Research, 33, 2738. Dalal, N. M., & Malik, R. K. (2003). Solution multiplicity in multicomponent distillation. A computational study. In A. Kraslawski, & I. Turunen (Eds.), European symposium on computer aided process engineering, Vol. 13 (p. 617). Finland: Elsevier. DOE/EIA Report. (2001). The transition to ultra-low-sulfur diesel fuel: Effects on prices and supply. In Office of integrated analysis and forecasting. Washington: U.S. Department of energy., p. 13. Froment, G. F., Depauw, G. A., & Vanrysselberghe, V. (1994). Kinetic modeling and reactor simulation in hydrodesulfurization of oil fractions. Industrial and Engineering Chemistry Research, 33, 2975. Furimsky, E. (1998). Selection of catalysts and reactors for hydroprocessing. Applied Catalysis A: General, 171, 177. Gani, R., & Jørgensen, S. B. (1994). Multiplicity in numerical solution of non-linear models: Separation processes. Computers and Chemical Engineering, 18(Suppl. 18), S55. Granados-Aguilar, A. S., Viveros-García, T., & Pérez-Cisneros, E. S. (2008). Thermodynamic analysis of a reactive distillation process for deep hydrodesulfurization of diesel: Effect of the solvent and operating conditions. Chemical Engineering Journal, 143, 210. Guckenheimer, J., & Holmes, P. (1983). Nonlinear oscillations, dynamical systems and bifurcations of vector fields. New York: Springer–Verlag. Hauan, S., Hertzberg, T., & Lien, K. M. (1995). Why methyl tert-butyl ether production by reactive distillation may yield multiple solutions. Industrial and Engineering Chemistry Research, 34, 987. Hauan, S., Hertzberg, T., & Lien, K. M. (1997). Multiplicity in reactive distillation of MTBE. Computers and Chemical Engineering, 21, 1117. Houalla, M., Nag, N. K., Sapre, A. V., Broderick, D. H., & Gates, B. C. (1978). Hydrodesulfurization of dibenzothiophene catalyzed by sulfided CoO-MoO3 /␥-Al2 O3 : The reaction kinetics. AIChE Journal, 24, 1015. Hung, S. B., Lee, M. J., Tang, Y. T., Chen, Y. W., Lai, I. K., Hung, W. J., et al. (2006). Control of different reactive distillation configurations. AIChE Journal, 52, 1423. Jacobs, R., & Krishna, R. (1993). Multiple solutions in reactive distillation for methyl tert-butyl ether synthesis. Industrial and Engineering Chemistry Research, 32, 1706. Jacobsen, E. W., & Skogestad, S. (1991). Multiple steady states in ideal two-product distillation. AIChE Journal, 37, 499. Jacobsen, E. W., & Skogestad, S. (1994). Instability of distillation columns. AIChE Journal, 40, 1466. Joback, K. G., & Reid, R. C. (1987). Estimation of pure-component properties from group-contributions. Chemical Engineering Communications, 57, 233. Kannan, A., Joshi, M. R., Reddy, G. R., & Shah, D. M. (2005). Multiple-steady-states identification in homogeneous azeotropic distillation using a process simulator. Industrial and Engineering Chemistry Research, 44, 4386. Katariya, A. M., Kamath, R. S., Moudgalya, K. M., & Mahajani, S. M. (2008). Nonequilibrium stage modeling and non-linear dynamic effects in the synthesis of TAME by reactive distillation. Computers and Chemical Engineering, 32, 2243. Kienle, A., & Marquardt, W. (2003). Nonlinear dynamics and control of reactive distillation processes. In K. Sundmacher, & A. Kienle (Eds.), Reactive distillation status and future directions (p. 241). Weinheim: Wiley-VCH. Knudsen, K. G., Cooper, B. H., & Topsøe, H. (1999). Catalyst and process technologies for ultra low sulfur diesel. Applied Catalysis A: General, 189, 205. Krasnyk, M., Ginkel, M., Mangold, M., & Kienle, A. (2007). Numerical analysis of higher order singularities in chemical process models. Computers and Chemical Engineering, 31, 1100. Krishna, R. (2002). Reactive separations: More ways to skin a cat. Chemical Engineering Science, 57, 1491. Kumar, M. V. P., & Kaistha, N. (2008). Role of multiplicity in reactive distillation control system design. Journal of Process Control, 18, 692. Mangold, M., Kienle, A., Gilles, E. D., & Mohl, K. D. (2000). Nonlinear computation in DIVA—Methods and applications. Chemical Engineering Science, 55, 441. Mohl, K. D., Kienle, A., Gilles, E. D., Rapmund, P., Sundmacher, K., & Hoffmann, U. (1997). Nonlinear dynamics of reactive distillation processes for the production of fuel ethers. Computers and Chemical Engineering, 21(Suppl.), S989.

J. Carlos Cárdenas-Guerra et al. / Computers and Chemical Engineering 34 (2010) 196–209 Mohl, K. D., Kienle, A., Gilles, E. D., Rapmund, P., Sundmacher, K., & Hoffmann, U. (1999). Steady-state multiplicities in reactive distillation columns for the production of fuel ethers MTBE and TAME: Theoretical analysis and experimental verification. Chemical Engineering Science, 54, 1029. Mohl, K. D., Kienle, A., Sundmacher, K., & Gilles, E. D. (2001). A theoretical study of kinetic instabilities in catalytic distillation processes: Influence of transport limitations inside the catalyst. Chemical Engineering Science, 56, 5239. Monroy-Loperena, R., Pérez-Cisneros, E. S., & Álvarez-Ramírez, J. (2000). A robust PI control configuration for a high-purity ethylene glycol reactive distillation column. Chemical Engineering Science, 55, 4925. Müller, D., & Marquardt, W. (1997). Experimental verification of multiple steady states in heterogeneous azeotropic distillation. Industrial and Engineering Chemistry Research, 36, 5410. Peng, D. Y., & Robinson, D. B. (1976). A new two-constant equation of state. Industrial and Engineering Chemistry Fundamentals, 15, 59. Pérez-Cisneros, E. S., Granados-Aguilar, A. S., Huitzil-Meléndez, P., & Viveros-García, T. (2002). Design of a reactive distillation process for ultra-low sulfur diesel production. In J. Grievink, Schijndel, & J. van (Eds.), European symposium on computer aided process engineering, Vol. 12 (p. 301). The Netherlands: Elsevier. Sakanishi, K., Ando, M., Abe, S., & Mochida, I. (1991). Extensive desulfurization of diesel fuel through catalytic two-stage hydrotreatment. Journal of the Japan Petroleum Institute, 34, 553. Sakanishi, K., Ando, M., & Mochida, I. (1992). Extensive desulfurization of diesel fuel through catalytic two-stage hydrotreatment (Part 2). Influence of reaction pressure on desulfurization and reactivity of refractory sulfur compounds. Journal of the Japan Petroleum Institute, 35, 403. Singh, B. P., Singh, R., Kumar, M. V. P., & Naistha, N. (2005). Steady state analysis of reactive distillation using homotopy continuation. Chemical Engineering Research and Design, 83A, 959. Takatsuka, T., Inoue, S. I., & Wada, Y. (1997). Deep hydrodesulfurization process for diesel oil. Catalysis Today, 39, 69. Taylor, R., & Krishna, R. (2000). Modelling reactive distillation. Chemical Engineering Science, 55, 5183. Vadapalli, A., & Seader, J. D. (2001). A generalized framework for computing bifurcation diagrams using process simulation programs. Computers and Chemical Engineering, 25, 445.

209

van Hasselt, B. W., Lebens, P. J. M., Calis, H. P. A., Kapteijn, F., Sie, S. T., Moulijn, J. A., et al. (1999). A numerical comparison of alternative threephase reactors with a conventional trickle-bed reactor. The advantages of countercurrent flow for hydrodesulfurization. Chemical Engineering Science, 54, 4791. Van Parijs, I. A., & Froment, G. F. (1986). Kinetics of hydrodesulfurization on a CoMo/␥-Al2 O3 catalyst. 1. Kinetics of the hydrogenolysis of thiophene. Industrial and Engineering Chemistry Product Research and Development, 25, 431. Van Parijs, I. A., Hosten, L. H., & Froment, G. F. (1986). Kinetics of hydrodesulfurization on a CoMo/␥-Al2 O3 catalyst. 2. Kinetics of the hydrogenolysis of benzothiophene. Industrial and Engineering Chemistry Product Research and Development, 25, 437. Vanrysselberghe, V., & Froment, G. F. (1996). Hydrodesulfurization of dibenzothiophene on a CoMo/Al2 O3 catalyst: Reaction network and kinetics. Industrial and Engineering Chemistry Research, 35, 3311. Vanrysselberghe, V., Gall, R. L., & Froment, G. F. (1998). Hydrodesulfurization of 4methyldibenzothiophene and 4, 6-dimethyldibenzothiophene on a CoMo/Al2 O3 catalyst: Reaction network and kinetics. Industrial and Engineering Chemistry Research, 37, 1235. Viveros-García, T., Ochoa-Tapia, J. A., Lobo-Oehmichen, R., de los Reyes-Heredia, J. A., & Pérez-Cisneros, E. S. (2005). Conceptual design of a reactive distillation process for ultra-low sulfur diesel production. Chemical Engineering Journal, 106, 119. Wang, J., Chang, Y., Wang, E. Q., & Li, C. Y. (2008). Bifurcation analysis for MTBE synthesis in a suspension catalytic distillation column. Computers and Chemical Engineering, 32, 1316. Wang, S. J., Wong, D. S. H., & Lee, E. K. (2003). Effect of interaction multiplicity on control system design for a MTBE reactive distillation column. Journal of Process Control, 13, 503. Yang, B., Wu, J., Zhao, G., Wang, H., & Lu, S. (2006). Multiplicity analysis in reactive distillation column using ASPEN PLUS. Chinese Journal of Chemical Engineering, 14, 301.