CHEMICAL ENGINEERING TRANSACTIONS VOL. 57, 2017 A publication of The Italian Association of Chemical Engineering Online at www.aidic.it/cet Guest Editors: Sauro Pierucci, Jiří Jaromír Klemeš, Laura Piazza, Serafim Bakalis Copyright © 2017, AIDIC Servizi S.r.l. ISBN 978-88-95608- 48-8; ISSN 2283-9216 Application of Homotopy Continuation Method for Robust Simulation of Distillation Columns Fernando C. Zumach, Reginaldo Guirardello School of Chemical Engineering, State University of Campinas (UNICAMP), Av. Albert Einstein 500, 13083-852, Campinas- SP, Brazil fernandocassoli@hotmail.com The production method most commonly used for obtaining bioethanol is the fermentation of sugarcane, which has been proved as one of the best raw material for the ethanol production. The alcohol produced from this process must be distilled to attain the minimum specifications required by the regulations. The hydrated ethanol fuel must be concentrated between 95.1% (v/v) and 96.0% (v/v) to be commercialized, and distillation is the most frequently used method to carry out this separation process. Despite distillation being the most frequently used method to carry that separation, this process has a high energy consumption, in the separation process due to the excess water. Also, the modeling of the distillation column sometimes has numerical difficulties, due to convergence problems in solving the nonlinear equations. Therefore, this work has as its objective the development of a rigorous model of a distillation column for the production of the anhydrous ethanol. The models generated in this study are implemented in the modeling language EMSO simulator; a free simulator that has full graphical environment and allows users to work with both dynamic and stationary processes. Given the difficulty in the numerical resolution of the nonlinear equations, in this work it will be applied the method of Homotopic Continuation. This method consists in creating an easy resolution system and gradually getting better initial estimates by varying a parameter that continuously transforms the easy problem in the original problem. In such way, the method has better convergence properties and is globally convergent. 1. Introduction In view of the global concern of energy security, new policies were created around the world for the creation of renewable energy sources. Thus, many measures were implemented; such as the addition of biofuels to fossil fuels. Due to that the attention given to the production of biodiesel and bioethanol as renewable fuel sources, in order to replace or supplement gas and diesel based fuels, has been exponentially increasing, (Zanette, 2010). The largest biofuel producers are the United States and Brazil. In 2011 the combined production of both countries represented about 87% of the world production of these fuels. The production of biofuels has grown greatly in recent years, in 2000 it was produced little less than 18 billion liters/year, however in 2011, this production reached 107 billion liters/year, therefore, an increase of approximately 89 billion liters in a few years (REN21, 2012). Currently, biofuels supply about 3% of the total fuel consumed in transport in the world; these levels can vary greatly in some countries. In Brazil, for example, 21% of fuels used in transportation are biofuels. In other countries, like the United States, this percentage is 4%, and the European Union is 3% (IEA, 2011). The technology to obtain ethanol is one of the most important scientific fields in the development of biofuels. The production technique depends on the type of raw material used. Ethanol can be manufactured using any raw material, which has sugar, together with the action of bacteria and yeasts, or it may be obtained from polysaccharides, such as starch and cellulose, which are broken down into monosaccharides, to then be metabolized to ethanol (Dodic et al., 2012). The main raw material used in the ethanol manufacturing process in Brazil is sugarcane. This occurs due to the fact that sugar cane has several advantages, such as being highly efficient with a low cost, besides being DOI: 10.3303/CET1757178 Please cite this article as: Zumach F., Guirardello R., 2017, Application of homotopy continuation method for robust simulation of distillation columns, Chemical Engineering Transactions, 57, 1063-1068 DOI: 10.3303/CET1757178 1063 the raw material that has the easiest production method, since the fermentation process takes place directly, thus, other chemical or biological steps are not necessary (Goldemberg et al., 2008). Besides that, the sugarcane bagasse coming from the production of biofuels can be used in energy cogeneration. This makes it possible for the power plants to become self-sufficient in heat and electric energy, with the possibility of selling surplus power (Ensinas, 2008). The fermentation process of biomass containing conventional sugars and starch generates a must with ethanol content of about 10% (v/v) (Haelssig et al., 2011). The alcohol produced from that process must go through a purification method to meet the minimum specifications required by the regulations. According to Resolution No. 7 of the ANP (2011), the anhydrous fuel ethanol must have a minimum concentration of 99.6% (v/v), while the hydrated ethanol fuel should be between 95.1% (v/v) and 96.0% (v/v). Projecting a distillation tower can be a difficult task, since many operating parameters affect the process of separation. Furthermore, the construction of a pilot plant is costly and consumes much time. An alternative to the problem is to use computer simulation of the distillation process. Nowadays there is a large number of simulators on the market with the capacity of simulating a distillation tower for separation of ethanol-water; and in this study we used the EMSO. The EMSO is a free simulator that has full graphical environment and allows users to work with both dynamic and stationary processes. It is possible to simply select and connect the template blocks, or insert your own programming language. Moreover, the software already has a library with various templates that can be used (Soares, 2007). Regarding the difficulty of efficiently making the modeling of a distillation column with many interactions, in this project the Homotopy continuation method will be used. This method creates an easy system resolution through the convex combination with the original system. It is possible to get gradually better initial estimates by the variance of parameter of the combination. Thus, the Homotopy continuation method has better convergence properties. 2. Methodology The utilization of simulation software to portray the industrial chemical processes is a technique that has been increasing in its use, because it facilitates the choice of optimum parameters of the process, reducing costs and improving efficiency, among other factors. Currently in the market there are several commercial programs that can be used to perform the simulation, such as: II Pro, HYSYS, ASPEN e EMSO. At this study we used the EMSO (Environment for Modelling, Simulation and Optimization) to carry out simulations. For the calculation of the distillation column, it is necessary to measure the temperature, the flow rate and compositions of the currents and the heat transfer rate in each stage. To achieve this determination the resolution of the material balance, energy balance and the equilibrium relationships for each stage is required. In developing the mathematical model for the distillation column, the stage of balance principle is utilized. The expressions utilized to describe the column are: mass balance (M), equilibrium relationships, (E), summation of the molar fractions (S), and energy balances (H). The combination of these equations is called MESH equations. In the figure below, there is a schematic representation of a j tray. Figure 1: Schematic representation of a generic tray. 1064 Based on Figure 1, it is possible to obtain the following basic equations: Equation M - Mass balance for component i on the tray j: 0)()( ,,,1,11,,   jijjjijjjijjijjiijji yWVxULzFyVxLM (1) Equation E - Phase equilibrium relations: 0 ,,,,  jijijiji xKyE (2) Equation S - Sum of mole fractions:   c jijy yS 1 , 01)( (3)   c jijx xS 1 , 01)( (4) Equation (H) - Energy Balance 0)()( ,,,1,11,   jVjjjLjjjFjjVjjLijj HWVHULHFHVHLH (5) where, F feed, L liquid flow, V steam flow, HV steam enthalpy, HL liquid enthalpy, x mole fraction of component i in the liquid, y mole fraction of component i in the steam, W steam purge, U liquid purge, T temperature and P pressure. For the development of a distillation column model, several peripheral equipment models were designed. To accomplish the necessary modeling it is necessary to enumerate stages of the distillation column, being that: the condenser will be counted as first stage and the reboiler as the last. In order to obtain the enthalpy values, fugacity and other thermodynamic parameters during the simulation of the distillation process (EMSO), a thermodynamic package called VRTherm was installed. With the installation of this software, it becomes possible to obtain physical and thermodynamic properties of the compounds, assisting in modeling, simulation and optimization of processes. The VRTherm has in its library properties of pure substances and mixtures, besides using correlation models as: Redlich-Kwong, Peng-Robinson, Van der Waals, UNIFAC, among others (Smith et al., 2007). Based on the thermodynamic models available in VRTherm, the following models will be used to predict the VLE of the mixture: UNIFAC/SRK, (for the calculation of the liquid phase), and Redlich-Kwong (to calculate the steam phase). Separation of ethanol-water during the production of hydrous ethanol. The model of the distillation column was constructed using equations described in previous sections. The equations were inserted into the EMSO, using the model, which is one of the EMSO programming languages. In addition, the deviation equilibrium was introduced through the concept of efficiency in the equilibrium relations. Murphree’s efficiency relates a real stage with a ideal, considering the liquid is completely mixed. In fact, it can be defined either by gas phase or by liquid phase using the equations above: Gas phase: η𝑖,𝑗 𝑉 = 𝑦𝑖,𝑗− 𝑦𝑖,𝑗+1 𝑦𝑖,𝑗 ∗ − 𝑦𝑖,𝑗+1 (6) Liquid phase: η𝑖,𝑗 𝐿 = 𝑋𝑖,𝑗−1− 𝑋𝑖,𝑗 𝑋𝑖,𝑗−1 ∗ − 𝑋𝑖,𝑗 (7) Noted that: 𝛈 efficiency, x molar fraction of liquid phase, x* molar fraction of the liquid in equilibrium with y, y molar fraction of gas phase, y* molar fraction of gas phase in equilibrium with x. For the simulation, a constant value of efficiency of Murphree was initially considered for all stages equal to 1.According to the process in Figure 2, all units of the column, (the equilibrium stage, reboiler, condenser, etc.), were modeled. Therefore, a set of submodels was created that connected together compose the model of the entire separation system under analysis. 1065 Figure 2: Flowchart of the process of distillation in EMSO. Method was applied to increase the degree of convergence, thus, increasing the robustness of the simulation. Such result was achieved by, reformulating the MESH equations and adding the Homotopy parameters and auxiliary homotopic functions for the solution of nonlinear equations systems. The great importance of good initial estimates is to decrease the computational effort and mainly ensure the convergence of the simulation. The traditional method of resolution of the EMSO (Newton-Raphson) for the solution of nonlinear equations shows high efficiency of convergence speed. However, it requires a specific initial estimate for the iterative process. Good initial estimates can lead to the solution quickly; however, poor initial estimates and / or the existence of multiple solutions may lead to divergence. The Homotopy Continuation method can circumvent such problems and ensure convergence (Wu, 2005; Tanskanen & Pohjola, 2000). In the original MESH system of equations, the feed molar fraction is considered to vary continuously according to EQ (6), (7) and (8): Feed column: Design i .zz i  1,.....,1  Ni (6)     1 1 Design .1 N i in zz  Ni  (7) Design i i z d dz   (8) However, since in this work the modelling is based on the transient (MESH) equations, the continuation variable ( ) is made dependent with time, given by equation (9): t e . 1     (9) Where the value of  must be defined to an adequate value according to the problem. Thus, the ODE that must be solved is: Design. ).( i ti ze d dz      (10) Which can be rearranged to: ).( Design ii i zz d dz    (11) and must be solved together with the transient MESH equations. 1066 3. Results Initially, dynamic state simulations of ethanol-water mixtures distillation were implemented in EMSO. In addition, a simulation was performed in the steady state through the results of the dynamic simulation. First, simulations on the model without the Homotopy Continuation method were developed to validate such model; then, these results found were compared with those obtained from AspenPlus® V.7.3 simulator. In order to validate the use of this method in the model, coarse initial conditions were selected (Table 1). All simulations of this work were established using a computer with the following configuration: Intel core I7 CPU 3.4 GHZ 3.4 GHZ and 12 GB RAM. Table 1: Initial conditions for simulation of a conventional distillation process. EMSO EMSO (Homotopy) Number of stages 25 25 Feed Stage 6 6 RR 1.5 1.5 Bottom stream 162 kmol/h 162 kmol/h Pressure (P) 1 atm 1 atm Ef. Murphee 1 1 Feed stream (F) 300 kmol/h 300 kmol/h Feed pressure (Pf) 1 atm 1 atm Feed temperature (Tf) 350 K 350 K Composition (F) / water 0.75 0.75 Composition (F) / ethanol 0.25 0.25 Based on the initial conditions from Table 1 in the two models (simulation time of 10000 seconds), the following results were verified (Table 2). In the model without the method, convergence failure. Nevertheless, the model with the Continuation method achieved good results. The model with the Homotopy Continuation method uses a larger computational time to get answers; however, only this model converged. Table 2: Results of conventional simulation. EMSO EMSO (Homotopy) Top stream (water) - 0.24 Top stream (ethanol) - 0.76 Bottom stream (water) - 0.89 Bottom stream (ethanol) - 0.11 Condenser heat (Kw) - -3700.96 Heat reboiler (Kw) - 3610.79 Machine processor (s) 40.960 46.919 Table 3: Initial conditions to simulate the extraction process. EMSO EMSO (Homotopy) Number of stages 35 35 Feed Stage (watel/ ethanol) 6 6 Feed Stage ethylene glycol 4 4 RR 0.5 0.5 Bottom stream 125 kmol/h 125 kmol/h Pressure (P) 1 atm 1 atm Ef. Murphee 1 1 Feed stream ethylene glycol 100 kmol/h 100 kmol/h Feed stream (F) 300 kmol/h 300 kmol/h Feed pressure (Pf) 1 atm 1 atm Feed temperature (Tf) 350 K 350 K Composition (F) / water 0.15 0.15 Composition (F) / ethanol 0.85 0.85 1067 A model for an extractive distillation column was built to separate ethanol-water using ethylene glycol as solvent. Simulations were made in AspenPlus® V.7.3 simulator with the same initial conditions as presented in Table 1 to validate the model. After the Homotopy Continuation method was applied, simulations for extractive distillation were performed with and without the method using the initial conditions in Table 3. As showed in Table 4, both simulations converged, and identical results were obtained. The difference between them is the time of simulation. Indeed, the models used found excellent results. The column top stream was able to achieve 0.925 (molar percentage), reaching the purity required by the market to sell ethanol as anhydrous ethanol. Besides that, the top stream did not presented significant amount of ethylene glycol. Moreover, the bottom stream left aside a small amount of ethanol, approximately 0.005 (molar percent), which is an acceptable loss for this process. Table 4: Results of extractive simulation. EMSO EMSO (Homotopy) Top stream (water) 0.075 0.075 Top stream (ethanol) 0.925 0.925 Top stream (ethylene glycol) 0 0 Bottom stream (water) 0.195 0.195 Bottom stream (ethanol) 0.005 0.005 Bottom stream (ethylene glycol) 0.8 0.8 Condenser heat (Kw)) -4500.42 -4500.42 Heat reboiler (Kw) 5085.22 5085.22 Machine processor (s) 31.176 81.736 4. Conclusions The robustness of the model of conventional distillation and extractive distillation was improved by Homotopy Continuation method because the convergence was accomplished even with poor initial conditions. However, models with this method requires a higher computational effort to converge. The algorithm used was efficient, considering the results obtained with the computational implementation are similar to those available in the literature for both, conventional and extractive distillation. References ANP (2011) - National Agency of Petroleum, natural gas and biofuels. Available at: . Accessed on: August 19, 2014. Ensinas, Av, (2008). Thermal integration and optimization thermoeconomic applied to industrial processes of sugar and ethanol from sugarcane. Doctoral Thesis, Unicamp, Campinas-SP. Goldemberg, J., Nigro, F, Rabbit, S. (2008). Bioenergy in São Paulo: Current State, Prospects, Barriers and Proposals. Official Press of the State of São Paulo, 152 p. Haelssig, J.B. et al (2011). Membrane Dephlegmation. A hybrid membrane separation process for efficient ethanol recovery. Journal of Membrane Science 381, p. 226–236. IEA (2011) – International Energy Agency. Technology Roadmap – Biofuels for Transport; REN21 (2012) – Renewable Energy Policy Network for the 21st century. Renewables 2012 – Global Status Report. Soares, R. P. (2007), “Manual EMSO”, avaliable at: http://www.enq.ufrgs.br/alsoc (accessed 28 October 2014). Tanskanen, J.; Pohjola, V. J (2000). A robust method for prediction state in a reactive distillation. Computers and Chemical Engineering, v. 24, n. 1, p. 81-88. WU, T.-M (2005). A study of convergence on the Newton-homotopy continuation method. Applied Mathematics and Computation, v. 168, n. 2, p.1169-1174. Zanette, A.F. (2010). Study of the transesterification of Jatropha oil using heterogeneous catalysts. Masters dissertation. Toledo, Paraná, Brazil. 1068