Introduction
Azeotropic or close boiling-point mixtures widely exist in the chemical industries such as ethanol/water, acetone/water and ethane/CO2. Due to the same or close volatilities of the components in the mixtures, they are extremely difficult to separate using the conventional distillation column. The extractive distillation (ED) is one of the widely used technologies for separation of azeotropic or close boiling-point mixtures due to less energy consumption and having more flexibility in the selection of solvents. An entrainer (or solvent) is introduced to change the relative volatility of the components through interacting differently with the components of the mixture [1]. This entrainer (or solvent) is generally non-volatile, has a high boiling point and is miscible with the mixture, but does not form an azeotropic mixture with the components. The original component with less interaction with the entrainer separates out as the top product of the extractive distillation column (EDC). The bottom product of the EDC consists of a mixture of the solvent and the other component, which can again be separated easily with a conventional distillation column (i.e., solvent recovery column, SRC) because the entrainer (or solvent) does not form an azeotrope with another component and the entrainer is usually much heavier. A typical conventional ED process for separation of a binary azeotropic mixture (A and B) is illustrated in Fig. 1. The entrainer and raw material are fed into the EDC with the component A drawn as the top distillate. The component B and entrainer leave from the EDC bottom and then are fed to the SRC. In the SRC, the pure component B is obtained on the top and the solvent is recovered at the bottom for recycle.
The conventional ED process is usually energy-intensive coupled with high capital cost. The dividing wall column (DWC), one of the most promising PI technologies, can reduce energy cost by 35% and decrease capital cost by 25% simultaneously [2]. Therefore, it has been used for the conventional ED process, leading to the development of the extractive dividing-wall column (EDWC). The applications of the EDWC for separation of azeotropic or close boiling-point mixtures include but not limited to ethanol/water separation in the bioethanol purification process [3] and CO2/ethane separation in the natural gas processing [4]. A significant reduction in energy and capital cost by 51.6% and 41.2% respectively is achieved compared to the conventional ED [4]. More importantly, a commercial EDWC has been established for separation of 2,6-lutidine/4-picoline [5].
Optimal design of the EDWC has received considerable attention during the last decade [6–8]. Several approaches have been proposed for the optimal design, including sequential sensitivity analysis-based optimisation methods, sequential quadratic programming (SQP) algorithms and genetic algorithms (GA). The sequential sensitivity analysis-based method conducts a series of process simulations using rigorous equilibrium stage models through changing one or a few design variables at each time while the other design variables are fixed [3,4,6,9]. This method can fail to converge in process simulation of the EDWC due to the entrainer recycle, complex configuration of the EDWC and strong nonideality of the mixture system. Furthermore, this method cannot identify an optimal design in most cases because the influence of the correlated design variables is not investigated. It often requires a large number of simulations to generate a best design. Although some heuristics are introduced to improve computational efficiency [6], they can identify a feasible design solution without guaranty of the solution optimality [10].
The SQP algorithm [11] built in a commercial process simulator such as Aspen Plus [12] is used to generate an optimal design of the EDWC [7,13,14]. The built-in SQP algorithm often fails to determine an optimal solution if a good initial point is not provided [15]. However, identification of a good initial point is not trivial. In addition, valid lower and upper bounds of the design variables are also required before optimisation, which are determined through a time-consuming sensitivity analysis [14]. The GA generates values of the independent variables and then passes them to a commercial process simulator such as Aspen Plus [12]. After each simulation, the GA evaluates the fitness function. Although this method is able to execute automatically, it requires a large number of fitness function evaluations (corresponding to process simulations), leading to large computational effort required [15]. This method can provide a good design solution without solution optimality guaranty.
All the aforementioned methods for optimal design of the EDWC conduct optimisation in the sequential modular simulation environment, which is inefficient due to either inaccurate derivative evaluation or not using derivative information. It is demonstrated that the gradient-based optimisation algorithm implemented in an equation-oriented (EO) environment is much more efficient for process optimisation compared to those in the sequential modular environment due to accurate derivative information used [16]. Although the gradient-based optimisation algorithms have been used for optimal design of thermally coupled ternary ED in the EO environment [17,18], they fail to find a feasible design of the EDWC for separation of azeotropic mixtures acetone/chloroform and ethanol/water. Waltermann et al. [19] generated an optimal design of the EDWC process using rigorous models in the EO optimisation environment General Algebraic Modeling System (GAMS). Although their solution method is shown to be efficient, complicated initialisation procedures and external functions are required to guarantee convergence and efficiency, which may not be suitable for optimal design of a complex EDWC such as the heat pump assisted EDWC.
To the best of our knowledge, it is still an open field to develop a systematic and efficient design method for optimal design of the EDWC using rigorous models due to high complexity of the EDWC configuration and strongly non-convex nonlinear physical properties of the azeotropic or close boiling-point mixture. Recently, Ma et al. [20] proposed a hybrid steady-state and time-relaxation feasible path optimisation algorithm where the original nonlinear programming (NLP) problem is decomposed into two levels. In the outer level, the SQP algorithm is used to navigate a group of decision variables (usually less than 1000), while in the inner level, process simulation based on the values of decision variables determined in the outer level is performed to generate values of other variables in the process. The pseudo-transient continuation (PTC) modelling approach is employed to resolve the difficulty in finding a feasible solution during process simulation in the EO environment. The PTC simulation using the tolerances-relaxation integration method [18] is performed only if the steady-state simulation fails to improve the computational efficiency. The computational results demonstrate the advantages of the hybrid optimisation algorithm in solving complex design optimisation problems. In this paper, we extend the hybrid steady-state and time-relaxation feasible path optimisation algorithm for optimal design of the EDWC. The tolerances-relaxation integration method [18] is further refined to ensure that the steady state is reached even at the initial integration phase with a large tolerance. The convergence of the hybrid algorithm is improved through allowing a larger tolerance for simulation in sensitivity analysis mode (algorithmic differentiation) when it diverges under small tolerance. In addition, valid lower bounds on the vapour flow rates at some stages and on the purity of the recovered entrainer are imposed to improve the computational efficiency when applying it to the optimisation of EDWC. Three examples from the literature [21–23] are used to illustrate the capability of the extended hybrid feasible path optimisation algorithm. The computational results demonstrate that the extended hybrid algorithm has good convergence performance and can find an optimal design of the EDWC within one hour in most cases. All the examples obtain much more economic designs than the existing literature reports. In addition, optimal design of the heat pump assisted EDWC is achieved for the first time.
Problem statement
Figure 2 illustrates a typical EDWC, which is used to separate an azeotropic or close boiling-point mixture A and B using an entrainer S. A dividing wall lies at the top of the column, leading to two condensers on the top of the column. The A and B mixture and the entrainer S are fed into the column from the left side of the dividing wall. The product A which has less interaction with the entrainer than B is drawn from the top left side of the column. The other product B is obtained from the top right side of the column. The entrainer S is obtained at the bottom of the column, which is recycled back to the entrainer feed. With this, the entire design problem can be stated as follows.
Given: an azeotropic or closing boiling-point mixture to be separated, the identified entrainer and column configuration, feed conditions and specified separation requirements.
Determine: 1) the entrainer make-up flow rate (FE) and its feed stage (NFE); 2) the feed stage of the raw material (NFR); 3) total number of stages (NM) in the main column and total number of stages (NS) in the side column; 4) reflux ratios of the main and side columns (RRM, RRS); 5) split fraction of the vapour stream to side column (SF) and the column bottom flow rate (FB); 6) connection stage between the main column and side column (NC).
Assumptions: 1) Vapor-liquid phase equilibrium is reached at each stage in the distillation column, i.e., the equilibrium stage model is used, which is usually accurate enough to simulate distillation column; 2) The column diameter is determined by the maximum diameter of the stages in the bottom section of the EDWC as the bottom section usually has larger vapour flow rate; 3) The heat transfer through the dividing wall is neglected. This assumption is usually given in literature. It should be interesting to investigate the influence of heat transfer through the dividing wall in the future; 4) The pressure drop on each stage is neglected. The operating pressure of the column is at atmospheric pressure. The pressure drop is usually small and does not have a large influence on the optimal design.
The objective is to minimize total annualised cost (TAC) which is calculated as total annualised capital cost plus the operating cost. The operating cost consists of energy and entrainer costs.
Modelling of EDWC
The column section model proposed by Pattison & Baldea [24] for modelling the normal DWC is extended to modelling the EDWC as shown in Fig. 3. There are five column sections in total. While the left side of the dividing wall in the EDWC is composed of sections 1–3, the right side of the dividing wall corresponds to the section 4. The bottom section in the EDWC is section 5. The make-up entrainer is mixed with the recycle entrainer before fed into the column. The raw material mixture to be separated is fed into the column through a different feed location, which is below that of the entrainer. This column section model can be extended for modelling other EDWC configurations. Using this column section model, the feed stage locations, connection stage location and the total number of stages in main and side columns can be determined easily from the number of stages in each column section. Some efforts [22,25] suggested that the number of stages on the left- and right-hand sides of the dividing wall should be identical to ensure the same pressure at the connection stages on both sides of the wall when trayed columns are used. However, such column may not be the most efficient. In the standard EDWC as shown in Fig. 2, the entrainer is often much heavier than other components, which is much easier to separate from the system. As a result, a much smaller number of stages is required on the right-hand side of the dividing wall in the EDWC, leading to some discrepancy in the number of stages on both sides. In this work, the number of stages on both sides of the dividing wall is allowed to be different in order to obtain better economic profit [8]. The equal pressure on the connection stages can be achieved easily through adjustment of the condenser pressure in the side column.
Each column section can be mathematically modelled using the bypass efficiency method [26], which is very similar to that of a conventional distillation column. In the bypass-efficiency method, the existence of a stage in a column section is represented by a bypass efficiency εj, which is treated as a continuous variable between 0 and 1 rather than a binary variable. This avoids the introduction of binary variables. As a result, only an NLP optimisation problem needs to be solved. The actual total number of stages (denoted as Nt) in a column section is given below,
Although the bypass efficiencies are treated as continuous variables between 0 and 1, they tend to be 0 or 1 at the globally optimal point because partial bypass (i.e., the values between 0 and 1) is thermodynamically inefficient [26]. However, some bypass efficiencies may be fractional in our solution due to locally optimal solutions identified using our hybrid algorithm. In addition, numerical errors may also cause the bypass efficiencies to be fractional. These fractional bypass efficiencies can be rounded to be 0 or 1 and a new optimal solution can be generated. It is demonstrated that the optimal solution does not change noticeably after the fractional bypass efficiencies round off in case studies of this work and the existing efforts [17,18].
The optimisation problem for optimal design of the EDWC (denoted as P0) can be represented in a general form as follows:
where x are process model variables to be optimised, f is the objective function representing TAC, h represents all equations (equality constraints) including mass balance, phase equilibrium, summation and heat balance (MESH) equations, process stream connection equations and some equations related with design requirements, and g refers to all inequality constraints including design specifications such as product purity and recovery requirements. The MESH equations with the bypass efficiency are presented in Appendix A (cf. Electronic Supplementary Material, ESM). The economic evaluation model is provided in Appendix B (cf. ESM). It should be noted that solving problem (P0) directly is extremely challenging due to the strongly non-convex nonlinear and coupled MESH equations, and the rigorous enthalpy, density and phase equilibrium models used.
Solution approach
The feasible path optimisation algorithm has demonstrated its capabilities for solving large-scale strongly non-convex nonlinear optimisation problem [27]. It decomposes the original problem P0 into two sub-problems including a small-scale optimisation problem (denoted as P1) in the outer level and a process simulation problem in the inner level [28]. The problem P1 is generated through elimination of all dependent variables using the corresponding equations in problem P0.
where xI represents independent variables (or decision variables) to be optimised, which are nI-dimensional, y and z are the process variables related with design requirements such as product purities and temperature difference in the heat exchanger. Note that y and z are dependent variables and determined by xI through an implicit function that is defined by the process simulation in the inner level. Similarly, the complex process model equations of the EDWC are also included in the simulation problem. As a result, the outer problem P1 is a small-scale optimisation problem. More details can be referred to Ma et al. [20]. The main challenge is how to solve the process simulation problem efficiently and reliably.
The feasible path optimisation algorithm often fails to provide a feasible solution when it is applied for optimal design of the complex EDWC using rigorous MESH distillation models. The main possible reason is due to the non-convergence of the steady-state process simulation in the EO environment [24]. The hybrid steady-state and time-relaxation feasible path optimisation algorithm from Ma et al. [20] is extended to solve the optimisation problem (P0). This hybrid algorithm is an improvement to the existing feasible path optimisation algorithm in which the combination of steady-state simulation and PTC simulation is used to improve the convergence of process simulation.
PTC modelling approach
The PTC modelling approach reformulates algebraic equations in the steady-state model into differential algebraic equations (DAE) with an equivalent steady-state solution of the original models [24]. The derived DAE model is the so-called PTC model. It is reported that process simulation using the PTC model can lead to better convergence than steady-state simulation due to the continuation on the pseudo time [24]. This is of great importance for process optimisation using complex rigorous unit operation models. The PTC distillation model from Ma et al. [17] is used for simulation of the EDWC, which is provided in Appendix C (cf. ESM).
Hybrid steady-state and time-relaxation optimisation algorithm
The hybrid steady-state and time-relaxation feasible path optimisation algorithm from Ma et al. [20] is depicted in Fig. 4. In this hybrid algorithm, the PTC simulation is performed only when the steady-state simulation fails as the PTC simulation is much more time-consuming than the steady-state simulation. The hybrid algorithm starts from an initial point x0 which must lead to a converged process simulation. At this converged point, the derivative information is evaluated and passed to a gradient-based optimiser such as the SQP optimiser. The optimiser checks whether the Karush-Kuhn-Tucker (KKT) conditions are satisfied at the current point. If it is satisfied, the optimisation is completed successfully with an optimal solution. Otherwise, new values of the decision variables xk+1 are returned by the optimiser, and a new steady-state simulation is performed. If the new simulation converges, the next iteration (k+1) starts. Otherwise, a PTC simulation is executed to generate a converged solution instead before the next iteration (k+1) starts. The iterations continue until the KKT conditions are satisfied, indicating the successful identification of an optimal solution, or the iterations may terminate prematurely without finding the optimal solution. There are two main reasons that may cause failure of an optimisation. One is the maximum number of line searches is reached without an acceptable step length found. The other one is that the directional derivative along the derived descent direction is positive. To guarantee the robustness of the algorithms, existing SQP optimisers may use some tricks to alleviate the above problems.
It should be noted that the tolerances-relaxation integration method proposed by Ma et al. [18] is also used in the hybrid algorithm to reduce computational effort from the PTC simulation. This integration method imposes a large integration tolerance in the beginning and adopts a tight tolerance later to guarantee the accuracy of the final steady-state solution. In addition, we found that the SQP optimiser used in the current work may fail due to the positive directional derivative obtained, which is usually caused by an indefinite or negative semi-definite approximate Hessian matrix. To improve the robustness of the hybrid algorithm, a new optimisation with a cold start is performed from the point where the last optimisation fails. The cold start forces the Hessian matrix to a default well-posed approximate Hessian matrix so that a descent direction leading to negative directional derivative can be obtained.
Improvements
Some improvements are made to improve computational efficiency of the hybrid algorithm in this work. First, the tolerances-relaxation integration method used in the PTC simulation from Ma et al. [18] is refined. In the previous efforts [18,20], the PTC simulation used a large integration tolerance such as 1 × 10–3 in the beginning and performed for a presumably long enough integration time T to ensure the solution close to the steady-state solution before the required tolerance was used. This is of great importance to reduce computational time of the PTC simulation, especially when a small tolerance such as 1 × 10–10 is used. However, it is often difficult to ensure that the integration time T provided is always long enough during the optimisation of the EDWC possibly due to large ratio (usually larger than 100) between the entrainer recycle flow rate and the fresh entrainer flow rate. To avoid such situation, the tolerances-relaxation integration method is improved through increasing the integration time T by 10 times until the steady-state simulation with a large tolerance is converged, as highlighted in Fig. 5. The improved tolerances-relaxation integration method is described below.
Step 1: Given initial values of the decision variables x, a required convergence tolerance tol_req, an integration period T and a list of tolerances tol in a descent order with the last one being the required tolerance tol_req. The outer and inner iteration indicators i and j are set as zero.
Step 2: The PTC simulation is integrated for 10j·T (pseudo time) with the initial value of j = 0. If the PTC simulation does not converge, the whole PTC simulation fails and then goes to step 8.
Step 3: The current converged point D0 is recorded and a steady-state simulation with toli is performed. If the steady-state simulation converges, then go to step 6. Otherwise go to the next step.
Step 4: If j>J_max, the entire PTC simulation fails and goes to step 8. Otherwise go to step 5.
Step 5: The initial point used for the simulation is reset to the point D0 and j is increased by 1, then return to step 2.
Step 6: If the current tolerance toli is equal to tol_req, then the simulation completes successfully and go to step 8. Otherwise, if toli>tol_req, the current steady-state solution S0 is recorded and a steady-state simulation with the tolerance of tol_req is performed, then go to the next step.
Step 7: If the steady-state simulation converges, the whole simulation completes and go to step 8. Otherwise, the initial point used for process simulation is reset to the point S0. The outer iteration indicator i is increased by 1 and return to step 2.
Step 8: Return.
Second, sensitivity analysis is performed after a successful process simulation to generate accurate derivative information for the gradient-based optimiser. The sensitivity analysis is also realized through process simulation but with derivative information calculated simultaneously together with process variables. We observe that the steady-state simulation in the sensitivity analysis mode can fail at some points where the steady-state simulation without sensitivity analysis converges. The possible reason is due to the contradiction between the machine epsilon of the double precision floating-point number adopted in the Aspen Custom Modeler (ACM) and the fairly small simulation tolerance used. For example, when the simulation tolerance is set as 1 × 10–10, the order of magnitude of the terms in the sensitivity equations should not exceed 6 in order to guarantee convergence under the machine epsilon of around 1.1 × 10–16 [29]. However, the order of magnitude of such terms may exceed 6 in some cases, causing unreliable convergence performance. The failure of simulation in the sensitivity analysis mode leads to no derivative information available for the gradient-based optimiser to generate a new search direction. As a result, the optimisation terminates prematurely. To avoid such case, a relaxed tolerance which is larger than the required tolerance is used for the steady-state simulation in the sensitivity analysis mode. Although the relaxed tolerance may reduce the accuracy of derivatives, it is used only in several optimisation steps and it is often still small enough such as 1 × 10–7, 1 × 10–8, and 1 × 10–9 based on our observation. It is worthwhile to generate derivative information to enable the optimisation to continue with sacrifice of some accuracy instead of failure in the optimisation. The improved procedure is illustrated in Fig. 6.
In Fig. 6, we first conduct the steady-state simulation using the required tolerance (tols = tol_req) in the sensitivity analysis mode. If the steady-state simulation converges, the derivative information evaluation completes successfully. Otherwise, if the current tolerance is smaller than the first tolerance tol0 in the list tol, the tolerance is magnified to 10 times and then go to the next iteration. If tols≥tol0, then the derivative information evaluation terminates with failure because no derivative information is available to generate new values of the decision variables.
To further improve the computational performance of the extended hybrid feasible path optimisation algorithm, valid lower bounds of some decision variables are imposed to remove infeasibility regions. First, a reasonable lower bound on the recycled entrainer purity (e.g., 0.90) is enforced to eliminate the optimisation region with highly impure recycled entrainer. However, this lower bound cannot be set too high (e.g., 0.95) because the optimal solution may be eliminated. In this work, the lower bound on the recycled entrainer purity is set to be 0.90 as the optimal purity of the recovered entrainer is usually higher than this value. Second, process simulation often becomes slow or even stuck when a single-phase region is reached at some stages in the distillation column since MESH equations are not applicable for a single-phase fluid at each stage. We also observe that the vapour flow rate may become zero in sections 2 and 3 during optimisation because the feed and the entrainer streams are usually at sub-cooled liquid pahse. To avoid such infeasible case, we impose a lower bound on the vapour flow rates (e.g., 1 or 10 kmol∙h–1) in sections 2 and 3 without loss of the optimal solution as the single-phase zone with no mass transfer rarely appears on a stage in an efficient distillation column.
Case studies
Three case studies are used to illustrate the capability of the extended hybrid feasible path optimisation algorithm for optimal design of the EDWC. The first case study separates acetone/chloroform mixture using dimethylsulfoxide (DMSO) as the entrainer. The second and the third case studies separate ethanol/water mixture with ethylene glycol (EG) as the entrainer. While the first two cases use the EDWC, the third case uses the heat pump assisted EDWC. The objective of all three case studies is to minimize TAC including total annualised capital cost and operating cost. Detailed correlations for TAC calculation are provided in Appendix B. In order to identify the globally optimal solution, six initial points are randomly generated and used to initialise the extended hybrid algorithm in all case studies. This allows to test the convergence performance of the extended algorithm. In all case studies, the simulation tolerance is set as 1 × 10–10 with optimisation tolerance of 1 × 10–5. The list of tolerances used in the tolerances-relaxation integration method is [1 × 10–3, 1 × 10–6, 1 × 10–10]. The initial integration time T is 2000 h. The PTC model of the EDWC is implemented in ACM [30]. The improved hybrid feasible path optimisation algorithm is implemented in Python [31] with the built-in solver slsqp [32] used as the NLP solver. All case studies are solved on a desktop with 3.20 GHz intel core i7 processor and 16 GB RAM running 64-bit windows operating system.
Case study 1: acetone/chloroform separation
Acetone and chloroform can form a maximum boiling point azeotrope. The azeotropic temperature is 67.51 °C and the azeotropic molar composition is 33.91% acetone and 66.09% chloroform. The feed rate is 100 kmol∙h–1 with molar compositions of acetone and chloroform being around 0.5 and 0.5, respectively. While the purity of acetone product has to reach at least 99.5% (mole), it should be at least 99.55% for chloroform. Luyben [21] designed a conventional ED process where DMSO was used as the entrainer. High-pressure steam and cooling water were used as hot and cold utilities, respectively. To reduce energy and capital costs, the EDWC is employed in this work. The entire optimisation problem involves 9399 equations, 16 inequality constraints, and 10570 variables including 90 decision variables. The initial values, and lower and upper bounds of decision variables are listed in Table 1.
Tab.1 Initial values and lower and upper bounds of the decision variables for case study 1 |
| Decision variable | Initial value | Lower bound | Upper bound |
|---|---|---|---|
| FE/(kmol∙h–1) | 0.01 | 0.0001 | 10 |
| RRM/(kmol∙kmol–1) | 1.0 | 0.1 | 10 |
| RRS/(kmol∙kmol–1) | 1.0 | 0.01 | 10 |
| SF/(kmol∙kmol–1) | 0.5 | 0.1 | 0.99 |
| FB/(kmol∙h–1) | 100 | 10 | 500 |
| N1 | – | 0 | 10 |
| N2 | – | 0 | 20 |
| N3 | – | 0 | 30 |
| N4 | – | 0 | 15 |
| N5 | – | 0 | 10 |
The computational performance of the extended hybrid algorithm starting from six initial points are provided in Table 2. The corresponding optimal solutions are provided in Table 3. As seen from Table 2, the computational time required is within 1 h in most cases except that starting from εj = 0.3 which requires 2 h. The smallest computational time required is about 20 min using the initial point εj = 0.1. From this initial point, the best locally optimal solution of 6.081 × 105 $∙year–1 can be obtained, as demonstrated in Table 3. A minor problem in this optimal solution is that the number of stages in the column section 4 is fractional due to several bypass efficiencies in this section being fractional values between 0 and 1. This is possibly because only the tray cost is involved in this column section and it is much smaller compared to the TAC. As a result, the TAC is much less sensitive to the values of bypass efficiencies in this section than those in other sections. In addition, it is not difficult to recover the entrainer in this column section 4 due to large relative volatility between chloroform and DMSO. After rounding off the fractional bypass efficiencies to be 0 or 1, the TAC slightly increases to 6.082 × 105 $∙year–1 with satisfied product purity requirements.
Tab.2 Computational performance of the extended hybrid algorithm for case study 1 from six different initial points |
| εj | 0.1 | 0.3 | 0.5 | 0.7 | 0.9 | 1.0 |
|---|---|---|---|---|---|---|
| Convergence | Yes | Yes | Yes | Yes | Yes | Yes |
| Iteration step | 239 | 302 | 409 | 509 | 380 | 262 |
| Simulation time/s | 1259 | 7722 | 2263 | 3698 | 1650 | 1714 |
| TAC/(105 $∙year–1) | 6.081 | 6.116 | 6.106 | 6.124 | 6.137 | 6.097 |
Tab.3 Optimal design of the EDWC for case study 1 from six different initial points |
| εj | 0.1 | 0.3 | 0.5 | 0.7 | 0.9 | 1.0 |
|---|---|---|---|---|---|---|
| FE/(kmol∙h–1) | 0.0008 | 0.0003 | 0.0007 | 0.0016 | 0.0006 | 0.0009 |
| RRM/(kmol∙kmol–1) | 0.70 | 0.74 | 0.74 | 0.51 | 0.81 | 0.68 |
| RRS/(kmol∙kmol–1) | 0.12 | 0.11 | 0.11 | 0.12 | 0.11 | 0.11 |
| SF/(kmol∙kmol–1) | 0.31 | 0.31 | 0.31 | 0.31 | 0.30 | 0.31 |
| FB/(kmol∙h–1) | 115.37 | 110.69 | 110.33 | 136.69 | 106.89 | 114.11 |
| N1 | 2.0 | 2.5 | 2.0 | 2.1 | 2.0 | 2.0 |
| N2 | 10.0 | 10.0 | 10.0 | 11.0 | 9.1 | 11.0 |
| N3 | 22.0 | 23.9 | 24.0 | 17.0 | 25.0 | 23.0 |
| N4 | 5.9 | 11.0 | 7.8 | 5.4 | 8.5 | 8.4 |
| N5 | 4.0 | 4.0 | 4.0 | 5.0 | 4.0 | 4.0 |
| TAC/(105 $∙year–1) | 6.081 | 6.116 | 6.106 | 6.124 | 6.137 | 6.097 |
The best optimal design with TAC of 6.081 × 105 $∙year–1 is illustrated in Fig. 7. As can be seen from Fig. 7, the dividing wall in the optimal design does not reach the column top due to a smaller number of stages on the right side of the dividing wall required, which is different from the one illustrated in Fig. 2. However, these two configurations were thermodynamically equivalent.
For a fair comparison, the optimal design of the conventional ED is also generated using the extended hybrid algorithm, which is illustrated in Fig. 8. The economic comparison among the optimal design of the conventional ED, the design of the conventional ED from Luyben [21], and the optimal design of the EDWC is provided in Table 4. It should be noted that the costs of the conventional ED in Table 4 are different from those reported in Luyben [21] due to the use of different formula for column diameter calculation.
Tab.4 Economic comparison between EDWC and conventional ED for case study 1 |
| Item | Conventional ED (Luyben [28]) | Conventional ED (optimal) | EDWC |
|---|---|---|---|
| Annualised capital cost/(105 $∙year–1) | 2.226 | 2.326 | 2.320 |
| Energy cost/(105 $∙year–1) | 5.025 | 3.951 | 3.754 |
| Entrainer cost/(105 $∙year–1) | 0.039 | 0.012 | 0.008 |
| TAC/(105 $∙year–1) | 7.290 | 6.289 | 6.082 |
As can be seen from Table 4, the TAC in the optimal design of the conventional ED is 6.289 × 105 $∙year–1, which is reduced by 13.7% compared to the TAC (7.290 × 105 $∙year–1) in the heuristic-based design from Luyben [21] due to the significant reduction in energy cost, which is about 21.4%. The capital cost is a bit higher (4.5%) than that of Luyben [28] due to the total number of stages increased from 34 to 49 although the diameters of both columns in our design decreased evidently (1.37 vs. 1.50 m and 1.02 vs. 1.14 m). It is also observed that the purity of the recycle entrainer in the optimal design of the conventional ED is 0.9986 which is different from that of 0.9999 determined using heuristics in the work of Luyben [28]. This result is consistent with the result from Ma et al. [17] for separation of the benzene/cyclohexane/toluene mixture, where the recycle entrainer purity was only 0.85.
On the other hand, the TAC of the optimal EDWC is 6.082 × 105 $∙year–1, which is reduced by only 3.3% compared to that of the optimal conventional ED. This is because the diameter of EDWC is much larger than that of the ED although one column is required in the EDWC, leading to a slight reduction in the annualised capital cost. The reboiler duty in the EDWC is 1970 kW which is reduced by 4.5% compared to that of 2062 kW in the optimal ED. Such small reduction in TAC indicates that the EDWC may not be attractive for acetone/chloroform separation as it is much more difficult to construct and control compared to the ED.
Case study 2: bioethanol purification
Bioethanol is an attractive alternative to traditional fossil fuels, and hence has received heated attention from both industries and academics [3]. Usually, it is produced from biomass fermentation and ethanol in the fermented broth is dilute (e.g., 5–12 wt-%) [34]. As ethanol and water form an azeotrope, the diluted bioethanol is separated through a three-column process including the pre-concentration column followed by a double-column ED process (EDC+ SRC). As this process is energy-intensive, many researchers have attempted to reduce its high energy consumption using the EDWC [3,7,13,23,35]. While most work generated the EDWC through a simple combination of the EDC column and the SRC column, Kiss and Ignat [22] proposed an innovative EDWC that combined the three columns together as illustrated in Fig. 9. The feed rate is 6516.03 kmol∙h–1 with the molar fraction of ethanol being 0.042. The target product purities are 99.9 wt-% ethanol and 99.8 wt-% water. The required recoveries for ethanol and water are 99.81% and 99.99% respectively. Medium-pressure (MP) steam is used in the reboiler on the left-hand side, while high-pressure (HP) steam is used in the reboiler on the right-hand side (see Fig. 9). The operating pressure is 1.013 × 105 Pa.
We use the extended hybrid algorithm to generate the optimal configuration of the EDWC proposed by Kiss and Ignat [22] with the same feed conditions and design requirements. The optimisation problem involves 14038 equations, 20 inequality constraints, 15805 variables including 125 design variables. The initial values, and lower and upper bounds of the decision variables are given in Table 5. Here, VF is the reboiler vaporisation fraction of the main column.
Tab.5 Initial values and lower and upper bounds of the decision variables for case study 2 |
| Decision variable | Initial value | Lower bound | Upper bound |
|---|---|---|---|
| FE/(kmol∙h–1) | 0.01 | 0.0001 | 10 |
| RRM/(kmol∙kmol–1) | 1.0 | 0.1 | 100 |
| VF/(kmol∙kmol–1) | 0.5 | 0.1 | 0.5 |
| SF/(kmol∙kmol–1) | 0.5 | 0.1 | 0.99 |
| FB/(kmol∙h–1) | 300 | 10 | 3000 |
| N1 | – | 0 | 10 |
| N2 | – | 0 | 60 |
| N3 | – | 0 | 20 |
| N4 | – | 0 | 20 |
| N5 | – | 0 | 10 |
The computational performance of the extended algorithm from six different initial points is provided in Table 6. The corresponding optimal solutions are provided in Table 7.
Tab.6 Computational performance of the extended hybrid algorithm for case study 2 from six different initial points |
| εj | 0.1 | 0.3 | 0.5 | 0.7 | 0.9 | 1.0 |
|---|---|---|---|---|---|---|
| Converged | Yes | Yes | Yes | Yes | No | Yes |
| Iteration | 760 | 908 | 591 | 577 | 398 | 601 |
| Simulation time/s | 2121 | 3742 | 2309 | 2204 | 1501 | 3327 |
| TAC/(M$∙year–1) | 5.382 | 5.383 | 5.388 | 5.405 | 5.388 | 5.382 |
As seen from Table 6, the extended algorithm fails to converge from the initial point εj = 0.9, whilst it successfully identifies locally optimal solutions from the other five initial points. Note that the optimisation starting from the initial point εj = 0.9 terminates at a point close to the optimal solution as its TAC is quite close to the optimal value and the product purities and recoveries are also close to the required. The number of iterations varies from 591 to 908 with the optimisation time ranging from 2121 to 3742 s. The smallest TAC of 5.382 M$∙year–1 is found from the initial points εj = 0.1 and εj = 1.0. The best optimal design is illustrated in Fig. 9. From Fig. 9, the total number of stages is 76 with the reboiler duty of 19536 kW. Note that we do consider the entrainer cost in the objective function to reduce the loss of the entrainer, which was not considered in the design of Kiss and Ignat [22]. As the entrainer EG which is mainly lost in water product as impurity is quite expensive, the optimal design tries to reduce its consumption. This is the main reason why the mass purity of water product in our design is 0.9987 that is higher than the required one (0.9980). In other words, the purity requirement constraint is inactive. This indicates that it may exclude optimal deign to enforce product purity or recovery requirements as equalities, which although may alleviate solution difficulty.
Tab.7 Optimal solutions for case study 2 from six different initial points |
| εj | 0.1 | 0.3 | 0.5 | 0.7 | 0.9 | 1.0 |
|---|---|---|---|---|---|---|
| FE/(kmol∙h–1) | 1.90 | 1.87 | 1.93 | 1.97 | 1.85 | 1.87 |
| RRM/(kmol∙kmol–1) | 1.89 | 1.87 | 1.90 | 1.93 | 1.87 | 1.88 |
| VF/(kmol∙kmol–1) | 0.16 | 0.16 | 0.16 | 0.16 | 0.16 | 0.16 |
| SF/(kmol∙kmol–1) | 0.25 | 0.24 | 0.26 | 0.27 | 0.24 | 0.25 |
| FB/(kmol∙h–1) | 112.95 | 110.65 | 115.89 | 119.67 | 109.16 | 111.07 |
| N1 | 2.0 | 2.0 | 2.0 | 2.0 | 2.0 | 2.0 |
| N2 | 54.0 | 59.0 | 49.0 | 44.0 | 63.0 | 58.0 |
| N3 | 12.0 | 12.0 | 12.0 | 12.2 | 12.0 | 12.0 |
| N4 | 12.0 | 12.0 | 12.0 | 12.2 | 12.0 | 12,0 |
| N5 | 6.0 | 6.0 | 6.0 | 6.0 | 6.0 | 6.0 |
| TAC/(M$∙year–1) | 5.382 | 5.383 | 5.388 | 5.405 | 5.388 | 5.382 |
For a fair comparison with the design of Kiss and Ignat [22], we use the extended hybrid algorithm to solve the optimisation problem again but without considering the entrainer cost in the objective function. The optimal TAC is about 3.969 M$∙year–1. The optimal design is illustrated in Fig. 10. From Fig. 10, our optimal design requires 49 column stages, which is similar to that (42) of Kiss and Ignat [22]. The total reboiler duty is 19480 kW, which is reduced by 24.4% than that (25775 kW) from Kiss and Ignat [22] due to a much smaller reflux ratio (1.7 vs. 3.4) required. The reduction in the reflex ratio indicates that a smaller reflux rate (481 vs. 923 kmol∙h–1) and a smaller vapour flow rate (753 vs. 1194 kmol∙h–1) are required in the main column. As a result, the heat duty of the main reboiler is 6478 kW in our optimal design, which is reduced by 48.4% compared to that in the design of Kiss and Ignat [22] (12563 kW). The reboiler duty in the side column is similar in both designs (13002 vs. 13028 kW). Furthermore, the entrainer feed rate is reduced from 335 to 153 kmol∙h–1 in our optimal design. While the number of stages increases from 14 to 29 in the section 2 to enhance the separation effectiveness of the ED section, it decreases from 17 to 12 in both sections 3 and 4. To make up the decrease in the number of stages in the section 3 (i.e., the pre-concentration section), more vapour from the section 5 is fed into this section 3 with a larger split fraction (0.78 in our design vs. 0.4 in the design of Kiss & Ignat [22]).
Finally, we compare the optimal design with and without consideration of the entrainer cost, which are denoted as the optimal designs A and B respectively for convenience. The optimal design B requires much smaller number of stages than that in the optimal design A (49 vs. 76) due to much more entrainer (153 vs. 115 kmol∙h–1) fed to the column in the optimal design B, making the ED easier. Also, the design A requires a little more heat duty than the design B (19536 vs. 19480 kW). Therefore, when only capital cost and energy cost are considered, the TAC of optimal design A is higher than that of optimal design B (4.16 vs. 3.97 M$∙year–1). However, 3.25 kmol∙h–1 of the entrainer make-up is required in the optimal design B, which is 71.1% more than that (1.90 kmol∙h–1) required in the optimal design A. As a result, the TAC increases to 6.07 M$∙year–1 when the entrainer make-up cost is also added in the optimal design B, which is increased by 12.8% compared to the TAC of 5.38 M$∙year–1 in the optimal design A. The entrainer make-up cost accounts for 22.8% of the TAC in the optimal design A whilst it accounts for 34.6% of the TAC in the optimal design B. This indicates the significance of taking the entrainer cost into account in the optimal design of the EDWC.
Case study 3: heat pump assisted EDWC
Luo et al. [23] proposed a heat pump assisted EDWC for purification of bioethanol, as shown in Fig. 11. The vapour stream from the top is compressed and then used to vaporise the bottom stream from the side column in a heat exchanger HX1 in order to decrease the reboiler duty of the side column. The outlet hot stream is fully condensed. Some of them return to the column as a reflux stream, whilst the remaining are used as the top product. The entrainer recovered from the column bottom is used to preheat the feed in another heat exchanger HX2.
Besides all decision variables discussed in the EDWC, three additional decision variables need to be determined including outlet pressure of the compressor (P) and areas of the two heat exchangers HX1 and HX2 (A1 and A2). In order to avoid getting stuck in some worse locally optimal solution, the outlet hot stream of the heat exchanger HX2 should be in a two-phase zone or liquid-phase zone through imposing vapour fraction to be less than 1. This allows the latent heat released from the condensation to be used for heating the stream from the column bottom. As a result, the reboiler duty required is reduced significantly. At the same time, the cooling duty in the condenser that is used to cool this outlet hot stream is also reduced.
The initial values and lower and upper bounds of the decision variables are provided in Table 8. The isentropic compression model is used for the compressor with the isentropic efficiency of 0.72 and mechanical efficiency of 1.0.
Tab.8 Initial values and lower and upper bounds of decision variables for case study 3 |
| Decision variable | Initial | Lower | Upper |
|---|---|---|---|
| FE/(kmol∙h–1) | 0.01 | 0.0001 | 10 |
| RRM/(kmol∙kmol–1) | 1.0 | 0.1 | 100 |
| VF/(kmol∙kmol–1) | 0.5 | 0.1 | 0.5 |
| SF/(kmol∙kmol–1) | 0.5 | 0.1 | 0.99 |
| FB/(kmol∙h–1) | 300 | 10 | 3000 |
| N1 | – | 0 | 10 |
| N2 | – | 0 | 60 |
| N3 | – | 0 | 20 |
| N4 | – | 0 | 20 |
| N5 | – | 0 | 10 |
| P/(1.013 × 105 Pa) | 3 | 1.001 | 10 |
| A1/m2 | 100 | 0 | 1 × 106 |
| A2/m2 | 10 | 0 | 1 × 106 |
The computational performance of the improved hybrid algorithm using six initial points is provided in Table 9. The optimal decision variables are provided in Table 10.
Tab.9 Computational performance of the improved hybrid algorithm for case study 3 starting from six initial points |
| εj | 0.1 | 0.3 | 0.5 | 0.7 | 0.9 | 1.0 |
|---|---|---|---|---|---|---|
| Converged | Yes | Yes | Yes | No | Yes | Yes |
| Iteration step | 1519 | 2175 | 2234 | 332 | 2392 | 1235 |
| Simulation time/s | 6566 | 6146 | 7238 | 1146 | 7624 | 3916 |
| TAC/(M$∙year–1) | 4.879 | 4.829 | 4.851 | 5.715 | 4.834 | 4.847 |
Tab.10 Optimal solutions for case study 3 from six initial points |
| εj | 0.1 | 0.3 | 0.5 | 0.7 | 0.9 | 1.0 |
|---|---|---|---|---|---|---|
| FE/(kmol∙h–1) | 1.94 | 1.71 | 1.83 | 1.99 | 1.79 | 1.87 |
| RRM/(kmol∙kmol–1) | 2.02 | 2.18 | 2.16 | 3.56 | 2.23 | 2.02 |
| VF/(kmol∙kmol–1) | 0.15 | 0.15 | 0.15 | 0.15 | 0.15 | 0.15 |
| SF/(kmol∙kmol–1) | 0.27 | 0.32 | 0.31 | 0.64 | 0.33 | 0.27 |
| FB/(kmol∙h–1) | 118.65 | 113.50 | 120.35 | 294.78 | 120.87 | 114.61 |
| N1 | 2.0 | 2.0 | 2.0 | 6.4 | 2.0 | 2.0 |
| N2 | 46.0 | 56.0 | 45.0 | 25.6 | 45.0 | 52.0 |
| N3 | 12.0 | 12.1 | 12.0 | 10.3 | 12.0 | 12.0 |
| N4 | 12.0 | 12.1 | 12.0 | 9.9 | 12.0 | 12.0 |
| N5 | 6.0 | 6.0 | 6.0 | 6.7 | 6.0 | 6.0 |
| P/(1.013 × 105 Pa) | 3.26 | 3.08 | 3.08 | 2.20 | 3.08 | 3.15 |
| A1/m2 | 854.10 | 1003.50 | 986.63 | 3904.22 | 1018.83 | 914.72 |
| A2/m2 | 37.04 | 35.48 | 37.68 | 5.04 | 37.84 | 35.75 |
| TAC/(M$∙year–1) | 4.879 | 4.829 | 4.851 | 5.715 | 4.834 | 4.847 |
As can be seen from Table 9, the computational time for optimal design of the heat pump assisted EDWC is almost double or triple the time required for optimal design of a normal EDWC in case studies 1–2. For instance, the computational time for optimal design of the EDWC staring from εj = 0.1 is 2121 s (see Table 6), which is increased to 6566 s for optimal design of the heat pump assisted EDWC. This is because the heat integration configuration makes the process more complex. The extended hybrid algorithm fails to generate an optimal solution using the initial point εj = 0.7. The best TAC of 4.829 M$∙year–1 is identified starting from the initial point εj = 0.3. As seen From Table 10, the numbers of stages in sections 3 and 4 are fractional, not integer due to fractional values of some bypass efficiencies obtained. Specifically, there is only one fractional bypass efficiency obtained in the column section 3, whilst five bypass efficiencies are fractional in the column section 4. After the fractional bypass efficiencies are rounded to 0 or 1, the TAC becomes 4.828 M$∙year–1, which is almost the same as the original TAC of 4.829 M$∙year–1. The slight decrease of TAC is mainly because the outlet temperature difference of the HX2 changes from 10.00000 °C to 9.99991 °C, which does not exactly match the minimum temperature approach (10 °C). However, such violation is negligible and should be acceptable. The optimal design with TAC of 4.828 M$∙year–1 is illustrated in Fig. 11. Similar to the case study 2, the mass purity of water is 0.9988 which is also higher than the required specification (0.9980).
We compare the optimal design of the heat pump assisted EDWC with that of the normal EDWC. The comparative results are shown in Table 11. As can be seen from Table 11, the capital cost of the optimal heat pump-assisted EDWC is 4.064 M$, which is and two heat exchangers used in the heat-pump assisted EDWC. However, utility cost is significantly reduced from 3.316 to 1.757 M$∙year–1, which is reduced by 47.0% benefited from heat integration. Although electricity cost is increased from 0 to 0.614 M$∙year–1, the overall operating cost still decreases from 4.541 to 3.473 M$∙year–1, which is reduced by 23.5%. Finally, the overall TAC is reduced by 10.3% as the reduction in operating cost outweighs the increase in annualised capital cost.
Tab.11 Cost breakdown in the optimal design of the normal EDWC in case study 2 and the optimal design of the heat pump assisted EDWC in case study 3 |
| Cost | Normal EDWC | Heat pump-assisted EDWC |
|---|---|---|
| Capital cost/M$ | ||
| Column cost | 1.707 | 1.746 |
| Heat exchanger cost | 0.818 | 1.185 |
| Compressor cost | 0 | 1.134 |
| Operation cost/(M$∙year–1) | ||
| Utility cost | 3.316 | 1.757 |
| Electricity cost | 0 | 0.614 |
| Entrainer cost | 1.225 | 1.102 |
| Total cost | ||
| Total capital cost/M$ | 2.524 | 4.064 |
| Total operation cost/(M$∙year–1) | 4.541 | 3.473 |
| TAC/(M$∙year–1) | 5.382 | 4.828 |
For a fair comparison with the design of Luo et al. [23], we use the same specifications as those in Luo et al. [23] to generate another optimal design. The optimal design is shown in Fig. 12. As seen from Fig. 12, the reflux ratio is only 1.84 in our optimal design, whilst it is 4.20 in the design of Luo et al. [23]. A higher reflux ratio is used in Luo et al. [23] because they want to ensure that the heat released from the condenser is enough for the side reboiler. Intuitively, this seems to be reasonable. However, the large reflux ratio also increases the heat duty of the reboiler in the main column and the burden of the compressor. Note that the side reboiler only requires MP or even low pressure steam as hot utility, while the main reboiler needs more expensive HP steam. Electricity cost required by the compressor is even more expensive than HP steam. Therefore, in our optimal design, a small reflux ratio is used and the condenser only provides a part of heat required by the side reboiler in order to minimize TAC. The specific energy consumption in our optimal design and the design of Luo et al. [23] are listed in Table 12.
As can be seen from Table 12, the design of Luo et al. [23] requires much more HP steam and electricity compared to our optimal design, although their design does not require any MP steam, whilst our design consumes 3999 kW of MP steam. The equivalent energy requirement which is defined as equivalent HP steam consumption required for per kg ethanol production in Luo et al. [23] is 1.24 kW∙kg–1 for the design of Luo et al. [23], while it is 1.06 kW∙kg–1 ethanol in our optimal design, reducing by 14.3%.
The capital cost of the two designs are calculated using the economic evaluation correlations from Luo et al. [23]. The results are listed in Table 13. As seen from Table 13, the total capital cost of our design is 3.069 M$, which is 31.4% smaller than that of Luo et al. [23] (4.477 M$). The most significant reduction in the capital cost appears in the total cost of heat exchangers (i.e., HX1 and HX2), which is reduced by 0.706 M$ from 1.503 to 0.797 M$ due to much smaller area of HX1 in our design. Although one additional side reboiler is required in our design, the area of this side reboiler is much smaller compared to the decrease in the area of HX1 due to a large temperature difference (around 75 °C) in the side reboiler. Therefore, the capital cost of reboilers in our design increases only 0.123 M$.
The capital cost of the compressor in our design also decreases significantly due to much lower compression work (1041 vs. 1808 kW). Obvious reduction in the capital cost of the distillation column can also be observed in Table 13 because our column has a much smaller column diameter benefitting from a smaller reflux ratio (1.84 vs. 4.20). Note that for the capital cost evaluation, the column diameter (2.87 m) is calculated by the internals-sizing procedure built in Aspen Plus to be comparable with Luo et al. [23], while the diameter (3.33 m) shown in Fig. 12 is got from the column diameter correlation used in our optimisation model.
Conclusions
In this work, we extended the hybrid steady-state and time-relaxation feasible path optimisation algorithm developed by Ma et al. [20] for optimal design of the EDWC that is used for separation of azerotropic or close boiling-point mixtures. The tolerances-relaxation integration method for the PTC simulation from Ma et al. [20] was refined to ensure that the steady state was reached even at the initial integration phase using large tolerances. Allowing a larger tolerance for the steady-state simulation in the sensitivity analysis mode successfully enhanced the converge performance of the hybrid algorithm. Although the accuracy of Jacoain matrix and gradient informaion might be slightly sacrificed during a few optimisation steps, it is more important to enable the optimisation to continue instead of failure. The additional valid lower bounds on the vapour flow rates and purity of the recovered entrainer improved the computational efficiency as they eliminated some infeasible region.
Three case studies from the literature [21–23] have demonstrated good convergence performance of the extended hybrid algorithm, and most optimal solutions can be generated within one hour for the optimisation of EDWC, indicating the efficiency of the algorithm. For all the examples, the optimal designs obtained in current work can lead to a decrease in TAC by more than 10% compared to those reported in the literature. The importance of considering the entrainer cost in the objective function has been demonstrated in the case study 2, which obtained a design with much less entrainer make-up and achieved a significant reduction in TAC by around 12.8% compared to that without consideration of the entrainer cost. It is also found that product purity constraint was inactive in the optimal solution, indicating that fixing product purity and/or recovery to their lower bounds may eliminate the optimal design. Finally, the optimal design of the heat pump assisted EDWC was achieved for the first time by using the extended hybrid algorithm.