Development and Validation of a Reduced-Index Dynamic Model of

Subscriber access provided by READING UNIV

Article

Development and Validation of a Reduced-Index Dynamic Model of an Industrial High-Purity Column Jose O.A. Matias, Misagh Ebrahimpour, and Galo A. Carrillo Le Roux Ind. Eng. Chem. Res., Just Accepted Manuscript • DOI: 10.1021/acs.iecr.7b02789 • Publication Date (Web): 15 Dec 2017 Downloaded from http://pubs.acs.org on December 16, 2017

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a free service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are accessible to all readers and citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

Industrial & Engineering Chemistry Research is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 45 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Industrial & Engineering Chemistry Research

Development and Validation of a Reduced-Index Dynamic Model of an Industrial High-Purity Column Jose O. A. Matias, Misagh Ebrahimpour, and Galo A. C. Le Roux∗ Department of Chemical Engineering, University of Sao Paulo, Sao Paulo E-mail: [email protected]

Abstract A depropanizer column is a high purity distillation unit, whose modeling is challenging due to its high nonlinear behavior and intricate dynamics. Also, depending on the assumptions, a differential algebraic equation (DAE) system of high index can arise in the modeling of the column. Such systems create significant difficulties to numerical simulation. Phenomenological hydraulic relationships can be used to avoid the index problems. However, these relationships may rely on unknown and arbitrarily assumed parameters and project details. This paper presents an alternative modeling approach called Dynamic Design Using Proportional Action (DDPA), which is used to simplify the hydraulic and pressure drop relation for a distillation unit. DDPA approach generates a reduced index dynamic model by adding equations similar to a proportional controller with arbitrarily large gain. A depropanizer column was modeled using DDPA and the model was validated with real process data obtained in an industrial-scale depropenizer column located at Refinaria de Paul´ınia, Paul´ınia, S˜ ao ∗

To whom correspondence should be addressed

1

ACS Paragon Plus Environment

Industrial & Engineering Chemistry Research 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Paulo, Brazil and owned by Petrobras S.A.. The results show that the model with DDPA is able to predict the actual column dynamic behavior properly.

Introduction Dynamic simulation and modeling are tools that are useful throughout the entire plant lifetime: from plant design (e.g. operability and controlability of a plant; testing operating procedures; environmental studies; etc.) to plant operation (e.g. computer based operator training; validation of safety procedures, etc.) 1,2 Dynamic modeling of distillation columns can be carried out by using different approaches with distinct levels of complexity to take into account the hydrodynamics and vapor-liquid equilibrium (VLE). Many contributions involving model simplification strategies are available in the literature, i.e: linearization; 3 compartmental models; 4,5 polynomial approximation; 6–8 neural network for adjusting experimental data; 9 to name a few. However, a more detailed description of column’s thermo and hydrodynamics may be necessary. In this case, a rigorous model should be used. Rigorous modeling strategies use a tray-by-tray model with mass- and energy balances, and equilibrium relations on each tray. Also, additional equations are applied to describe pressure and flow dynamics. 10 Nevertheless, even rigorous modeling may also include simplifications, such as perfect mixing in both phases, thermal and thermodynamic equilibrium between phases, vapor flows dynamics negligible, constant molar flows, etc.. Clearly, each of these simplifications has consequences for the model. Table 1 summarizes how different assumptions about pressure and holdup influence the differential algebraic equation (DAE) model. A DAE model is basically a combination of differential and algebraic equations. The differential equations are used to describe the dynamic balances of mass, components and energy. In turn, the algebraic constraint set usually describes pressure and temperature equilibrium between liquid and vapor phases and also thermodynamic

2


Page 2 of 45

Page 3 of 45 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


equilibrium relations. 11 Fundamentally, DAE systems are distinct from ordinary differential equation (ODE) systems. 12 The index of a DAE system provides a conceptual measurement of this difference. The system index can be defined as the minimum number of times the algebraic equations need to be differentiate in order to obtain an equivalent ODE system. 13 Table 1: VLE Modeling strategies - Distillation Columns Molar Holdup

Pressure

Equation or Assumption

Consequence

- Francis weir formula for flow through weirs - gas pressure drop through plate

Variable

Variable

Fixed

Fixed

- liquid and vapor flow constant in sections - no heat of mixing

Fixed

Fixed

- differential equations added to the model: dM/dt = 0 and dP/dt = 0

DAE model is Index 1 Validity of the correlation and parameter accuracy are not guaranteed. DAE model is Index 1 McCabe-Thiele hypothesis. No energy balances in trays. Oversimplification of heat relations. DAE model is Index 2 McCabe-Thiele hypothesis are not necessary. Liquid and vapor flow change through the column.

Examples Gani et al. 14 , Lockett 15 , and most commercial simulators McCabe and Thiele 16 Jogwar and Daoutidis 17

Albet et al. 18

As seen in Table 1, different assumptions for the description of the hydraulic behavior lead to DAE models of different indexes. If the DAE system is index one, it can be easily rearranged and treated as an ODE system. Then, it can be solved with ODE solvers or conventional DAE solvers. 13 If the system is index two or more, the model index can be either reduced to one by algorithms that reduce the index automatically 19,20 , or efficiently solved by stiffly stable solvers, like IDAS solver included in SUNDIALS. 21 Nevertheless, for modeling distillation columns, both options need higher order derivatives of the thermodynamic properties, which are not normally provided by conventional thermodynamic packages and are quite difficult to obtain manually. Additionally, high index DAE systems bring numerical challenges for the simulation, e.g. finding a set of consistent initial conditions 22 ; and numerical inaccuracies in enforcing the algebraic constraints. 11 Hence, in practice, it is interesting to have an index one model. According to Table 1, for the three most commonly used models in the literature, index one models can be obtained either by using inaccurate algebraic relations for pressure drop and tray hydraulics,

3



like Francis weir formula, or by assuming ideal solution behavior and constant holdup. In the first option, the use of phenomenological relations for tray hydraulic and pressure drop equations give rise to errors as a result of unknown project details and parameters. For the second option, the model assumptions make liquid and vapor flows to be constant in each column section such that the energy balance in the stages is not needed. This assumption applies only if there is no heat of solution, and it is unrealistic for many challenging problems in literature. A different approach, which is proposed in this paper, is to use a simple linear relation between pressure/vapor flowrate and molar holdup/liquid flowrate offering a simple equation to relate tray hydraulics and pressure drop. This linear relation is similar to a proportional controller with arbitrarily large gain. The proposed approach was called Dynamic Design Using Proportional Action (DDPA). The main advantage of the DDPA is that the tray model with the DDPA equations is still index 1, which allows the use of thermodynamic package in order to calculate complex equilibrium relations accurately. To highlight the commitments involved on choosing the modeling approach, the steps are summarized in Figure 1. The paper is organized as follows. In the next section mathematical considerations to the DDPA equations is presented. Next, two case studies are described, a simple flash tank and, the process of interest, a vapor recompression distillation unit. In the Case Studies sections, both are analyzed in order to tune the parameters of the linear relations included in the model and the approach is validated using actual process data. Next, some conclusions of the performance of the approach are drawn in the last section.

Mathematical Considerations To exemplify the DDPA strategy and the sources of the index problem, the following system can be considered: 12

4


Page 4 of 45

Page 5 of 45 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Figure 1: Modeling options based on available solvers and information

M(x)x˙ = f (x)

(1)

In which, M(x) is a state-dependent matrix called mass matrix. If this matrix is nonsingular, which indicates an ODE system, the model can be easily rearranged as:

x˙ = M−1 f (x)

(2)

However, if the mass matrix is singular (i.e. the system is composed by a combination of differential and algebraic equations), the rearrangement strategy is different. Considering the following system:

y˙ = F(y, z) (3) 0 = G(y, z) Then, the equations may be expressed in the form 1. 5



Page 6 of 45



     I 0 y˙   F(y, z)  Mx˙ =    =   = f (x) 0 0 z˙ G(y, z)

(4)

If (∂G/∂z T ) is nonsingular, it is easy to determine a set of consistent initial conditions, i.e. one that satisfies all algebraic equations. Simply by setting y, and computing z using Newton’s method to solve: 12

0 = G(y, z)

(5)

Also, this nonsingularity allows the system to be rearranged in order to obtain a ODE system. 



d  y˙ = F(y, z)  ¨=  ⇒y dt 0 = G(y, z)

∂F ∂yT

y˙ +

∂F ∂zT

˙ 0= z,

∂G ∂yT

    I 0 F(y, z)    y˙  =        ∂G ∂G 0 − F(y, z) ˙ z ∂zT ∂yT

y˙ +

∂G ∂zT

z˙

(6)



(7)

The difficulty in solving the DAE system depends on the singularity of ∂G/∂zT , which relates to the index problem defined previously. 23 If ∂G/∂zT is singular, the system has an index 2 or higher. ODEs are considered to have index 0. For a more concrete example, consider the modeling of distillation column trays, given in Table 1. Assuming pressure and liquid holdup constant is the same as adding ML = MLref and P = P ref to the algebraic equation set. The other algebraic equations are equilibrium equation and thermodynamic relations. However, with this assumption, vapor (V ) and liquid (L) outlet flowrates do not appear in the algebraic set. Because both are algebraic variables (i.e. they do not explicit differential equations), ∂G/∂zT becomes singular, generating the index problem. It can be removed by not assuming ML and P constant, and introducing L = f (ML ) and V = f (P ) to the 6


Page 7 of 45 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


algebraic set. For further discussion please refer to Ebrahimpour 24 . As discussed previously, the correlations for calculating the pressure drop (vapor flow rate) and liquid flow rate in the trays outlet have this form and will also decrease the system index to 1. But, as explained previously, they usually rely on unknown and arbitrarily assumed parameters and project details. Therefore, the motivation to develop DDPA appears. DDPA intends to approximate the tray hydraulics by a ”perfect control” strategy in the form of:

L − Lref = KM (ML − MLref )

V − V ref = KP (P − P ref )

(8)

Note that in the cases that the gain, KP or KM , is large, the changes in L in relation to Lref do not affect the difference (ML − MLref ) to a great extent. However, the inverse statement is not true. If ML differs from MLref , the value of (L − Lref ) is drastically affected. Hence, P and ML tend to their reference value, satisfying the DDPA equations, while the liquid and vapor flows can change according to the thermodynamic calculations. Another advantage of the DDPA approach is that the relationship between the properties involved appears explicitly in the model, highlighting the phenomenological relationship between them.

Case study: Flash Tank To illustrate the index problem, its numerical solution difficulties and the advantages of the methodology proposed, a flash tank example was chosen, Figure 2. For the flash model, the following points were considered: Phases are well mixed and adiabatic operation. Three different flash models were proposed, all the equations are described in the Appendix. The models contain: mass, components, and energy balances; mechanical, thermal and thermodynamic equilibrium of the streams leaving the flash; and the component fraction 7



Figure 2: Flash tank summation. The difference between the models lies in the way the pressure (P ) and liquid holdup (ML ) are calculated. In the first model, both variables are explicitly specified (ML = MLref and P = P ref ), resulting in an index-2 model. In the second formulation, instead of directly specifying the holdup and the pressure, valve pressure drop equations as well as a level controller are added to the model, which generates an index-1 model. Finally, to illustrate the DDPA formulation, Equations 8 are used to replace the explicit specification of pressure and holdup. Taking advantage of the equation structure, ML and P become indirectly specified. However, in contrast to the first model, this approach results in a index-1 model, avoiding the inherent higher index problems.

Flash Tank Results The methodology was applied to a simple equilibrium system, a flash tank, to provide some insights about tuning the K values. To assess the impact of different gain on the solution, the flash model with DDPA equations was compared to: a high-index model (simulated with the DASSLC solver 25 embedded in EMSO); and the model with the phenomenological

8


Page 8 of 45

Page 9 of 45 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


relationships plus controller (for brevity, named as correlation model). All the models are described in Appendix 1. The flash tank, Figure 2, was simulated in two different scenarios. In the first one, the simulation starts at a steady state operation point and a step-like disturbance increased the feed flowrate from its nominal value (F ref = 20kmol/h) to 21kmol/h (variation of 5%). The second scenario is similar, a step-like disturbance occurred in the light component fraction of the feed, resulting in a increase of 20 % in the steady state nominal value (xref = 0.5). f Both disturbances occurred at t = 2.77h (103 s). The assessment of both scenarios is based on comparing the simulated profiles of: liquid holdup/liquid outlet; pressure/vapor outlet; and compositions. In the case of the pressure and liquid molar holdup (i.e.: values that should be fixed in the reference value by the DDPA equations), the analysis was extended to the comparison of average error. The error is calculated by dividing the absolute value of the deviation from the index-2 model by the current value of the variable. Then, the values are averaged over time. These results are shown in Tables 2 and 3. Additionally, since the variables were normalized previously, the difference between KP and KM was not tested. Both gain values were specified at the same value and the subscript was removed (KP = KM = K). For the first scenario, it is possible to notice that the behavior of the liquid and vapor outlet flows is similar for all models except for a small overshooting after the disturbance (Figure 3). In the liquid outlet, the overshooting happens due to level controller action of the correlation model. Similarly in the vapor outlet, this behavior happens due to the fact that the DDPA equation with small gain behaves as a proportional controller. The most remarkable result comes from the analysis of the pressure and liquid holdup. For all simulations, deviations from the index-2 model solution can be noticed. However, analyzing Table 2, the order of magnitude of this imprecision is small compared to the actual value of the property, especially in the DDPA cases with gains greater than ≥ 5 (deviation 1%) and most importantly, the deviation decreases as K grows, which supports the previous

9



statement that the model with DDPA equation with large gain behaves as an index-2 models. In relation to the pressure deviation, the deviations of the correlation model and the DDPA model with K = 1 have the same order of magnitude. That happens due to the form of the equations. Instead of approaching the index-2 model, both equations emulate the behavior of proportional controllers. Hence, a small offset from the P ref (P ext in the correlation model) is expected in the steady-state value after the disturbance. Flow disturbance

Flow disturbance

17.6

6.78 6.76

17.2

Liquid holdup [kmol]

Liquid outlet [kmol/h]

17.4

Index 2 K=1 K=5 K = 10 K = 20 K = 40 Correlations

17

16.8

16.6

16.4 0

6.74 Index 2 K=1 K=5 K = 10 K = 20 K = 40 Correlations

6.72 6.7 6.68 6.66

2

4

6 Time [h]

8

10

6.64 0

12

2

4

(a) Liquid outlet

6 Time [h]

8

10

12

(b) Liquid holdup

Flow disturbance

3.7

Flow disturbance

302

3.65 301.5 Index 2 K=1 K=5 K = 10 K = 20 K = 40 Correlations

3.6

3.55

Pressure [psia]

Vapor outlet [kmol/h]

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 10 of 45


301

300.5 3.5

3.45

300 0

2

4

6 Time [h]

8

10

12

0

2

(c) Vapor outlet

4

6 Time [h]

8

10

12

(d) Pressure

Figure 3: Dynamic response of the flash model to disturbance case 1 In the second scenario (Figure 4), the behavior of the liquid holdup, liquid/vapor outlet flow and pressure is very similar to the first scenario and the same remarks are valid. However, due to the nature of the disturbance, composition profiles changed in time in this scenario 10


Page 11 of 45 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Table 2: Percentage of deviation from high-index solution, scenario 1 Liquid holdup [%] Pressure [%] K=1

1.11

0.26

K=5

0.23

5.26 e−2

K = 10

0.11

2.63 e−2

K = 20

5.66 e−2

1.32 e−2

K = 40

2.83 e−2

6.58 e−3

Correlations

1.23 e−2

0.27

and also the temperature profiles, which are not shown because they are similar to the composition profiles. In relation to the composition, the method performance is satisfactory even in cases with small values of K. The deviations from the index-2 model simulation for pressure and liquid holdup are presented in Table 3. They were calculated as in the previous table. The results are also similar to the first scenario, but a small increase in the deviation of the pressure to its reference is observed. However, no substantial differences were found in any of the variables analyzed. The comparison of the different models in both scenarios provides additional support for the DDPA approach. They show that the model with the DDPA equations can reproduce the behavior of the high-index model with an average error of less than 1% when the gain values are greater 5. On the other hand, large values of K are prone to inflate rounding errors in the model. For example, a small imprecision of the calculation of pressure can change the value of vapor flow in one or two orders of magnitude, implying in solver failure. Additionally, the fact that the augmented system with DDPA equations is equivalent to an index-2 model should be taken in account while choosing the solver. The solution of the augmented model with an inappropriate solver can be inaccurate, unstable and cause oscillations in the solution. Hence, careful attention must be paid on the trade-off caused by the K value. For

11



the processes of interest (system that present a time-scale separation between flows and composition, e.g. VRD), the authors experience suggests that K = 20 is a satisfactory initial value for the gain and that K values larger than 40 do not show a noteworthy increase in the precision. Table 3: Percentage of deviation from high-index solution, scenario 2 Liquid holdup [%] Pressure [%] K=1

2.67

2.77

K=5

0.52

0.57

K = 10

0.26

0.29

K = 20

0.13

0.14

K = 40

6.51 e−2

7.20 e−2

Correlations

2.86 e−2

3.8

Vapor Recompression Distillation (VRD) Unit Modeling Vapor recompression distillation (VRD) is an energy integrated distillation configuration 17 and its applications are widely discussed since the decade of 1980 (Muhrer et al. 26 and the references therein) due to the significant savings in energy consumption that the VRD offers for close-boiling liquid separations. The energy savings come mostly from the reboiler-condenser combination. The vapor that comes from the top of the column is used to vaporize the bottom stream. A compressor is needed to increase the pressure of the top vapor stream. The increase in the pressure allows the vapor phase-change, facilitating the heat transfer in the combined reboiler-condenser. 17 This configuration leads to an integrated dynamic between energy addition (reboiler-duty) or removal (condenser duty), i.e. they cannot be varied independently. 26 Another characteristic of VRD systems is the large number of trays. This particular 12


Page 12 of 45

Page 13 of 45

Composition disturbance


18

6.7

17.5 6.65

16.5

Liquid holdup [kmol]

Liquid outlet [kmol/h]

17


16 15.5 15 14.5 14

6.6


6.55

6.5

6.45

13.5 13 0

2

4

6 Time [h]

8

10

6.4 0

12

2

4

(a) Liquid outlet

8

10

12


320

5

315 Index 2 K=1 K=5 K = 10 K = 20 K = 40 Correlations

4.5

4

Pressure [psia]

Vapor outlet [kmol/h]

6 Time [h]

(b) Liquid holdup


5.5

3.5


310

305

300

3

295

0

2

4

6 Time [h]

8

10

12

0

2

4

6 Time [h]

(c) Vapor outlet

8

10

12

(d) Pressure



0.64

0.0165

0.63

0.016

0.62 0.61

Fraction C2H6 [−]

0.0155 Fraction nC4 [−]

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60



0.6 0.59 0.58


0.015 0.0145 0.014

0.57 0.0135

0.56 0.55 0

2

4

6 Time [h]

8

10

0.013 0

12

(e) Heaviest component on liquid outlet, molar fraction

2

4

6 Time [h]

8

10

12

(f) Lightest component on vapor outlet, molar fraction

Figure 4: Dynamic response of the flash model to disturbance case 2

13



Figure 5: Vapor recompression distillation column with controllers

14


Page 14 of 45

Page 15 of 45 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


design leads to time constants of composition loops that are much slower than pressure loops. 26 Consequently, composition dynamics are nearly unaffected by the flow dynamics. 27 A modeling peculiarity of VRD modeling is that the vapor hold up must be considered due to high operating pressure. In columns operating at moderate pressures (less than 10 atm), ignoring vapor hold up is a good assumption, because vapor densities are much lower than liquid densities. This means that most of the material on the trays is in the liquid phase, despite the fact that the volume of vapor is typically a factor of 10 greater than the volume of liquid. However, as pressure increases toward the critical pressures of the components being distilled, the difference between liquid and vapor densities decreases, and vapor hold up becomes more important. 28 In propylene-propane system, vapor hold up represents about 30% of the total molar hold up and has to be considered.

Modeling The methodology to obtain a reduced index model was applied to an actual process and validated with actual process data. In this work, rigorous tray by tray model based on vapor and liquid thermodynamic equilibrium was developed. Distillation trays were modeled dynamically. Sump and accumulator dynamics follow the same model, however top and bottom products control the levels of accumulator and sump respectively. Dynamics of expansion valves, trim heat exchanger, compressor and pump were considered fast and were neglected. 29 For the reboiler-condenser, the dynamic in the cold and hot side was also neglected, but the dynamic of the heat transfered between both sides was considered and modeled as a first order process. 30 The system was modeled in EMSO 31 and the equations are listed in Appendix . Distillation columns dynamic models have extra degrees of freedom in comparison to stationary models due to the inclusion of the holdup calculations. As seen in Results and Discussion Section, phenomenological relations can be included or one may specify the pressure and liquid holdup in the column trays to meet this extra degrees of freedom. The last 15



Page 16 of 45

alternative is simpler to model, however it leads to higher index models. The simulation software used solves higher index problems. But, the thermodynamic package is unable to supply higher order derivatives of the thermodynamic properties. According to Figure 1, an alternative is to develop a reduced order dynamic model of the column (i.e. higher order derivatives will not be necessary). VRD has some characteristics that justify the specification of pressure and liquid holdup, i.e. tightly controlled pressure and composition dynamics that are much slower than the flow dynamics. 27 The equations added to the model are:

ref Lnorm − Lref norm = KM (LevelL,norm − LevelL,norm )

ref ref Vnorm − Vnorm = KP (∆Pnorm − ∆Pnorm )

(9)

(10)

In which, V and L are the vapor and liquid molar flows leaving the tray; LevelL is the liquid level of the tray; ∆P is the difference of pressure between the tray and the one above; The subscript norm indicates that the variable was normalized (the normalization method is explained in Appendix , Equation 29); The superscript ref is related to the reference values, i.e. the specification of pressure difference and level; KP and KM are the gains of DDPA, as explained above. Choosing large values of K in Equations 9 and 10 is equivalent to keep the value of pressure and liquid holdup constant. These large values generate a stiff system that tends to an index-2 system when K tends to infinity. In addition to the gain adjustment, data reconciliation and parameter estimation were performed. These steps allow the model to represent correctly the operational points. The experimental data used in these steps was obtained from an industrial-scale depropenizer column owned by Petrobras SA. The dynamic simulation was performed for more than two days of operation in order to verify the quality of model response compared to real plant

16


Page 17 of 45 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


data. Due to the fact that data reconciliation and parameter estimation use a stationary model, steady state detection was applied in order to identify the stationary experimental point. The SS detection approach created by Jiang et al. 32 was applied. This routine is based on wavelet transform and determines, for multivariable processes, one global steady state index with a value between 0 and 1. In the period of times in which the index value is equal to one, the process is assumed to be at steady state. In this case study the proposed dynamic model uses the adjusted vapor efficiency values from steady state parameter estimation and a modified value of the flow rate of stream 12 (all the stream names used throughout the paper are related to Figure 5), modified based on the steady state data reconciliation of the initial point. For solving data reconciliation and parameter estimation problems, a simplex optimization routine was applied. Both steps were carried out using EMSO.

Dynamic Simulation For the dynamic simulation, simultaneous variations in reflux flow rate (stream 5 ), flow rate of reboiler hot side (stream 12 ), feed flow rate (stream 1 ) and composition of feed were considered. Flow rates of bottom product (stream 3) and distillate product were specified as the plant data, since during the data acquisition they were manipulated manually in the plant so as to control the sump and accumulator levels.

Vapor Recompression Distillation Results The dynamic behavior of the VRD model with the DDPA equations was studied in two cases. In the first case, step changes in the manipulated variables were carried out. In the second one, the plant data was used as input to the model. In both studies, the dynamic response of key variables of the model was analyzed.

17



Model validation: case 1 The dynamic system has two manipulated variables, reflux flow rate (stream 5) and flow rate of the hot stream through the reboiler side (stream 12). In the plant, the specification of propylene in the top product is controlled by manipulating the reflux flow rate. In turn, the loss of propylene in bottom is regulated by the manipulation of reboiler/reflux flow ratio. Hence, plant data does not allow to check the effect of changing each manipulated variable separately, which is crucial for the model validation. To analyze the individual effects of the manipulated variables in the bottom and top compositions, a qualitative analysis was carried out, i.e. the dynamic simulation was not compared to actual plant data, instead process knowledge was used to inspect the simulation results For the model validation, distillate (stream 4) and bottom product (stream 3) flow rates are manipulated applying PI controllers to keep the sump and accumulator levels close to their set points values. The model specification and controller parameters are shown in the support Supporting Information Section. The model is shown in Appendix . The DDPA equations were tuned based on the analysis of the flash, the values of K used in Equations 10 (pressure equation) and 9 (level equation) are listed in Table 4. The normalization of the variables (based on Equation 29) was carried out using the equipment data in order to determine the maximum and minimum value of the pressures, flow rates and levels. Table 4: K values used for simulation K value KP KM KP Sump KM KP Accumulator KM Column

20 20 40 20 40 5

The qualitative analysis was divided in two parts. First, a step of −4.2t/h was applied 18


Page 18 of 45

Page 19 of 45

to stream 5 (flow rate of the hot stream) and the dynamic behavior of the top and bottom composition was analyzed, Figure 6.

Bottom propylene composition xbottom [%mol]

45 40 35 30 25 20 0

5

10

15

20

25

30

35

40

45

35

40

45

Time [h] Top propane composition 1800 1600

xtop [ppm]

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


1400 1200 1000 800 0

5

10

15

20

25

30

Time [h]

Figure 6: Simulation results, step of -4.2 t/h in reboiler hot side stream Reduction of the flow rate of the hot stream through the reboiler side implies in less heat entering the column. The immediate effect is the decrease of internal vapor flow through the column and, consequently, the total molar flow rate at the bottom product increases (while top molar flow rate decreases). In turn, molar composition of propylene in both products grows. At the bottom, this increase is justified by the lower rate of vaporization of the feed, which is mainly composed by propylene (around 75%). Hence, more propylene is exiting through the bottom. At the top, despite the overall decrease in the molar flow, propane is heavier than propylene, hence the proportion of propane at the top product decreases as shown. For the second part, a step with the same value was applied in stream 12 (reflux flow 19



rate) and the composition behavior of top and bottom composition was also analyzed, Figure 7.

Bottom propylene composition xbottom [%mol]

25 20 15 10 5 0 0

5

5

#10

10

15

20

25

30

35

40

45

35

40

45

Time [h] Top propane composition

4

4

xtop [ppm]

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 20 of 45

3 2 1 0 0

5

10

15

20

25

30

Time [h]

Figure 7: Simulation results, step of -4.2 t/h in reflux stream The reduction of reflux flow rate decreases the top product purity as expected. However, it also causes an increase in the total flow rate of the top product, because more vapor is flowing up the column, as the there is less liquid to be heated up. Such increase implies that the total molar flowrate of propylene leaving the column at top grows despite the fact that composition decreases. In addition, less propylene is entering the column with the reflux. These two effects result in the reduction of the propylene composition at the bottom. Additionally, comparing both parts of this case study shows that the effect of the hot stream flow rate (stream 12) reduction is faster than the effect of the decrease of the reflux flow rate, specially on the top product composition. For an additional analysis, in the case in which reflux flow rate changes (Figure 8) high20


Page 21 of 45

lights the difference between column internal composition and flowrate dynamics. Composition curve shows a much slower dynamics compared to the flow rate. This effect is in accordance with the hypothesis necessary to apply the extra equations (instead of specifying the pressure and holdup values) explained in Methods Section.

Internal liquid flow in tray #100

L100 [kmol/h]

8080 8075 8070 8065 8060 8055 0

5

10

15

20

25

30

35

40

45

Time [h] Internal composition of propylene in liquid in tray #100 0.958 0.956

x100 [-]

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


0.954 0.952 0.95 0.948 0.946 0

5

10

15

20

25

30

35

40

45

Time [h]

Figure 8: Simulation results, step of -4.2 t/h in reflux stream. Comparison between internal flow and composition dynamic responses

Model validation: case 2 In the second case, the dynamic simulation was performed for more than two days of operation in order to verify the quality of model response compared to real plant data. The K values used were the same as in the previous studies. Steady state detection, data reconciliation and parameter estimation To compare the behavior of the DDPA model with actual plant data, a previous step was 21



Page 22 of 45

carried out in order to fit the model to the measurements obtained from a real process. This data set contains all the available measurements of the vapor recompression unit during a month, with frequency of 1 minute. To identify the stationary periods, the data set was analyzed off-line using the method proposed in. 32 The variables chosen to represent the steady state were the flow rate and temperature of: feed stream, (stream 1) 1 , and both process outputs streams: the top product, (stream 4), and the bottom product, (stream 3). After detecting the stationary periods, data reconciliation was performed based on an weighted least square objective function

Fob =

X (Ai − Ai,measured )2 σA2 i

i

(11)

In which, the residual between measured and predicted values is minimized. The variables chosen were: flow rates of streams 1, 3, 5 and 12 ; the summation of ethane and propane composition in the top of the column, stream 4 ; and the propylene composition in the bottom of the column, stream 3. The variance of the measured variables is chosen as the weight factor for two main reasons. First, the errors present heteroscedasticity due to the different nature of measurements and instruments, i.e. the variance of the errors is not constant, and this effect should be taken in account in the reconciliation. Second, measurements with small variance contain more information and, as the weight is inversely proportional to the variance, observations with small variance should have a large weight in the objective function. The data was adjusted to allow the accurate representation of the real plant data. The data reconciliation results are shown in Table 5. The only major adjustment was in the reflux hot side stream flow rate (stream 12 ). To avoid subsequent reconciliations, only the correction of the flow value obtained at the initial point was kept constant, by implying a constant bias during the simulation. The constant bias can be seen in Figure 9d. 1

all the stream names in this section are related to Figure 5

22


Page 23 of 45 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Table 5: Reconciliated and planta data Plant data Reconciliated data Stream Stream Stream Stream

1 Flow 3 Flow 5 Flow 12 Flow

rate rate rate rate

[t/h] [t/h] [t/h] [t/h]

36.212 11.267 364.94 350.85

34.401 10.711 364.26 360.39

In order to reproduce bottom and top product specifications and compressor discharge temperature, parameter estimation was performed based on experimental data. Vapor Murphree efficiency of stripping and rectification zones, and compressor efficiency were adjusted. The vapor efficiencies were chosen as adjustable parameters due to non ideal behavior of propylene-propane mixture in high pressures. 26 Mass transfer characteristics in distillation processes are commonly taken into account using vapor Murphree efficiencies. These efficiencies depend on flow regime, tray layout, and physical properties of the mixture. The adjusted parameters values are shown in the Supporting Information. The manipulated variables were approximated by a series of step functions as seen in Figure 9. Additionally, the flow rate of stream 12 (Figure 9d) was modified based on the reconciled data of the dynamic simulations initial point. Temperature variations of boil up, temperature of trays 17, 85 and 171, composition of propylene in distillate (stream 4) and propane in the bottom product (stream 3) were compared with real plant data in Figures 10 and 11. Analyzing the behavior of propylene bottom composition (Figure 11a) and entrances (feed flow rate, Figure 9a, and feed propylene composition, Figure 9b), a significant correlation is observed. The feed flowrate step around 200 minutes causes the first considerable increase in the mole fraction of propylene in the bottom. Then, the mole fraction of propylene in the feed drops around 300 minutes, the top purity increases, and consequently the bottom composition decreases. After, the joint effect of the decrease of feed flowrate and propylene composition in the feed causes a steep decrease in the bottom propylene composition between 800 23



Feed flow rate

Mole fraction of propylene in feed stream 0.78

42

0.77

40

0.76 Fraction [-]

Flow [t/h]

44

Plant data Approximation

38


0.75

36

0.74

34

0.73 0.72

32 0

1000

2000 Time [min]

0

3000

(a) Feed flow rate, stream 1, t/h

1000

2000 Time [min]

3000

(b) Mole fraction of propylene in feed stream 1

Reflux flow rate

370

Reboiler hot side flow rate

500

368

450

366

400

364

Flow [t/h]

Flow [t/h]

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 24 of 45


362


350 300

360

250

358 356

200 0

1000

2000 Time [min]

3000

0

(c) Reflux flow rate, stream 5, t/h

1000

2000 Time [min]

3000

(d) Reboiler hot side flow rate, stream 12, t/h

Figure 9: Manipulated variables approximation in comparison to actual plant data

24


Page 25 of 45

Temperature of tray #17

295

297

294

296

293 Plant data Simulation data

292


298

Temperature [K]

Temperature [K]

296

291


294 293

290

292

289

291

288

290 0

1000

2000 Time [min]

3000

0

(a) Temperature of tray 17, K

1000

2000 Time [min]

3000

(b) Temperature of tray 85, K


Temperature of boil up 310

302

300

Temperature [K]

Temperature [K]

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Plant data Simulation data

298

296


300

295 294 0

1000

2000 Time [min]

3000

0

(c) Temperature of tray 171, K

1000

2000 Time [min]

3000

(d) Temperature of boil up, K

Figure 10: Dynamic response of the VRD model in comparison to actual plant data

25



Mole fraction of propylene in bottom product 0.35

Concentration of propane in top product 6000

0.3 Concentration [ppm]

5000

0.25 Fraction [-]

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 26 of 45

Plant data Simulation data

0.2 0.15 0.1


3000 2000 1000

0.05

0 0

1000

2000 Time [min]

3000

0

(a) Mole fraction of propylene in bottom product stream 3

1000

2000 Time [min]

3000

(b) Propane in top product stream 4, ppm

Figure 11: Dynamic response of the model bottom and top compositions in comparison to actual plant data and 1200 minutes. The behavior in bottom propylene until, approximately, 2500 minutes can also be explained by the entrances variation. However, after changing the reflux flow rate and flow rate of reboiler hot side, their effects on the dynamic responses dominated the feed effects. Hence, after 2500 minutes, the dynamic response of the column bottom temperature (Figure 10d), and propylene concentration in two products (bottom - Figure 11a, top Figure 11b) can be better explained according to the changes in reflux flow rate (Figure 9c) and thereafter flow rate of reboiler hot side (Figure 9d). First, the decrease of both flow rates led to a significant increase in concentration of propane in the top product. However, the decrease of reflux rate decreased the propylene amount (mole flow) inside the column and consequently at the bottom product. Next, as the flow rate through reboiler hot side continued to go down, its effect on the products composition showed up. It led to less heat entering the column, increasing the propylene concentration at the bottom product and decreasing the propane concentration at the top product, which characterized the peak seen in Figure 11a. This behavior of the composition is in accordance to the results obtained in the previous validation case.

26


Page 27 of 45 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Due to the fact that temperature is intrinsically connected with composition, the small deviation seen in the temperatures profile (Figure 10) is very important indeed because it is an indicator of changes in the concentrations of top impurities in a high-purity distillation column. The top tray temperatures were approximately constant during the studied period. On the other hand, the trays closer to the bottom are more influenced by entrances and manipulated variables changes. The peak seen between 2500 and 3000 in Figures 10b, 10c, and 10d is related to the abruptly change in the reflux flow rate and flow rate of reboiler hot side. A second peak can be seen in Figure 10d, which is connected to the effect of the decrease of feed flowrate and propylene composition between 800 and 1200 minutes. Due to the fact that more propane enters the system, this heavier component goes down the column, causing a slight decrease of the propylene mass fraction in the bottom and, consequently, decreasing the equilibrium temperature in the reboiler (the pressure is kept constant). As shown in Figures 10 and 11, the model can predict the dynamic behavior of the plant properly. The main differences between model predicted values and plant data are: in propane composition at the top product, however, before 2500 minutes, plant data was clearly affected by instrument failure; and in the bottom composition, in which the model trajectory have a qualitative difference from the plant data. The simulated profile increases while the plant data decreases. One possible reason is the linear approximation of the dynamics of the reboiler-condenser. In other periods, there are several sources that could cause this loss in the model accuracy, manly regarding unmeasured disturbances and thermodynamics, for example: uncertainties in composition measurements and in flow rate measurements; neglecting some minor components in the feed; and inaccuracies in thermodynamic parameters, like binary coefficient for the pair propylene/propane. However, apart from the lack of accuracy, the process dominant time constants of both

27



plant and model match within an adequate precision. Especially for control purposes, the response curves and the steady-state changes should be represented rather than the absolute value actual process variable. We believe that this lack of accuracy does not affect the main conclusions of the paper, which is the validation of the VRD model using the DDPA approach.

Conclusion In a simple system (the flash tank), the simulated trajectories of the reduced-index model obtained with the DDPA strategy were compared to the ones simulated by higher index solvers, and those obtained from the model with phenomenological correlations. The trajectories were similar in the cases in which the gains of the DDPA equations were large, above 10. A VRD column was modeled based on the DDPA strategy to produce a reduced-index model. Its capability for representing the column dynamics is demonstrated in two situations. First, by analyzing the dynamic behavior of the column by applying step changes in manipulated variables. Second, by comparing the simulated data with the real plant data during a long operation period. In both situations the DDPA approach showed promising results. Future works could involve the application of the dynamic model as a virtual plant for tuning a Real Time Optimization (RTO) software and also for analyzing its performance. Additionally, Moving Horizon Estimation (MHE) could be implemented using the dynamic model for estimating the system states.

28


Page 28 of 45

Page 29 of 45 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Appendix: Flash Model implemented in the equationoriented environment This example was extracted from, 33 Exercise 7.21. The model was implemented in EMSO and the components considered were: C2 H4 , C2 H6 , C3 H6 , C3 H8 , nC4 , iC4 . The feed condition was fixed (F = 1 kmol/h, hf = 13.210 BTU/mol and z = [0.02, 0.03, 0.05, 0.1, 0.6, 0.2]). The flash volume was arbitrarily chosen as 1.6 m3 The nominal values of liquid holdup and pressure were 6.6 kmol and 300.2 psia, respectively. The model equations are: d(ML + MV ) =F −V −L dt

d(ML xi + MV yi ) = F zi − V yi − Lxi dt

(i = 1, ..., ncomp )

d(ML hL + MV hV − MV RT ) = F hf − V hV − LhL dt

ˆ 1/2 = c1,i + c2,i T + c3,i T 2 , h L,i

with T [o R] and i = 1, . . . , ncomp

(12)

(13)

(14)

(15)

ncomp

X

ˆ L,i = hL xi h

(16)

i=1

ˆ 1/2 = e1,i + e2,i T + e3,i T 2 , h V,i

with T [o R] and i = 1, . . . , ncomp

(17)

ncomp

X

ˆ V,i = hv yi h

(18)

i=1

ncomp

X

xi = 1

i=1

29


(19)


Page 30 of 45

yi = ki xi , i = 1, . . . , ncomp

(20)

1/3 ki = a1,i + a2,i T + a3,i T 2 + a4,i T 3 , with T [o R] and i = 1, . . . , ncomp T

(21)

V ltotal = V lV + V lL

(22)

where, F , L and V are respectively the feed, liquid outlet, and vapor outlet molar flows; ML and MV are the liquid and vapor holdups; z, x, and y are the molar composition of the feed, ˆ i is the specific liquid outlet, and vapor outlet; T is the temperature; P is the pressure; h vapor enthalpy of the component i; h is the specific enthalpy of the stream, subscripts f , L and V indicate feed, liquid outlet, and vapor outlet; ki is the equilibrium constant of component i; Vl , Vv , and Vt are the liquid, vapor and total volume of the flash; c1 , c2 , c3 , e1 , e2 , e3 , a1 , a2 , a3 , and a4 are constants and their values are shown in Tables 6, 7 and 8. Table 6: Equilibrium constant coefficients at 300 psia Component

a1 × 102

a2 × 105

C2 H4 C2 H6 C3 H6 C3 H8 nC4 iC4

-5.1779950 -9.8400210 -25.098770 -14.512474 -14.181715 -18.967651

62.124576 67.545943 102.39287 53.638924 36.866353 61.239667

a3 × 108

a4 × 1012

-37.562082 8.0145501 -31.459290 -9.0732459 -75.221710 153.84709 -5.3051604 -173.58329 16.521412 -248.23843 -17.891649 -90.855512

For the first approach, liquid holdup and the flash pressure are fixed in their reference values:

P = P ref

(23)

ML = MLref

(24)

30


Page 31 of 45 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Table 7: Liquid enthalpy coefficients at 300 psia Component

c1

c2 × 10

c3 × 105

C2 H4 C2 H6 C3 H6 C3 H8 nC4 iC4

-7.2915 -8.4857 -12.4279 -14.50006 -20.29811 -16.553405

1.5411962 1.6286636 1.8834652 1.9802223 2.3005743 2.1618650

-1.6088376 -1.9498601 -2.4839140 -2.9048837 -3.8663417 -3.1476209

Table 8: Vapor enthalpy coefficients at 300 psia Component C2 H4 C2 H6 C3 H6 C3 H8 nC4 iC4

e2 × 10

e1

e3 × 105

56.79638 615.93154 2.4088730 61.334520 588.7543 11.948654 71.828480 658.5513 11.299585 81.795910 389.81919 36.47090 152.66798 -1153.4842 140.64125 147.65414 -1185.2942 152.87778

In the second approach, the pressure drop equations for the liquid and vapor outlet valve were used to replace Equations 23 and 24. Also, a PI controller was implemented in order to keep the level constant. The valve pressure drop equations that were applied are the standard models provided by EMSO: r ∗ V = Cv,V

(P − Pext ) G

∗ L = Cv,L xovalve

p

LevelL

(25)

(26)

∗ In which, P ext is fixed at 290 psia; G is the specific gravity of the gas (G = 0.03); Cv,L ∗ and Cv,V are the valve sizing coefficients with values 0.5 and 30, respectively; LevelL is the

flash level; and xvalve is the valve opening, which was determined by a level controller. Next, the proposed approach, DDPA, used simple linear relations to replace Equations 23 and 24. As a consequence, a pseudo index-2 model was obtained, in which liquid holdup

31



Page 32 of 45

and the flash pressure are still fixed in their nominal values, but the model is actually an index 1 model.

ref ref Vnorm − Vnorm = KP (Pnorm − Pnorm )

(27)

ref Lnorm − Lref norm = KM (ML,norm − ML,norm )

(28)

The pressure, flowrates and pressure applied in the DDPA equations were normalized as shown below for the pressure variable, enabling to use similar values of K, even when relationships between different properties are studied.

Pnorm =

P − Pmin Pmax − Pmin

(29)

Appendix: VRD Model implemented in the equationoriented environment This section summarizes the equations used to model the vapor recompression distillation column. The dynamic model is presented in this appendix, the stationary model used for reconciliation and parameter estimation is simple derivation of it. The model was also implemented in EMSO (the codes in EMSO language can be seen in Ebrahimpour 24 ) and the components considered were: propylene and propane. The distillation column (Equations 30 - 40), expansion valves (Equations 44 - 49), compressor (Equations 51 - 54), and reboiler/condenser (Equations 55 - 67) equations are listed. The column sump and the accumulator were also modeled, but due to the similarities with the tray model, their models were omitted. All the equipments and their arrangement are shown in Figure 5. The adiabatic distillation column was modeled tray-by-tray, numbered from top, 1, to 32


Page 33 of 45 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


bottom, ntray = 197. The tray model comprises 3ncomp + 6 equations (where ncomp is the number of components, which is equal to 2): mass, component and energy balances (Equations 30, 31, and 32); thermodynamic, mechanical and temperature equilibria of the streams leaving the system (Equations 33, 34, and 35); tray efficiency (Equation 37 and 36); summation (Equation 38); tray volume constraint (Equation 40); and also the equations relative to the DDPA. For the DDPA the reference values of pressure and pressure difference are shown in the Supporting Information. The reference flow values were not shown. In order to obtain these values, the stationary solution needs to be calculated based on the initial state shown in the Supporting Information. The holdup of the trays was assumed to be 1 kg mol−1 . The notation used: F , V and L are respectively the molar flow rate of the feed, and internal vapor and liquid streams. For the trays that do not have a feed stream the value of F is equal to 0. All the variables related to the feed properties have the subscript f . x and y are the mole fractions of liquid and vapor phases, k the equilibrium constant. hL and hV are the liquid and vapor enthalpies and ν is the specific volume of the liquid. The latter values were calculated with the thermodynamic package embedded in EMSO. Additionally, the values of Murphree Efficiencies of the stripping section, EsM V , and rectification section ErM V were estimated based on plant data. The tray area, Atray , and volume, V ltray , were obtained directly from plant documentation. The values of the sump area/volume and accumulator area/volume were also obtained from plant documentation. The parameter values are shown in the Supporting Information.

d(ML,j + MV,j ) = FV,j + FL,j + Vj+1 + Lj−1 − Vj − Lj dt

(j = 1, ..., ntray )

(30)

d(ML,j xj,i + MV,j yj,i ) = FV,j yf,j,i + FL,j xf,j,i + Vj+1 yj+1,i + Lj−1 xj−1,i − Vj yj,i − Lj xj,i dt (31) (j = 1, ..., ntray ), (i = 1, ..., ncomp ) 33



Page 34 of 45

d(ML,j hL,j + MV,j hV,j − Pj V ltray,j ) = FV,j hV,f,j + FL,j hL,f,j + Vj+1 hV,j+1 + Lj−1 hL,j−1 − Vj hV,j − Lj hL,j dt (j = 1, ..., ntray ) (32)

? kj,i xj,i = yj,i

ErM V =

EsM V =

(j = 1, ..., ntray ), (i = 1, ..., ncomp )

(33)

TV,j = TL,j

(j = 1, ..., ntray )

(34)

PV,j = PL,j

(j = 1, ..., ntray )

(35)

yj,i − yj+1,i ? yj,i − yj+1,i

yj,i − yj+1,i ? yj,i − yj+1,i

(j = 1, ..., nf eed − 1), (i = 1, ..., ncomp )

(36)

(j = nf eed − 2, ..., ntray ), (i = 1, ..., ncomp )

(37)

ncomp

X

? yj,i =1

(j = 1, ..., ntray )

(38)

yj,i = 1

(j = 1, ..., ntray )

(39)

i=1

ncomp

X i=1

V ltray,j = ML,j νL,j + MV,j νV,j

Leveltray,j =

ML,j νL,j Atray,j

34

(j = 1, ..., ntray )

(j = 1, ..., ntray )


(40)

(41)

Page 35 of 45 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


ref Lnorm − Lref norm = KM (LevelL,norm − LevelL,norm )

(42)

ref ref ) = KP (∆Pnorm − ∆Pnorm Vnorm − Vnorm

(43)

The expansion valves were modeled as a adiabatic vapor-liquid equilibrium process. The equation used are: mass, component and energy balances, Equations 44, 45, and 46; phase equilibrium, Equations 47, 48, and 49; and summation, Equation 50. The equations used to model the valves comprise 2ncomp + 6. Due to the faster dynamic of the valves, stationary models were used. The auxiliary cooler is also modeled as the expansion valves.

F −V −L=0

F zi − V yi − Lxi = 0

(j = 1, ..., ncomp )

F hF − V hV − LhL = 0

ki xi = yi

(i = 1, ..., ncomp )

(44)

(45)

(46)

(47)

TV = TL

(48)

PV = PL

(49)

ncomp

X

yi = 1

i=1

35


(50)


Page 36 of 45

The compressor was modeled based on isentropic efficiency. First, the enthalpy of the isentropic process, hisen is calculated (Equation 54). Then, using the efficiency, the enthalpy of process is calculated finding the actual outlet temperature (Equation 54). The other equations involved in the compressor modeling are: material balances and pressure change through compressor. Due to the faster dynamic of the compressor in comparison to the column dynamics, stationary models were applied.

F in − F out = 0

yiin − yiout = 0

(i = 1, ..., ncomp )

(51)

(52)

Pdrop = P out − P in

(53)

hisen − hin hout − hin

(54)

ηCP =

The reboiler-condenser modeling has some special issues. Both sides were modeled as stationary process, however the dynamic of the equipment was modeled by a first order lag of the heat transfer between the cold and the hot side (Equation 67). As shown in Figure 5, the cold side is the liquid coming from the bottom of the column, indicated by subscript C. The hot side is the overhead vapor coming from the top of the column, denoted by the subscript H. The cold side outlet streams is modeled as saturated vapor while the hot side can leave the reboiler as saturated or subcooled liquid depending on the specification of the temperature. Both cold and hot side equations sets comprise material balances, equilibria, summation and energy balances.

FHin − FHout = 0

36


(55)

Page 37 of 45 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


in − xout yH,i H,i = 0

(i = 1, ..., ncomp )

(56)

Pdrop,H = PHout − PHin

(57)

vf B FCin − VCout = 0

(58)

(1 − vf B )FCin − Lout C = 0

(59)

out out out out FCin xin C,i − VC yC,i − LC xC,i = 0

out kC,i xout C,i = yC,i

(j = 1, ..., ncomp )

(i = 1, ..., ncomp )

(60)

(61)

out out TV,C − TL,C =0

(62)

out out PV,C − PL,C =0

(63)

ncomp

X

xout C,i = 1

(i = 1, ..., ncomp )

(64)

i=1

out out QB,H + FHin hin H,V − FH hH,L = 0

(65)

out out QB,C + FCin hf,C − VCout hout C,V − LC hC,L = 0

(66)

37



τRB

dQB,C + QB,C = −QH,C dt

Page 38 of 45

(67)

Finally, the first case of the VRD model validation uses controllers to maintain the level of accumulator and sump on their setpoints. Both are PI controllers. The sump level controller was named as bottom controller and the accumulator level controller as top controller. The sump and the accumulator were both modeled as trays.

Laccum =

Lsump =

SP Biastop + Kp,top (Levelaccum − Levelaccum )+

SP Biasbottom +Kp,bot (Levelsump −Levelsump )+

1

Z

Ki,top

1 Ki,bottom

SP (Levelaccum − Levelaccum ) (68)

Z

SP (Levelsump −Levelsump ) (69)

Nomenclature The symbols were kept constant throughout the paper. However, different units were used depending on the case study. In the variable specification, the first unit was used in the flash example and the second in the VRD example. A= transversal tray area, m2 a = equilibrium constant coefficients, units change depending on the coefficient Bias = controller bias c = liquid enthalpy coefficients, units change depending on the coefficient ∗ Cv,V = vapor valve sizing coefficient, kmol h− 1 psia−1/2 ∗ Cv,L = vapor valve sizing coefficient, kmol h− 1 m−1

e = vapor enthalpy coefficients, units change depending on the coefficient ErM V = Murphree Efficiency of rectification section

38


Page 39 of 45 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


EsM V = Murphree Efficiency of stripping section F = feed molar flow, kmol h−1 Fob = parameter estimation objective function G = specific gravity of the gas leaving the flash h = specific enthalpy, BT U mol−1 and kJ kmol−1 ˆ i = specific enthalpy of the component i, BT U mol−1 h k = equilibrium constant K = DDPA equation gain - subscript indicates the variable Kp = controller proportional gain Ki = coefficient of integral term L = liquid molar flow, kmol h−1 LevelL = liquid level, m M = molar holdup, kmol ncomp = number of components nf eed = feed tray ntray = number of trays P = pressure, psia and atm Pdrop = pressure drop, atm Pext = external pressure, psia ∆P = pressure difference between , atm QB = reboiler heat load, kJ R = gas constant, BT U kgmol−1 o R−1 s = specific entropy, kJ kmol−1 K −1 t = time, h T = temperature, o R and K vf B = vapor fraction in the reboiler hot side V = vapor molar flow, kmol h−1

39



V l = volume, m3 x = liquid molar composition xmass = liquid mass composition xovalve = valve opening y = vapor molar composition z = feed molar composition (liquid + vapor) ηCP = compressor efficiency τRB = reboiler heat exchange time constant, h ν = specific volume molar, m3 kmol−1

Subscripts and superscripts accum

= variables relative to the accumulator model

bottom

= variables relative to the bottom controller

C

= reboiler cold side - bottom of the column

f

= feed related variables = reboiler hot side - top of the column

H i

= component index

j

= tray index

L

= variables related to liquid state

max

= variable upper bound

min

= variable lower bound

norm

= normalized value of a variable - according to Equation 29

sump

= variables relative to the sump model

V

= variables related to vapor state

in

= related to variables entering the control volume

isen out

- isentropic operation = related to variables entering the control volume

40


Page 40 of 45

Page 41 of 45 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


ref

= reference value

SP

= setpoints

?

= equilibrium composition

Acknowledgement The authors thank CNPq and Petrobras S.A. for the financial support.

References (1) Seider, W. D.; Seader, J. D.; Lewin, D. R. Product and process design principles: synthesis, analisys and evaluation; John Wiley and Sons, 2009. (2) Luyben, W. L. Process modeling, simulation and control for chemical engineers; McGraw-Hill Higher Education, 1989. (3) Gilles, E. D.; Retzbach, B. Reduced models and control of distillation columns with sharp temperature profiles. Decision and Control including the Symposium on Adaptive Processes, 1980 19th IEEE Conference on. 1980; pp 865–870. (4) Benallou, A.; Seborg, D. E.; Mellichamp, D. A. Dynamic compartmental models for separation processes. AIChE J. 1986, 32, 1067–1078. (5) Musch, H. E.; Steiner, M. Order reduction of rigorous dynamic models for distillation columns. Comp. Chem. Eng. 1993, 17, S311–S316. (6) Stewart, W. E.; Levien, K. L.; Morari, M. Simulation of fractionation by orthogonal collocation. Chemical Engineering Science 1985, 40, 409–421. (7) Pinto, J. C.; Biscaia, E. Order reduction strategies for models of staged separation systems. Comp. Chem. Eng. 1988, 12, 821–831.

41



(8) Huss, R. S.; Westerberg, A. W. Collocation methods for distillation design. 1. Model description and testing. Ind. Eng. Chem. Res. 1996, 35, 1603–1610. (9) Osorio, D.; Pérez-Correa, R.; Belancic, A.; Agosin, E. Rigorous dynamic modeling and simulation of wine distillations. Food Control 2004, 15, 515–521. (10) Skogestad, S. Dynamics and control of distillation columns-a critical survey. Modeling, identification and control 1997, 18, 177–217. (11) Daoutidis, P. Surveys in Differential-Algebraic Equations II ; Springer, 2015; pp 69–102. (12) Beers, K. J. Numerical methods for chemical engineering: applications in Matlab; Cambridge University Press, 2006. (13) Brenan, K. E.; Campbell, S. L.; Petzold, L. R. Numerical solution of initial-value problems in differential-algebraic equations; Siam, 1996. (14) Gani, R.; Ruiz, C. A.; Cameron, I. T. A generalized model for distillation columns I: Model description and applications. Comp. Chem. Eng. 1986, 10, 181–198. (15) Lockett, M. J. Distillation tray fundamentals. 1986, (16) McCabe, W. L.; Thiele, E. W. Graphical design of fractionating columns. Ind. Eng. Chem. Res. 1925, 17, 605–611. (17) Jogwar, S. S.; Daoutidis, P. Vapor recompression distillation: Multi-scale dynamics and control. American Control Conference, 2009. ACC’09. 2009. (18) Albet, J.; Le Lann, J. M.; Joulia, X.; Koehret, B. Operational policies for the start-up of batch reactive distillation column. INSTITUTION OF CHEMICAL ENGINEERS SYMPOSIUM SERIES. 1994; pp 63–63.

42


Page 42 of 45

Page 43 of 45 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(19) Pantelides, C. C.; Gritsis, D.; Morison, K. R.; Sargent, R. W. H. The mathematical modelling of transient systems using differential-algebraic equations. Comp. Chem. Eng. 1988, 12, 449–454. (20) Louren¸co, M.; Secchi, A. R. Solu¸caõ de sistemas de equa¸co˜es algébrico-diferenciais ordinárias de ´ındice superior. Revista Liberato 2008, 9 . (21) Hindmarsh, A. C.; Brown, P. N.; Grant, K. E.; Lee, S. L.; Serban, R.; Shumaker, D. E.; Woodward, C. S. SUNDIALS: Suite of nonlinear and differential/algebraic equation solvers. ACM Transactions on Mathematical Software (TOMS) 2005, 31, 363–396. (22) Pantelides, C. C. The consistent initialization of differential-algebraic systems. SIAM Journal on Scientific and Statistical Computing 1988, 9, 213–231. (23) Ascher, U. M.; Petzold, L. R. Computer methods for ordinary differential equations and differential-algebraic equations; Siam, 1998. (24) Ebrahimpour, M. Dynamic and Steady-State Modeling of VRD Column in EquationOriented Environment. M.Sc. thesis, Universidade de São Paulo, 2015. (25) Costa, E.; Secchi, A.; Biscaia, E. Automatic integration of high-index dynamic systems. Computer Aided Chemical Engineering 2002, 10, 865–870. (26) Muhrer, C. A.; Collura, M. A.; Luyben, W. L. Control of vapor recompression distillation columns. iecres 1990, 29, 59–71. (27) Levy, R. E.; Foss, A. S.; Grens, E. A. Response modes of a binary distillation column. Industrial & Engineering Chemistry Fundamentals 1969, 8, 765–776. (28) Choe, Y. S.; Luyben, W. L. Rigorous dynamic models of distillation columns. Ind. Eng. Chem. Res. 1987, 26, 2158–2161. (29) Luyben, W. L. Use of dynamic simulation to converge complex process flowsheets. Comments from the Editors 2004, 43



(30) Gross, F.; Baumann, E.; Geser, A.; Rippin, D. W. T.; Lang, L. Modelling, simulation and controllability analysis of an industrial heat-integrated distillation process. Comp. Chem. Eng. 1998, 22, 223–237. (31) Soares, R. P.; Secchi, A. R. EMSO: A new environment for modelling, simulation and optimisation. Computer Aided Chemical Engineering 2003, 14, 947–952. (32) Jiang, T.; Chen, B.; He, X.; Stuart, P. Application of steady-state detection method based on wavelet transform. Comp. Chem. Eng. 2003, 27, 569–578. (33) Henley, E. J.; Seader, J. D. Equilibrium-stage separation operations in chemical engineering; Wiley,, 1981.

44


Page 44 of 45

Page 45 of 45 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Graphical TOC Entry

45


Development and Validation of a Reduced-Index Dynamic Model of

Recommend Documents