- Email: [email protected]
Rigorous dynamic modeling and simulation of wine distillations
Food Control 15 (2004) 515–521 www.elsevier.com/locate/foodcont
Rigorous dynamic modeling and simulation of wine distillations Daniel Osorio, Ricardo Perez-Correa *, Andrea Belancic, Eduardo Agosin Department of Chemical and Bioprocess Engineering, P. Universidad Cat olica de Chile, P.O. Box 306––22, Santiago, Chile Received 23 November 2002; received in revised form 11 August 2003; accepted 11 August 2003
Abstract A novel simulation strategy for dynamic distillation of complex mixtures, such as wine, is proposed and evaluated in terms of computing efficiency and accuracy. The model developed describes wine distillation as a multicomponent reactive batch distillation process. The simulation approach transforms the system of differential algebraic equations (DAE) into a set of ordinary differential equations, by pre-solving the algebraic equations and replacing them with artificial neural networks. This new simulation strategy for wine distillation is 40% faster than the rigorous solution of the DAE system, compared at the same level of accuracy. The model can be applied to the distillation of other spirits or complex mixtures, as well as in other separation processes in which the recovery of aromas is essential. 2003 Elsevier Ltd. All rights reserved. Keywords: Reactive-batch distillation; Terpenes; Pisco
1. Introduction Distilling wine for the production of spirits is an early application of reactive batch distillation. Nowadays, challenges regarding this highly non-linear process, as design improvements, automatic control and optimization, all demand efficient and accurate simulation tools. Although various rigorous batch distillation models are available (Alejski & Dupart, 1996; Ruiz, Basualdo, & Scenna, 1995; Schneider, Noeres, Kreul, & G orak, 2001), solving these models requires the application of special purpose differential-algebraic equation (DAE) solvers, which are computationally intensive and may present initialization problems. Published work on the modeling of wine distillation is scarce and mainly focuses on the much simpler continuous system (Berna, Bon, Sanju an, & Mulet, 1995; Lora, Iborra, Perez, & Carbonell, 1992). Bosley and Edgar (1994) introduced a functional approach for batch distillation modeling called direct dynamic derivatives (DDD), which avoids the need for sophisticated and expensive DAE solvers. It was
*
Corresponding author. Tel.: +56-582-686-4258; fax: +56-562-6865803. E-mail address: [email protected] (R. Perez-Correa). 0956-7135/$ - see front matter 2003 Elsevier Ltd. All rights reserved. doi:10.1016/j.foodcont.2003.08.003
employed successfully to model water/ethanol batch distillation. However, the simple vapor/liquid equilibrium model and polynomial fitting they used is limited and unsuitable for simulating multicomponent distillation. In this work we explore the use of neural networks to improve the accuracy, performance and range of applicability of the DDD approach, making it suitable for wine distillation. A simplified multicomponent mixture was used representing wine.
2. Modeling wine distillation process Thermodynamically, wine is a highly non-ideal colloidal suspension that contains several hundred components. Its main constituents are water and ethanol, although the many volatile aromas that are found in remarkably small proportions define the quality of the wine and its distillate. Among these trace components, which we denote here as ‘‘minor components’’, the main positive quality factors in wine distillates are esters, terpenes and aromatic alcohols. In turn, methanol, fatty acids, C6-components and superior alcohols are the most deleterious to distillate quality. During the distillation process, volatile components are extracted from wine and concentrated in the spirit, with the aim to recover the maximum amount of ethanol and positive
516
D. Osorio et al. / Food Control 15 (2004) 515–521
Nomenclature Notations a heat transfer area, m2 c concentration, mg/l Cp heat capacity, J/C kg dw weir diameter, cm EA energy of activation EM Murphy efficiency FR coolant flowrate, l/min H molar enthalpy of vapor, J/mol h molar enthalpy of liquid, J/mol hw weir height, cm k0 Arrhenius constant kw Francis weir constant, mol/s cm5=2 L liquid flowrate, mol/s l liquid level on the tray, cm M molar holdup, mol MW molecular weight, g/mol NC number of components P total pressure, mmHg P0 pure component vapor pressure, mmHg Q heat transfer rate, J/s R universal gas constant, J/mol K r net reaction rate, mol/s S reaction order T temperature, C
characteristic aromas, while minimizing the off-flavors in the distillate. As a case study, this work focused on the distillation of Pisco brandy, a wine distillate from Chile and Peru with a distinctive Muscat and fruity aroma, whose character is mainly attributable to its high terpene content (Loyola, Almy, Martın-Alvarez, & Cabezudo, 1987). In Muscat wines, terpenes are found either free or bound to sugars as glycosides (Baumes, Bayonove, & Gunata, 1994). The latter are termed aroma precursors, as they are unable to express their aromatic character, although when hydrolysed during distillation these precursors do constitute a further source of free terpenes (Ohta, 1995). Distillation temperature and low pH conditions however oxidize terpenes, driving them to less positive or negative odorant forms. The behaviour of linalool, the most important terpene found in Pisco, is analysed in detail in this work.
t U V W x y z
time, s global heat transfer coefficient, J/s m2 vapor flowrate, mol/s polynomial coefficients matrix liquid phase composition, mol/mol vapor phase composition, mol/mol state variables vector
Greek characters c activity coefficient q density, g/ml U apparent liquid volume, ml / polynomial fit terms W mass matrix Subindex n tray i component eth ethanol out outlet partial condenser coolant stream in inlet partial condenser coolant stream Superscripts * ideal T transpose
• Constant pressure through the column; • Perfect mixing occurs in each tray; • Vapor phase accumulation is neglected. In wine distillates, water and ethanol account for over 96% of the total mass. Minor components, while key to spirit quality, have an insignificant influence upon the thermodynamic properties of the mixture. Hence, distillate density, enthalpy, molecular weight and boiling temperature were obtained from water–ethanol binary expressions. Furthermore, the heat of reaction and the heat of mixing were also neglected. This case study considers a simplified mixture including eight of the most relevant components in the aromatic profile of Pisco brandy: ethanol, water, methanol, octanoic acid, ethyl hexanoate and free, bounded and degraded linalool. 2.2. Fundamental equations
2.1. Model assumptions
Total mass, component, and energy balances (Eqs. (1)–(3)) are presented for a generic tray n.
Given the small size of the distillation column generally used in this process (6–10 trays, 0.4 m diameter), the following assumptions are considered applicable:
C X dMn ¼ Vnþ1 þ Ln1 ðVn þ Ln Þ þ rn;i dt i¼1
N
ð1Þ
D. Osorio et al. / Food Control 15 (2004) 515–521
Mn
Table 1 Parameters for the Francis hydraulic equation
dxn;i ¼ Ln1 ðxn1;i xn;i Þ Vn ðyn;i xn;i Þ dt NC X þ Vnþ1 ðynþ1;i xn;i Þ þ rn;i xn;i rn;i ;
Francis weir expression Ln ¼ kw dw ðln hw Þ1;5
i¼1
8i ¼ 1; . . . ; NC 1 Mn
ð2Þ
dhn ¼ Ln1 ðhn1 hn Þ Vn ðHn hn Þ dt NC X rn;i þ Vnþ1 ðHnþ1 hn Þ þ Qn hn
ð3Þ qn ¼
Since only NC 1 balances are independent, water balance was not considered. The net reaction rate for component i, rn;i , is the difference between generation and consumption rates. The heat transfer term, Qn , is positive for the reboiler and negative for the partial condenser while it is neglected for intermediate trays (the column is insulated).
The Raoult law was applied to find the theoretical equilibrium composition of the gas phase, yn;i , using Antoine expression to obtain vapor pressures of pure components, P 0 . As the sum of the gas fractions on each tray must be equal to one, the Raoult equation and the consistency equation can be combined into the following single algebraic equation: ð4Þ
Experimental and theoretical studies have shown that activity coefficients of organic minor components in water/ethanol mixtures depend on the ethanol content of the mixture only (Conner, Brikmyre, Patterson, & Piggott, 1998). The UNIFAC method (Hansen, Rasmussen, Fredenslund, Schiller, & Gmehling, 1991) was used to estimate the activity coefficients of minor components in terms of ethanol composition, not available in literature. Activity coefficients for water and ethanol were estimated using the more accurate Van Laar expression (Prauznitz, Lichtenthaler, & Gomes de Azevedo, 1999). The actual composition in the gas phase, yn , was calculated introducing the global efficiency Murphy factor, EM . yn;i ¼ ynþ1;i þ EM ðyn;i ynþ1;i Þ
6 32 3
xn;eth MWeth þ ð1 xn;eth ÞMWw 103 Un xn;eth þ ð1 xn;eth ÞMWw =qw
ð6Þ
where Un is the apparent liquid molar volume, that for water/ethanol mixture was obtained from the following expression (Franks, 1972) Un ¼ 5:1214 102 þ 6:549 103 xn;eth þ 7:406 105 Tn
ð7Þ
Enthalpy expressions for the water/ethanol mixture were taken from Neuburg and Perez-Correa (1994)
2.3. Vapor–liquid equilibrium
NC 1 X 0 c xn;i Pn;i ¼0 Pn i¼1 n;i
kw (mol/s cm5=2 ) dw (cm) hw (cm)
(1974) to calculate volumetric accumulation required in Francis equation,
i¼1
1
517
ð5Þ
2.4. Constitutive equations Liquid flowrates were calculated using Francis’ weir expression (Perry, Green, & Maloney, 1999), using parameters shown in Table 1. Liquid density was modeled using the method recommended by Murphy and Gaines
Hn ¼ 44765:71 6171:03yn;eth þ ð31:46 11:98yn;eth Þ Tn þ ð4:063 104 þ 0:073yn;et Þ Tn
ð8Þ
hn ¼ ð55:678 xn;eth þ 75:425Þ Tn 0:0057 xn;eth 0:00125
ð9Þ
2.5. Partial condenser Reflux ratio is generated by partially condensing the vapor outlet stream at the top of the column using a heat exchanger device. The heat transfer in this partial condenser (tray number 0) is assumed at pseudo-steadystate, therefore the energy accumulation term is neglected. As a result, the energy balance applied to the coolant stream is reduced to an algebraic equation that must be solved iteratively for Tout to obtain the partial condenser heat absorption. Q0 FR Cp ðTout Tin Þ ðTout Tin Þ ¼ aU lnððT0 Tin Þ=ðT0 Tout ÞÞ
ð10Þ
2.6. Reaction kinetics For the reactive components, in this case bounded, free and degraded linalool, a simple chain reaction model was employed. It is assumed that reactions take place in the liquid phase in each tray. Reaction rates are calculated using an experimentally fitted Arrhenius kinetic expression obtained in this work (see Figs. 1 and 2). Parameter values are shown in Table 2.
518
D. Osorio et al. / Food Control 15 (2004) 515–521
Fig. 1. Free linalool reaction kinetic at different temperatures.
Fig. 2. Bounded linalool reaction kinetic at different temperatures.
the implicit non-linear algebraic equations (Eqs. (4) and (10)) are pre-solved using standard non-linear optimization routines, and then polynomials are fitted to their solution. This avoids solving the algebraic system repeatedly at every step of the integration procedure. Moreover, it eliminates the requirement of consistent initial conditions, an inconvenient of DAE formulation that becomes harder to overcome as the number of variables increases. According to the vapor-liquid equilibrium (VLE) equations, the mixture composition depends only on equilibrium temperature and pressure. Consequently, the system of VLE equations may be solved independently for temperature for a wide range of composition and pressure values, and then an empirical non-linear expression can be fitted to the solution. As a result, the equilibrium temperature and enthalpy expressions can be evaluated explicitly in terms of tray pressure and composition using the adjusted function. The replacement allows the direct evaluation of internal vapor flowrates, since temperature derivatives are eliminated from the energy balance (see Bosley & Edgar, 1994, for details). Therefore, and grouping composition as vectors, (i.e. xn ¼ fx1 ; x2 ; . . . ; xn g),
oh T oh T PNC Ln1 hn1 hn ox ðxn1 xn Þ þ Vnþ1 Hnþ1 hn ox ðynþ1 xn Þ þ Qn hn i¼1 rn;i Vn ¼ oh T Hn hn ox ðyn xn Þ
3. Simulation The model described above (Eqs. (1)–(10)) corresponds to a DAE system of index two (Pantelides, Gritsis, Morrison, & Sargent, 1988). There is no generalpurpose stable algorithm to solve high index DAEs, and before any algorithm integrate the equations, a set of consistent initial conditions must be defined, hence not all differential variables may be freely chosen (Ascher & Petzold, 1998). Some state of the art numerical integration packages can solve DAE systems (DDASSL, gPROMS, MATLAB, for example), although this task is more difficult and numerically more expensive than solving ordinary differential equations (ODE). However, the model above can be transformed into a much simpler ODE system using the DDD approach proposed by Bosley and Edgar (1994). In this method, Table 2 Reaction kinetic parameters for linalool degradation Arrhenius kinetic expression r ¼ k 0 eEA =RT ðcÞS Free linalool Bounded linalool
k0
EA =R
S
1.205 0.058
12120.56 11000.0
2.052 1.540
ð11Þ
The whole simulation procedure is summarized in Fig. 3. A Gauss–Newton routine (FSOLVE) from Matlab 5.3 was used to pre-solve the set of non-linear algebraic equations, represented by the VLE model (Eq. (4)) and the partial condenser model (Eq. (10)). Two approaches were assessed for fitting the solution of the algebraic equations: polynomials and neural nets (NN). Polynomials can fit non-linear expressions with high accuracy if the prediction range is narrow. Bosley and Edgar (1994) originally proposed the use of polynomials, pðxÞ, of arbitrary degree in a quadratic form (Eq. (12)) where the vector /ðxÞ includes an arbitrary number of powers of the input variables, and W is a symmetrical matrix with coefficients that must be calibrated. pðxÞ ¼ /T ðxÞW/ðxÞ
ð12Þ
The main drawback of this approach is that polynomials quickly become intractable and difficult to calibrate if the system has many input variables. This is particularly critical in multicomponent distillation models. Artificial NN are a good alternative to polynomial fitting because they can accurately fit a wide range of non-linear functions with several input and output variables (Sartori & Panos, 1991). Sharma, Singhal, Ghosh, and Dwivedi (1999) have already demonstrated
D. Osorio et al. / Food Control 15 (2004) 515–521
519
Fig. 3. Simulation procedure diagram: algebraic equations systems are solved off-line and then replaced by an equivalent empirical function for the integration of differential equations.
the usefulness of NN as accurate and effortless computing tools for modeling VLE of highly polar multicomponent mixtures, where standard thermodynamic methods usually fail. Sigmoid activation functions and a back-propagation training algorithm were used in our work. The number of parameters in each model and the size of data sets used in model calibration are shown in Table 3. The prediction range of the polynomial model of the partial condenser was reduced in order to attain a reasonable fit. The drawback of both approaches is that training must be repeated when the process they represent is modified or if new experimental data becomes available. Table 3 Number of parameters and size of calibration data sets used in VLE and partial condenser models Polynomial fit
Artificial NN
VLE model (1 input)
Parameters Data
14 120
16 (5 neurons) 120
Partial condenser (3 inputs)
Parameters Data
50 750
55 (18 neurons) 1400
4. Results Three simulation strategies for wine distillation were compared in terms of computing efficiency and accuracy. In the first strategy, both VLE and partial condenser models were replaced by polynomials, while in the second strategy NN replaced these models. The third method solves the DAE distillation model directly, written in the matrix form W z0 ¼ f ðt; zÞ
ð13Þ
where empty rows in the diagonal mass matrix W cancel derivative terms to represent algebraic equations. The same stiff integration routine, ODE15S of MATLAB 5.3 based on numerical differentiation formulas (Shampine, Reichelt, & Kierzenka, 1999), and the same convergence criteria were used for the three simulation procedures. A 10 h distillation period under constant operating variables was used as the basis of the simulation. Physicochemical properties and initial concentrations of compounds in the mixture are listed in Table 4.
520
D. Osorio et al. / Food Control 15 (2004) 515–521
Table 4 Thermodynamic properties and initial concentration of mixture components Compound
Methanol Ethyl hexanoate Octanoic acid Linalool Degraded linalool Bounded linalool Ethanol Water
Concentration (mol/mol) 4
1.16 · 10 1.3 · 106 5.0 · 107 1.0 · 105 0 4.0 · 105 0.04 0.95983
Antoine Eqn. Parameters B log10 ðP Þ ¼ A T þC
UNIFAC activity coefficient P log10 ðcÞ ¼ 4k¼0 gk ðxeth Þk
A (mmHg)
B (mmHg Æ C)
C (C)
g0
8.07 7.00 7.39 7.20 7.20 – 8.21 5.08
1574.99 1496.34 1825.95 1629.57 1629.57 – 1652.05 1657.46
238.8 195.7 165.8 178.8 178.8 – 231.4 227.0
0.3591 )2.884 5.288 3.6519 )10.103 16.748 2.9444 )10.053 17.365 3.6717 )12.523 21.377 3.6717 )12.523 21.377 – – – K ¼ 1:689 (Van Laar) K ¼ 0:95174 (Van Laar)
4.1. Computing efficiency Each simulation strategy was assessed for an increasing number of trays to evaluate the computing efficiency at different problem scales. Fig. 4 illustrates the
g1
g2
g3
g4
)5.169 )5.151 )16.101 )9.347 )9.347 –
2.469 5.530 6.010 7.117 7.117 –
number of floating operations as a function of number of trays for the three simulation strategies. Our results show that the NN method was the most efficient, independently of the dimension of the equations system. The NN method required slightly fewer floating-point operations than the polynomial method, but only 58% of those demanded for solving the DAE distillation model, mainly due to the reduction in the size of the Jacobian. Such differences become particularly important when the simulator is used for control or optimization purposes (as is the case of this model), where a large number of simulation runs are required. 4.2. Computer accuracy
Fig. 4. Computing efficiency of assessed methods for a variable column size.
In terms of accuracy, the polynomial method presented the worst performance. The results shown in Fig. 5 depict the simulated evolution of the distillate composition in a six-tray column, setting Murphy efficiency
Fig. 5. Evolution of main components concentration in the distillate, as predicted by three different simulation methods.
D. Osorio et al. / Food Control 15 (2004) 515–521
to 80% for each tray. The evolution predicted by the polynomial method differs considerably from that given by the rigorous solution of the DAE system. By contrast, the NN approach followed the DAE simulation closely. Those differences are due to the limitation of polynomials to fit multivariable functions accurately. The moderate error introduced in the partial condenser model by the polynomial method has severe cumulative consequences in the accuracy of the dynamic simulation. Even after a polynomial with 50 coefficients was fitted for the partial condenser model, the mean quadratic error of the polynomial fitting was two orders of magnitude higher than it was for the NN model. The results indicate that it is highly convenient to combine neural networks along with DDD approach, as this simulation strategy maintains the accuracy while greatly reducing the computational time required.
5. Conclusions Our results showed that the neural networks strategy is more accurate than the polynomial fitting strategy, especially when dealing with multicomponent systems. In addition, the proposed neural network strategy is 40% more efficient than the DAE solving procedure compared, while keeping the same level of accuracy. The proposed model could be applied to the distillation of other spirits, as well as in other separation processes where aroma recovery is relevant, such as in juices, perfumes and throughout the food industry. Current work is considering the application of this model to the design of optimal operation policies in wine distillation, to maximize both the production and quality of Pisco brandy.
Acknowledgements The authors thank Lenka Torres for her assistance in laboratory work and Alex Crawford for editorial assistance.
References Alejski, K., & Dupart, F. (1996). Dynamic simulation of the multicomponent reactive distillation. Chemical Engineering Science, 51(18), 4237–4253.
521
Ascher, U. M., & Petzold, L. R. (1998). Computer methods for ordinary differential equations and differential-algebraic equations. Philadelphia: SIAM. Baumes, R., Bayonove, C., & Gunata, Z. (1994). Connaissances actuelles sur le potentiel aromatique des Muscats. Progres Agric et Vitic, 111(11), 251–256. Berna, A., Bon, J., Sanjuan, N., & Mulet, A. (1995). Concentraci on de aromas del vino: estudio mediante un simulador general (ProSim). Food Science and Technology International, 1, 117–127. Bosley, J. R., & Edgar, T. (1994). An efficient dynamic model for batch distillation. Journal of Process Control, 4(4), 195–204. Conner, J., Brikmyre, L., Patterson, A., & Piggott, J. (1998). Headspace concentrations of ethyl esters at different alcoholic strengths. Journal of the Science of Food and Agriculture, 77, 121–126. Franks, F. (1972). Water a comprehensive treatise (vol. 1). London: Plenum Press. Hansen, H. K., Rasmussen, P., Fredenslund, A., Schiller, M., & Gmehling, J. (1991). Vapor liquid equilibria by UNIFAC group contribution, 5. Revision and extension. Industrial Engineering and Chemical Research, 30, 2355–2358. Lora, J., Iborra, M., Perez, R., & Carbonell, I. (1992). Revista Espa~nola de Ciencia y Tecnologı a de Alimentos, 32(6), 621–633. Loyola, E., Almy, J., Martın-Alvarez, P., & Cabezudo, M. D. (1987). The flavor of Chilean Pisco. In Proceedings of the fifth international flavor conference, Porto Karras, Greece (pp. 729–742). Murphy, J. A., & Gaines, G. L., Jr (1974). Density and viscosity of aqueous hydrogen sulfide solutions at pressures to 20 atm. Journal of Chemical and Engineering Data, 19(4), 359–362. Neuburg, H. J., & Perez-Correa, J. R. (1994). Dynamic and steady state modelling of a pilot binary tray distillation column. Latin American Applied Research, 24, 1–15. Ohta, T., (1995). Transformations of geraniol, nerol and their gluco sides during Schochu distillation. In R. Cantagrel (Ed.), Elaboration et Connaissance des Spirituex. Vient de Paraitre sous l’egide du Bureau National Interprofessionel du Cognac (pp. 313–315). Pantelides, C. C., Gritsis, D., Morrison, K. R., & Sargent, R. W. (1988). The mathematical modeling of transient systems using differential and algebraic equations. Computers and Chemical Engineering, 12, 449. Perry, R. H., Green, D. W., & Maloney, J. O. (1999). Perry’s chemical engineers’ handbook (7th ed). New York: McGraw-Hill. Prauznitz, J. M., Lichtenthaler, R. N., & Gomes de Azevedo, E. (1999). Molecular thermodynamics of fluid––phase equlibria. New York: Prentice-Hall. Ruiz, C., Basualdo, M., & Scenna, N. (1995). Reactive distillation dynamic simulation. Transactions in Chemical Engineering A, 73, 363–378. Sartori, M. A., & Panos, J. A. (1991). A simple method to derive bounds on the size and to train multilayer neural networks. IEEE Transactions on Neural Network, 2(4), 467–471. Shampine, L. F., Reichelt, M. W., & Kierzenka, J. A. (1999). Solving Index-1 DAEs in MATLAB and Simulink. SIAM Review, 41, 538– 552. Schneider, R., Noeres, L., Kreul, L., & G orak, A. (2001). Dynamic modeling and simulation of reactive batch distillation. Computers and Chemical Engineering, 25, 169–176. Sharma, R., Singhal, D., Ghosh, R., & Dwivedi, A. (1999). Potential applications of artificial neural networks to thermodynamics: vapor–liquid equilibrium predictions. Computers and Chemical Engineering, 23, 385–390.
Rigorous dynamic modeling and simulation of wine distillations Daniel Osorio, Ricardo Perez-Correa *, Andrea Belancic, Eduardo Agosin Department of Chemical and Bioprocess Engineering, P. Universidad Cat olica de Chile, P.O. Box 306––22, Santiago, Chile Received 23 November 2002; received in revised form 11 August 2003; accepted 11 August 2003
Abstract A novel simulation strategy for dynamic distillation of complex mixtures, such as wine, is proposed and evaluated in terms of computing efficiency and accuracy. The model developed describes wine distillation as a multicomponent reactive batch distillation process. The simulation approach transforms the system of differential algebraic equations (DAE) into a set of ordinary differential equations, by pre-solving the algebraic equations and replacing them with artificial neural networks. This new simulation strategy for wine distillation is 40% faster than the rigorous solution of the DAE system, compared at the same level of accuracy. The model can be applied to the distillation of other spirits or complex mixtures, as well as in other separation processes in which the recovery of aromas is essential. 2003 Elsevier Ltd. All rights reserved. Keywords: Reactive-batch distillation; Terpenes; Pisco
1. Introduction Distilling wine for the production of spirits is an early application of reactive batch distillation. Nowadays, challenges regarding this highly non-linear process, as design improvements, automatic control and optimization, all demand efficient and accurate simulation tools. Although various rigorous batch distillation models are available (Alejski & Dupart, 1996; Ruiz, Basualdo, & Scenna, 1995; Schneider, Noeres, Kreul, & G orak, 2001), solving these models requires the application of special purpose differential-algebraic equation (DAE) solvers, which are computationally intensive and may present initialization problems. Published work on the modeling of wine distillation is scarce and mainly focuses on the much simpler continuous system (Berna, Bon, Sanju an, & Mulet, 1995; Lora, Iborra, Perez, & Carbonell, 1992). Bosley and Edgar (1994) introduced a functional approach for batch distillation modeling called direct dynamic derivatives (DDD), which avoids the need for sophisticated and expensive DAE solvers. It was
*
Corresponding author. Tel.: +56-582-686-4258; fax: +56-562-6865803. E-mail address: [email protected] (R. Perez-Correa). 0956-7135/$ - see front matter 2003 Elsevier Ltd. All rights reserved. doi:10.1016/j.foodcont.2003.08.003
employed successfully to model water/ethanol batch distillation. However, the simple vapor/liquid equilibrium model and polynomial fitting they used is limited and unsuitable for simulating multicomponent distillation. In this work we explore the use of neural networks to improve the accuracy, performance and range of applicability of the DDD approach, making it suitable for wine distillation. A simplified multicomponent mixture was used representing wine.
2. Modeling wine distillation process Thermodynamically, wine is a highly non-ideal colloidal suspension that contains several hundred components. Its main constituents are water and ethanol, although the many volatile aromas that are found in remarkably small proportions define the quality of the wine and its distillate. Among these trace components, which we denote here as ‘‘minor components’’, the main positive quality factors in wine distillates are esters, terpenes and aromatic alcohols. In turn, methanol, fatty acids, C6-components and superior alcohols are the most deleterious to distillate quality. During the distillation process, volatile components are extracted from wine and concentrated in the spirit, with the aim to recover the maximum amount of ethanol and positive
516
D. Osorio et al. / Food Control 15 (2004) 515–521
Nomenclature Notations a heat transfer area, m2 c concentration, mg/l Cp heat capacity, J/C kg dw weir diameter, cm EA energy of activation EM Murphy efficiency FR coolant flowrate, l/min H molar enthalpy of vapor, J/mol h molar enthalpy of liquid, J/mol hw weir height, cm k0 Arrhenius constant kw Francis weir constant, mol/s cm5=2 L liquid flowrate, mol/s l liquid level on the tray, cm M molar holdup, mol MW molecular weight, g/mol NC number of components P total pressure, mmHg P0 pure component vapor pressure, mmHg Q heat transfer rate, J/s R universal gas constant, J/mol K r net reaction rate, mol/s S reaction order T temperature, C
characteristic aromas, while minimizing the off-flavors in the distillate. As a case study, this work focused on the distillation of Pisco brandy, a wine distillate from Chile and Peru with a distinctive Muscat and fruity aroma, whose character is mainly attributable to its high terpene content (Loyola, Almy, Martın-Alvarez, & Cabezudo, 1987). In Muscat wines, terpenes are found either free or bound to sugars as glycosides (Baumes, Bayonove, & Gunata, 1994). The latter are termed aroma precursors, as they are unable to express their aromatic character, although when hydrolysed during distillation these precursors do constitute a further source of free terpenes (Ohta, 1995). Distillation temperature and low pH conditions however oxidize terpenes, driving them to less positive or negative odorant forms. The behaviour of linalool, the most important terpene found in Pisco, is analysed in detail in this work.
t U V W x y z
time, s global heat transfer coefficient, J/s m2 vapor flowrate, mol/s polynomial coefficients matrix liquid phase composition, mol/mol vapor phase composition, mol/mol state variables vector
Greek characters c activity coefficient q density, g/ml U apparent liquid volume, ml / polynomial fit terms W mass matrix Subindex n tray i component eth ethanol out outlet partial condenser coolant stream in inlet partial condenser coolant stream Superscripts * ideal T transpose
• Constant pressure through the column; • Perfect mixing occurs in each tray; • Vapor phase accumulation is neglected. In wine distillates, water and ethanol account for over 96% of the total mass. Minor components, while key to spirit quality, have an insignificant influence upon the thermodynamic properties of the mixture. Hence, distillate density, enthalpy, molecular weight and boiling temperature were obtained from water–ethanol binary expressions. Furthermore, the heat of reaction and the heat of mixing were also neglected. This case study considers a simplified mixture including eight of the most relevant components in the aromatic profile of Pisco brandy: ethanol, water, methanol, octanoic acid, ethyl hexanoate and free, bounded and degraded linalool. 2.2. Fundamental equations
2.1. Model assumptions
Total mass, component, and energy balances (Eqs. (1)–(3)) are presented for a generic tray n.
Given the small size of the distillation column generally used in this process (6–10 trays, 0.4 m diameter), the following assumptions are considered applicable:
C X dMn ¼ Vnþ1 þ Ln1 ðVn þ Ln Þ þ rn;i dt i¼1
N
ð1Þ
D. Osorio et al. / Food Control 15 (2004) 515–521
Mn
Table 1 Parameters for the Francis hydraulic equation
dxn;i ¼ Ln1 ðxn1;i xn;i Þ Vn ðyn;i xn;i Þ dt NC X þ Vnþ1 ðynþ1;i xn;i Þ þ rn;i xn;i rn;i ;
Francis weir expression Ln ¼ kw dw ðln hw Þ1;5
i¼1
8i ¼ 1; . . . ; NC 1 Mn
ð2Þ
dhn ¼ Ln1 ðhn1 hn Þ Vn ðHn hn Þ dt NC X rn;i þ Vnþ1 ðHnþ1 hn Þ þ Qn hn
ð3Þ qn ¼
Since only NC 1 balances are independent, water balance was not considered. The net reaction rate for component i, rn;i , is the difference between generation and consumption rates. The heat transfer term, Qn , is positive for the reboiler and negative for the partial condenser while it is neglected for intermediate trays (the column is insulated).
The Raoult law was applied to find the theoretical equilibrium composition of the gas phase, yn;i , using Antoine expression to obtain vapor pressures of pure components, P 0 . As the sum of the gas fractions on each tray must be equal to one, the Raoult equation and the consistency equation can be combined into the following single algebraic equation: ð4Þ
Experimental and theoretical studies have shown that activity coefficients of organic minor components in water/ethanol mixtures depend on the ethanol content of the mixture only (Conner, Brikmyre, Patterson, & Piggott, 1998). The UNIFAC method (Hansen, Rasmussen, Fredenslund, Schiller, & Gmehling, 1991) was used to estimate the activity coefficients of minor components in terms of ethanol composition, not available in literature. Activity coefficients for water and ethanol were estimated using the more accurate Van Laar expression (Prauznitz, Lichtenthaler, & Gomes de Azevedo, 1999). The actual composition in the gas phase, yn , was calculated introducing the global efficiency Murphy factor, EM . yn;i ¼ ynþ1;i þ EM ðyn;i ynþ1;i Þ
6 32 3
xn;eth MWeth þ ð1 xn;eth ÞMWw 103 Un xn;eth þ ð1 xn;eth ÞMWw =qw
ð6Þ
where Un is the apparent liquid molar volume, that for water/ethanol mixture was obtained from the following expression (Franks, 1972) Un ¼ 5:1214 102 þ 6:549 103 xn;eth þ 7:406 105 Tn
ð7Þ
Enthalpy expressions for the water/ethanol mixture were taken from Neuburg and Perez-Correa (1994)
2.3. Vapor–liquid equilibrium
NC 1 X 0 c xn;i Pn;i ¼0 Pn i¼1 n;i
kw (mol/s cm5=2 ) dw (cm) hw (cm)
(1974) to calculate volumetric accumulation required in Francis equation,
i¼1
1
517
ð5Þ
2.4. Constitutive equations Liquid flowrates were calculated using Francis’ weir expression (Perry, Green, & Maloney, 1999), using parameters shown in Table 1. Liquid density was modeled using the method recommended by Murphy and Gaines
Hn ¼ 44765:71 6171:03yn;eth þ ð31:46 11:98yn;eth Þ Tn þ ð4:063 104 þ 0:073yn;et Þ Tn
ð8Þ
hn ¼ ð55:678 xn;eth þ 75:425Þ Tn 0:0057 xn;eth 0:00125
ð9Þ
2.5. Partial condenser Reflux ratio is generated by partially condensing the vapor outlet stream at the top of the column using a heat exchanger device. The heat transfer in this partial condenser (tray number 0) is assumed at pseudo-steadystate, therefore the energy accumulation term is neglected. As a result, the energy balance applied to the coolant stream is reduced to an algebraic equation that must be solved iteratively for Tout to obtain the partial condenser heat absorption. Q0 FR Cp ðTout Tin Þ ðTout Tin Þ ¼ aU lnððT0 Tin Þ=ðT0 Tout ÞÞ
ð10Þ
2.6. Reaction kinetics For the reactive components, in this case bounded, free and degraded linalool, a simple chain reaction model was employed. It is assumed that reactions take place in the liquid phase in each tray. Reaction rates are calculated using an experimentally fitted Arrhenius kinetic expression obtained in this work (see Figs. 1 and 2). Parameter values are shown in Table 2.
518
D. Osorio et al. / Food Control 15 (2004) 515–521
Fig. 1. Free linalool reaction kinetic at different temperatures.
Fig. 2. Bounded linalool reaction kinetic at different temperatures.
the implicit non-linear algebraic equations (Eqs. (4) and (10)) are pre-solved using standard non-linear optimization routines, and then polynomials are fitted to their solution. This avoids solving the algebraic system repeatedly at every step of the integration procedure. Moreover, it eliminates the requirement of consistent initial conditions, an inconvenient of DAE formulation that becomes harder to overcome as the number of variables increases. According to the vapor-liquid equilibrium (VLE) equations, the mixture composition depends only on equilibrium temperature and pressure. Consequently, the system of VLE equations may be solved independently for temperature for a wide range of composition and pressure values, and then an empirical non-linear expression can be fitted to the solution. As a result, the equilibrium temperature and enthalpy expressions can be evaluated explicitly in terms of tray pressure and composition using the adjusted function. The replacement allows the direct evaluation of internal vapor flowrates, since temperature derivatives are eliminated from the energy balance (see Bosley & Edgar, 1994, for details). Therefore, and grouping composition as vectors, (i.e. xn ¼ fx1 ; x2 ; . . . ; xn g),
oh T oh T PNC Ln1 hn1 hn ox ðxn1 xn Þ þ Vnþ1 Hnþ1 hn ox ðynþ1 xn Þ þ Qn hn i¼1 rn;i Vn ¼ oh T Hn hn ox ðyn xn Þ
3. Simulation The model described above (Eqs. (1)–(10)) corresponds to a DAE system of index two (Pantelides, Gritsis, Morrison, & Sargent, 1988). There is no generalpurpose stable algorithm to solve high index DAEs, and before any algorithm integrate the equations, a set of consistent initial conditions must be defined, hence not all differential variables may be freely chosen (Ascher & Petzold, 1998). Some state of the art numerical integration packages can solve DAE systems (DDASSL, gPROMS, MATLAB, for example), although this task is more difficult and numerically more expensive than solving ordinary differential equations (ODE). However, the model above can be transformed into a much simpler ODE system using the DDD approach proposed by Bosley and Edgar (1994). In this method, Table 2 Reaction kinetic parameters for linalool degradation Arrhenius kinetic expression r ¼ k 0 eEA =RT ðcÞS Free linalool Bounded linalool
k0
EA =R
S
1.205 0.058
12120.56 11000.0
2.052 1.540
ð11Þ
The whole simulation procedure is summarized in Fig. 3. A Gauss–Newton routine (FSOLVE) from Matlab 5.3 was used to pre-solve the set of non-linear algebraic equations, represented by the VLE model (Eq. (4)) and the partial condenser model (Eq. (10)). Two approaches were assessed for fitting the solution of the algebraic equations: polynomials and neural nets (NN). Polynomials can fit non-linear expressions with high accuracy if the prediction range is narrow. Bosley and Edgar (1994) originally proposed the use of polynomials, pðxÞ, of arbitrary degree in a quadratic form (Eq. (12)) where the vector /ðxÞ includes an arbitrary number of powers of the input variables, and W is a symmetrical matrix with coefficients that must be calibrated. pðxÞ ¼ /T ðxÞW/ðxÞ
ð12Þ
The main drawback of this approach is that polynomials quickly become intractable and difficult to calibrate if the system has many input variables. This is particularly critical in multicomponent distillation models. Artificial NN are a good alternative to polynomial fitting because they can accurately fit a wide range of non-linear functions with several input and output variables (Sartori & Panos, 1991). Sharma, Singhal, Ghosh, and Dwivedi (1999) have already demonstrated
D. Osorio et al. / Food Control 15 (2004) 515–521
519
Fig. 3. Simulation procedure diagram: algebraic equations systems are solved off-line and then replaced by an equivalent empirical function for the integration of differential equations.
the usefulness of NN as accurate and effortless computing tools for modeling VLE of highly polar multicomponent mixtures, where standard thermodynamic methods usually fail. Sigmoid activation functions and a back-propagation training algorithm were used in our work. The number of parameters in each model and the size of data sets used in model calibration are shown in Table 3. The prediction range of the polynomial model of the partial condenser was reduced in order to attain a reasonable fit. The drawback of both approaches is that training must be repeated when the process they represent is modified or if new experimental data becomes available. Table 3 Number of parameters and size of calibration data sets used in VLE and partial condenser models Polynomial fit
Artificial NN
VLE model (1 input)
Parameters Data
14 120
16 (5 neurons) 120
Partial condenser (3 inputs)
Parameters Data
50 750
55 (18 neurons) 1400
4. Results Three simulation strategies for wine distillation were compared in terms of computing efficiency and accuracy. In the first strategy, both VLE and partial condenser models were replaced by polynomials, while in the second strategy NN replaced these models. The third method solves the DAE distillation model directly, written in the matrix form W z0 ¼ f ðt; zÞ
ð13Þ
where empty rows in the diagonal mass matrix W cancel derivative terms to represent algebraic equations. The same stiff integration routine, ODE15S of MATLAB 5.3 based on numerical differentiation formulas (Shampine, Reichelt, & Kierzenka, 1999), and the same convergence criteria were used for the three simulation procedures. A 10 h distillation period under constant operating variables was used as the basis of the simulation. Physicochemical properties and initial concentrations of compounds in the mixture are listed in Table 4.
520
D. Osorio et al. / Food Control 15 (2004) 515–521
Table 4 Thermodynamic properties and initial concentration of mixture components Compound
Methanol Ethyl hexanoate Octanoic acid Linalool Degraded linalool Bounded linalool Ethanol Water
Concentration (mol/mol) 4
1.16 · 10 1.3 · 106 5.0 · 107 1.0 · 105 0 4.0 · 105 0.04 0.95983
Antoine Eqn. Parameters B log10 ðP Þ ¼ A T þC
UNIFAC activity coefficient P log10 ðcÞ ¼ 4k¼0 gk ðxeth Þk
A (mmHg)
B (mmHg Æ C)
C (C)
g0
8.07 7.00 7.39 7.20 7.20 – 8.21 5.08
1574.99 1496.34 1825.95 1629.57 1629.57 – 1652.05 1657.46
238.8 195.7 165.8 178.8 178.8 – 231.4 227.0
0.3591 )2.884 5.288 3.6519 )10.103 16.748 2.9444 )10.053 17.365 3.6717 )12.523 21.377 3.6717 )12.523 21.377 – – – K ¼ 1:689 (Van Laar) K ¼ 0:95174 (Van Laar)
4.1. Computing efficiency Each simulation strategy was assessed for an increasing number of trays to evaluate the computing efficiency at different problem scales. Fig. 4 illustrates the
g1
g2
g3
g4
)5.169 )5.151 )16.101 )9.347 )9.347 –
2.469 5.530 6.010 7.117 7.117 –
number of floating operations as a function of number of trays for the three simulation strategies. Our results show that the NN method was the most efficient, independently of the dimension of the equations system. The NN method required slightly fewer floating-point operations than the polynomial method, but only 58% of those demanded for solving the DAE distillation model, mainly due to the reduction in the size of the Jacobian. Such differences become particularly important when the simulator is used for control or optimization purposes (as is the case of this model), where a large number of simulation runs are required. 4.2. Computer accuracy
Fig. 4. Computing efficiency of assessed methods for a variable column size.
In terms of accuracy, the polynomial method presented the worst performance. The results shown in Fig. 5 depict the simulated evolution of the distillate composition in a six-tray column, setting Murphy efficiency
Fig. 5. Evolution of main components concentration in the distillate, as predicted by three different simulation methods.
D. Osorio et al. / Food Control 15 (2004) 515–521
to 80% for each tray. The evolution predicted by the polynomial method differs considerably from that given by the rigorous solution of the DAE system. By contrast, the NN approach followed the DAE simulation closely. Those differences are due to the limitation of polynomials to fit multivariable functions accurately. The moderate error introduced in the partial condenser model by the polynomial method has severe cumulative consequences in the accuracy of the dynamic simulation. Even after a polynomial with 50 coefficients was fitted for the partial condenser model, the mean quadratic error of the polynomial fitting was two orders of magnitude higher than it was for the NN model. The results indicate that it is highly convenient to combine neural networks along with DDD approach, as this simulation strategy maintains the accuracy while greatly reducing the computational time required.
5. Conclusions Our results showed that the neural networks strategy is more accurate than the polynomial fitting strategy, especially when dealing with multicomponent systems. In addition, the proposed neural network strategy is 40% more efficient than the DAE solving procedure compared, while keeping the same level of accuracy. The proposed model could be applied to the distillation of other spirits, as well as in other separation processes where aroma recovery is relevant, such as in juices, perfumes and throughout the food industry. Current work is considering the application of this model to the design of optimal operation policies in wine distillation, to maximize both the production and quality of Pisco brandy.
Acknowledgements The authors thank Lenka Torres for her assistance in laboratory work and Alex Crawford for editorial assistance.
References Alejski, K., & Dupart, F. (1996). Dynamic simulation of the multicomponent reactive distillation. Chemical Engineering Science, 51(18), 4237–4253.
521
Ascher, U. M., & Petzold, L. R. (1998). Computer methods for ordinary differential equations and differential-algebraic equations. Philadelphia: SIAM. Baumes, R., Bayonove, C., & Gunata, Z. (1994). Connaissances actuelles sur le potentiel aromatique des Muscats. Progres Agric et Vitic, 111(11), 251–256. Berna, A., Bon, J., Sanjuan, N., & Mulet, A. (1995). Concentraci on de aromas del vino: estudio mediante un simulador general (ProSim). Food Science and Technology International, 1, 117–127. Bosley, J. R., & Edgar, T. (1994). An efficient dynamic model for batch distillation. Journal of Process Control, 4(4), 195–204. Conner, J., Brikmyre, L., Patterson, A., & Piggott, J. (1998). Headspace concentrations of ethyl esters at different alcoholic strengths. Journal of the Science of Food and Agriculture, 77, 121–126. Franks, F. (1972). Water a comprehensive treatise (vol. 1). London: Plenum Press. Hansen, H. K., Rasmussen, P., Fredenslund, A., Schiller, M., & Gmehling, J. (1991). Vapor liquid equilibria by UNIFAC group contribution, 5. Revision and extension. Industrial Engineering and Chemical Research, 30, 2355–2358. Lora, J., Iborra, M., Perez, R., & Carbonell, I. (1992). Revista Espa~nola de Ciencia y Tecnologı a de Alimentos, 32(6), 621–633. Loyola, E., Almy, J., Martın-Alvarez, P., & Cabezudo, M. D. (1987). The flavor of Chilean Pisco. In Proceedings of the fifth international flavor conference, Porto Karras, Greece (pp. 729–742). Murphy, J. A., & Gaines, G. L., Jr (1974). Density and viscosity of aqueous hydrogen sulfide solutions at pressures to 20 atm. Journal of Chemical and Engineering Data, 19(4), 359–362. Neuburg, H. J., & Perez-Correa, J. R. (1994). Dynamic and steady state modelling of a pilot binary tray distillation column. Latin American Applied Research, 24, 1–15. Ohta, T., (1995). Transformations of geraniol, nerol and their gluco sides during Schochu distillation. In R. Cantagrel (Ed.), Elaboration et Connaissance des Spirituex. Vient de Paraitre sous l’egide du Bureau National Interprofessionel du Cognac (pp. 313–315). Pantelides, C. C., Gritsis, D., Morrison, K. R., & Sargent, R. W. (1988). The mathematical modeling of transient systems using differential and algebraic equations. Computers and Chemical Engineering, 12, 449. Perry, R. H., Green, D. W., & Maloney, J. O. (1999). Perry’s chemical engineers’ handbook (7th ed). New York: McGraw-Hill. Prauznitz, J. M., Lichtenthaler, R. N., & Gomes de Azevedo, E. (1999). Molecular thermodynamics of fluid––phase equlibria. New York: Prentice-Hall. Ruiz, C., Basualdo, M., & Scenna, N. (1995). Reactive distillation dynamic simulation. Transactions in Chemical Engineering A, 73, 363–378. Sartori, M. A., & Panos, J. A. (1991). A simple method to derive bounds on the size and to train multilayer neural networks. IEEE Transactions on Neural Network, 2(4), 467–471. Shampine, L. F., Reichelt, M. W., & Kierzenka, J. A. (1999). Solving Index-1 DAEs in MATLAB and Simulink. SIAM Review, 41, 538– 552. Schneider, R., Noeres, L., Kreul, L., & G orak, A. (2001). Dynamic modeling and simulation of reactive batch distillation. Computers and Chemical Engineering, 25, 169–176. Sharma, R., Singhal, D., Ghosh, R., & Dwivedi, A. (1999). Potential applications of artificial neural networks to thermodynamics: vapor–liquid equilibrium predictions. Computers and Chemical Engineering, 23, 385–390.