Single-and multi-objective optimisation of hybrid distillation-pervaporation and dividing wall column structures

The separation of azeotropic mixtures is often energy intensive, thus process intensification (PI) becomes an attractive route to enhance energy efficiency. Two of the most commonly used separation intensifications are dividing wall columns and hybrid distillation-membrane processes. In this work, three typical hybrid distillation structures, distillation followed by pervaporation (D-P), pervaporation followed by distillation (P-D), and distillation followed by pervaporation then by distillation (D-P-D), are considered and compared with a hybrid dividing wall column (H-DWC) structure, which is a highly integrated process combining a dividing wall column and a pervaporation membrane network. The four structures are compared by both single-objective and multi-objective optimisation. It is shown that the D-P-D and H-DWC structures require significantly lower total annualized costs than the other two designs due to requiring smaller membrane area, as these two structures use the membrane only to help the mixture composition cross the azeotropic point.


Introduction
Distillation is the most technologically mature, and the most widely used, separation method in the chemical industry.Conventional distillation is very powerful in separating zeotropic mixtures, however, many separation problems of industrial interest involve azeotropic mixtures, such as waterethanol, nitric acid-water, and acetone-methanol-chloroform.When the separation process involves azeotropic or close-boiling mixtures, then separation within a single conventional distillation column is impossible.Some of the most commonly used methods to separate azeotropic mixtures are azeotropic distillation, extractive distillation, and pressureswing distillation, all involving multiple columns and therefore also a high energy consumption.The increase in environmental awareness over recent decades has called for more sustainable methods for the separation of azeotropic mixtures.This is where hybrid processes enter into the picture, where in the context of this work, a hybrid process is defined as a process combining at least one distillation column with at least one membrane process in an integrated manor.We will also consider a hybrid dividing wall column and demonstrate when such a structure may be superior to more standard hybrid separation processes.Pressly and Ng (1998) introduced 15 hybrid structures for the separation of binary mixtures, including common configurations such as the distillation-pervaporation (D-P) structure, the distillation-pervaporation-distillation (D-P-D) structure, and the pervaporation-distillation (P-D) structure.In the D-P structure, a single distillation column is used to pre-separate the azeotropic mixture until close to the azeotropic point, followed by a pervaporation membrane unit to obtain the desired product at the specified purity from either the permeate or rententate stream.Alternatively, the feed can be introduced to the pervaporation membrane unit first, to overcome the azeotropic point, followed by the distillation column for further separation (P-D structure).These hybrid designs can, however, be limited by low capacity and high capital cost of the membrane unit.Instead, another design with a membrane to cross the azeotropic point placed between two distillation columns, may be more economically beneficial as the use of the membrane unit can be minimized with this design (D-P-D structure).Pressly and Ng (1998) presented a screening method based on the membrane break-even cost, i.e. the maximum allowable membrane cost calculated by subtracting the cost of the distillation column in the hybrid process from the cost of the conventional distillation process, and this was used to quickly identify if a hybrid process was likely to be economically superior compared to its conventional distillation counterpart.Koczka et al. (2007) compared the performances of the D-P, D-P-D, and P-D structures for the dehydration of ethanol using azeotropic distillation as the base case.It was found that while the D-P-D structure had a lower total annualised cost (TAC) and lower energy requirement when compared to azeotropic distillation, the D-P structure offered the most savings in terms of TAC, and the P-D structure required the least energy.Skiborowski et al. (2014) proposed an optimisation framework where the first step was to decompose the hybrid process model into individual distillation column unit(s) and individual membrane network(s).Then each of the individual units was initialised and optimised separately using a simpler objective function, such as minimising energy and minimising membrane area.The (sub-)optimal designs for each of the individual units were then "recombined" into an initial hybrid process, and more rigorous optimisation with the desired objective function was performed based on this initial design to obtain the final optimal design.In one of the case studies for ethanol dehydration, with about 42 mol% ethanol in the feed, it was shown that the D-P structure was economically favourable when compared to the D-P-D structure and to a pressure-swing distillation structure.Furthermore, in the D-P-D structure, out of the 4.765 mil € TAC, the second column accounted for only 75 k€ (equivalent to about 1.57% of TAC), which indicated that the second column in the D-P-D design was almost negligible for the separation of the ethanol/water mixture in that case study.
In addition to comparison studies of different hybrid structures, studies on individual designs have also been presented.Luyben (2009) studied the design and control of a D-P structure for ethanol dehydration (7 mol% ethanol in the feed) using a lumped membrane model (i.e. the membrane unit was split into a fixed number of fragments).Novita et al. (2018) studied a hybrid extractive distillation process combined with a membrane unit (ED-P) process for various alcohol dehydrations, including ethanol dehydration with 89 mol% ethanol in the feed.In that study, an entrainer (glycerol) was present in the hybrid process.Therefore, instead of breaking the ethanol/water azeotrope using the membrane, the entrainer was used to break the azeotrope and the membrane unit used to separate the entrainer from water, thus replacing the role of the recovery column in a conventional extractive process.Wu et al. (2020) studied the ED-P process for dehydration of n-propanol with glycerol as the entrainer.Both Novita et al. (2018) and Wu et al. (2020) showed that the ED-P process could save around 21 − 25% of the TAC when compared to a conventional extractive process.Meng et al. (2020) extended the work of Wu et al. (2020) to separate a ternary ethyl acetate/ethanol/water mixture using a ED-ED-P structure.The performance of the ED-ED-P structure was compared to a three-column extractive distillation (TCED) process, a three-column-pressure-swingdistillation (TCPSD) without heat integration, and a TCPSD with partial heat integration.It was found that the ED-ED-P structure could save nearly 61% TAC when compared to the TCPSD without heat integration, and the ED-ED-P process also had the lowest carbon dioxide emission of all four processes.In both of the studies by Wu et al. (2020) and Meng et al. (2020), the membrane unit was modelled based on a lumped model.
Other than the more commonly explored D-P structure, there are a few other hybrids structures that have been explored in the literature.Kreis and Górak (2006) reported that for a hybrid structure where the membrane unit was located at a sidedraw stream from the column, with retentate and permeate streams recycled back to the column (termed as D- side structure), this structure could assist the separation of close boiling mixtures.If the D-side structure was modified to make, for example, the retentate stream as a product stream, then the D-side structure could be used for ternary mixture separation.The D-side structure (both retentate and permeate recycled, or only one of retentate or permeate recycled) has been studied and discussed by various researchers (González and Ortiz, 2002;Koch and Sudhoff, 2013;Lee et al., 2016).For the P-D structure, where the membrane is located before the distillation column and used for preliminary separation of the feed, Zarca et al. (2018) reported that the same structure could save up to 50% total operating cost compared to conventional distillation columns for the separation of an olefin/paraffin mixture.
The increase in awareness of sustainable operation is calling for new energy-efficient and greener processes, aligned with the concept of Process Intensification (PI).Following this concept, the two columns in the D-P-D structure can be integrated into a single column with a vertical wall installed in the middle, thus forming a hybrid dividing wall column (H-DWC).Although there have been a few studies on the design of hybrid reactive dividing wall columns (e.g.Holtbruegge et al. (2015); Li et al. (2020)), the features and potentials of the H-DWC have not yet been properly explored.Therefore, in this work, an H-DWC is compared with three base cases of hybrid processes including the P-D, D-P, and D-P-D structures.For the H-DWC, we will also consider the option of having the wall extended to the top or to bottom based on the thermodynamic characteristics of the mixture.Both single objective and multiobjective optimisation will be performed, not only to compare the economic performance of each structure, but also to properly identify the various relationships between the different key variables within each design.
Note that, although this work considers pervaporation as the membrane process in the hybrid distillation-membrane process, the pervaporation membrane can easily be replaced by a vapour permeation unit by changing the membrane feed from liquid phase to vapour phase (e.g., by adding a feed heater to vaporise the membrane feed).Various studies have considered comparisons between different hybrid processes based on vapour permeation (Moganti et al., 1994;Stephan et al., 1995;Holtbruegge et al., 2015).

Methodology
To fully explore the potential of the hybrid structures, one of the key requirements is to use a proper membrane unit model which adequately captures the key characteristics of the membrane performance.In this study, pervaporation is considered due to its successful applications in industry and great potential.A rigorous membrane model will be used, the details of which will be discussed in Section 2.1 (Fig. 1).Four different hybrid processes (Fig. 2), as well as the H-DWC, will be introduced in Section 2.2.The simulations considered in the case studies are carried out in gPROMS Process (Process Systems Enterprise, 2021) with the physical properties of the mixtures considered obtained from Multiflash (KBC Advanced Technologies, 2015).The initialisation strategy for the simulations will be discussed in Section 2.3.Finally, the optimisation of all the structures will be considered in Section 2.4.

Membrane and network models
The membrane model used in this work is a parallel-flow hollow-fibre pervaporation membrane, and the shell-side and tube-side models considered are based on the models originally developed by Marriott and Sorensen (2003a).The authors considered two models with different complexities, a simpler one-dimensional model (i.e., considering only variations along the axial axis), and a more complex twodimensional model (i.e., variations along both the axial and the radial axes).Both the one-and two-dimensional models are distributed models described by partial differential and algebraic equations (PDAEs).Using PDAEs, however, may lead to problems with initialisation and long computation times.For optimisation purposes, where a large number of simulations are required, one therefore generally avoids using PDAEs.In order to simplify the membrane model for easy initialisation, many researchers have used a lumped model (i.e., a membrane model which is divided into several segments) instead of a rigorous model (Luyben, 2009;Wu et al., 2020;Meng et al., 2020), and this is the approach also taken in this work.Fig. 1 shows the membrane network used in this work, which consists of N ms number of membrane stages (i.e., the larger boxes) connected in series, and in each membrane stage n there are N mm,n number of membrane modules (i.e., the smaller boxes with diagonal lines) connected in parallel, i.e., the number of membrane modules in each membrane stage can be different.This membrane network was originally proposed by Marriott and Sorensen (2003b), and the advantage of the network is that it is simple, representative, and computationally efficient.The key assumptions are: (1) recycling streams between membrane stages are not considered as they greatly increase the computational burden during optimisation and the benefit from recycling streams between membrane stages are quite minor; (2) the feeds into the membrane modules connected in parallel are assumed to be equally distributed as there are no clear benefits of an uneven feed distribution; (3) a heater may exist in front of each of the membrane stage and their existence is optimised, and if it exists, the feed heater will always operate at the maximum allowable temperature.(This is because a higher temperature will increase the permeate flux, and thus the separation performance, so it is favourable to operate the membrane at its maximum allowable temperature (Bausa and Marquardt, 2000).The permeate streams leaving the membrane stages are in the vapour phase and initially at low vacuum pressure, however, it is preferred to collect the final permeate stream in the form of liquid for storage and transportation purposes.Therefore, after combining the permeate streams into a single (total) permeate stream, a cooler is used to cool the vapour permeate stream into saturated liquid, followed by a pump to raise the pressure back to 1 bar which is the assumed operating pressure of the columns.
The mathematical model of the membrane network is built such that there is a maximum number of N ms membrane stages, and during optimisation, the optimal number of membrane stages will be determined (N ms opt ).This is because the numerical solver in gPROMS (used in this work) required that any array size should be fixed beforehand (i.e., the array size is a parameter with a fixed value which cannot be optimised).Therefore, an efficient way to overcome this problem is to perform the calculation with N ms stages, but to collect the retentate and permeate stream information only at the N ms opt stage, which is determined from the optimisation.
This strategy has been found to be highly efficient (Chia and Sorensen, 2022).

Hybrid process models
The hybrid structures considered in this work are the distillation-pervaporation structure (D-P structure, Fig. 2a), pervaporation-distillation structure (P-D structure, Fig. 2b), the distillation-pervaporation-distillation structure (D-P-D structure, Fig. 2c), and the hybrid dividing wall column (H-DWC, Fig. 2d).The location of the membrane network will depend on the type of azeotrope, such that it is located at the distillate for a minimum boiling azeotrope, and located at the bottom for a maximum boiling azeotrope.As minimum boiling azeotropes are far more common, these are considered in this work.For the H-DWC, given the minimum boiling azeotrope, the dividing wall is extended to the bottom of the column so that the two products can be removed from either side.The structure of the membrane network used in each of the hybrids was described in Section 2.1, and includes membrane modules and stages, membrane feed heaters, a permeate cooler, and a permeate pump.It should be noted that for the feed heater into the first membrane stage in the membrane network, some a priori knowledge may be required to set a good initial value.If the temperature of the feed to the membrane network is higher than the maximum allowable temperature of the membrane module, then a distillate cooler is required rather than a heater so that the feed can be cooled down to the maximum allowable temperature.(As discussed in Section 2.1, to simplify the optimisation problem it is assumed that the membrane module will operate at its maximum allowable temperature.) Consequently, the membrane stage feed heater of the first membrane stage can then be set as non-existent (i.e., set the binary variable to zero) as the feed to the first membrane stage (i.e., feed to the membrane network) is already at the maximum allowable temperature.Otherwise, if the feed to the membrane network has a temperature lower than the maximum allowable temperature, then the design of the membrane stage feed heater will be optimised.A membrane feed pump is required to raise the pressure so that the feed (at the retentate side) will not vaporise at the membrane operating temperature.Ideally, this pressure should be optimised together with all other variables.To reduce the optimisation burden, however, this pressure will be fixed at a relatively high pressure in which the retentate stream will always remain in the liquid phase even with a large flow rate and membrane area (i.e., setting the membrane feed flow rate and membrane area at their upper bounds).This simplification is reasonable and has often been applied by other researchers (González and Ortiz, 2002;Luyben, 2009) because this pressure has a minor impact on the membrane separation due to the retentate stream remaining in the liquid phase and the physical property barely changing with the pressure change (e.g., pressure changed from 1 bar to 5 bar).
For cases where the retentate stream should be sent into a distillation column for further separation, if the retentate pressure is higher than the operating pressure of the distillation column, then the stream will be sent through a valve to reduce its pressure before entering the column.Depending on the separation system, most of the time either the retentate or the permeate stream is the required product stream, rarely both.Then, depending on the purity of the relevant stream, the stream can either be directly collected as a product stream, potentially combined with product streams of the same component from the distillation column, or the stream can be recycled back into the distillation column for further separation.In the D-P structure (Fig. 2a), the distillation column will separate the mixture into one pure component (i.e., product) stream and the azeotrope.Then, the membrane network located after the column will separate the azeotropic stream into its (almost) pure components.For the other (P-D, D-P-D, and H-DWC) structures, the role of the membrane network is to aid in overcoming the azeotropic point.In the P-D structure (Fig. 2b), if one of the product streams is from the membrane network, the membrane network will separate the feed (a combination of fresh feed and the recycle stream which is at/ near the azeotropic point) into one pure component stream (i.e., the product stream) and the other stream, which has now crossed the azeotropic point, will be sent to the distillation column for further separation to obtain the other component.In cases where both the permeate and retentate streams from the membrane network are sent to distillation column(s), the membrane network will merely need to preseparate the feed and the distillation column(s) will further separate the retentate and permeate into pure product streams (not considered in this work).In the D-P-D and H-DWC structures, the distillation columns are responsible for purifying the products (i.e., at least one distillation column each for the purification of the mixture before and after the azeotropic point) and the membrane network will only need to help to "move" the system across the azeotropic point.

Initialisation procedure and shortcut design
Proper initialisation is essential for any simulation and optimisation.Although the initialisation of a design should not affect the simulation results, a poor initial design may lead to a failed optimisation or poor local optima.A shortcut method is usually applied to obtain an initial design, which is then used as the initial design for the model initialisation.However, current available shortcut methods for hybrid distillation-pervaporation systems (Stephan et al., 1995;Pressly and Ng, 1998;Bausa and Marquardt, 2000;Caballero et al., 2009), may either be limited to certain types of structures, only estimate the minimum membrane area required without taking into account the actual membrane structure, and/or cannot provide information about recycle streams.Also, even by applying these shortcut methods, there is still a high chance that the model cannot be initialised based on the shortcut design.A trial-and-error-based shortcut procedure is therefore often applied by the user before the optimisation.A typical approach is to split the structure into its individual units, i.e., only consider the distillation column(s) or only the membrane network, then for instance use the McCabe-Thiele method or FUG method for distillation columns, and often a trial-and-error method for the membrane network, to obtain the initial design.
Taking the D-P structure as an example, for the initial design of the distillation column, and for a minimum boiling azeotrope, the distillate composition can be assumed to be close to or at the azeotropic composition, and the bottom composition can be considered at the required product composition for the heavy product.Together with feed information, the minimum number of stages and minimum reflux ratio can then be obtained from the McCabe-Thiele diagram or the FUG method.Furthermore, the distillate flow rate can also be obtained by performing simple mass balances.The distillate from the column is then used as the feed for the membrane network.As described in Section 2.1, the membrane network considered in this work is formed by N ms number of membrane stages connected in series and N mm,n number of membrane modules connected in parallel within a membrane stage n.For a proper initial membrane network design, especially for optimisation purposes, the initial value in the shortcut design of the membrane network is set as N ms , i.e, the maximum/upper bound of the number of membrane stages.As for the number of membrane modules in a membrane stage n, N mm,n , these can be set as a middle value between the lower and upper bounds.In addition, initially, all heaters are assumed present (i.e., the binary variables for the existence of membrane feed stage heaters are initially all set to one).With these initial specifications, the design should achieve good product purities (i.e., close to the product requirements).If not, then the upper bounds of the membrane variables such as N ms and N mm,n can be increased to provide a larger search space for optimisation.
For the other three structures, the membrane network is used for crossing the azeotropic point, which means that the composition of the retentate stream cannot be determined directly.Therefore, a reasonable composition of the retentate stream can be estimated to be, e.g. the middle point between the pure component and the azeotropic point.Then, the same procedure as described above can be applied to initialise the distillation column(s) and membrane network.
In the presence of recycle streams, the recycle stream should initially be torn (i.e., the recycle stream is not recycled and a pseudo-recycle stream is added as an additional feed).The convergence between the recycle stream and pseudorecycle stream can slowly be achieved, starting with the assumption that there is no recycle (i.e., the recycle stream is directly removed from the system and the pseudo-recycle stream has zero flow rate) and then the pseudo-recycle flow rate can be gradually increased until it is the same as the recycle flow rate and the recycle loop is then closed.

Optimisation
The optimisation of the hybrid structures involves both continuous (e.g., reflux ratio, distillate and bottom flowrates) and discrete (e.g., total number of distillation stages, feed location, number of membrane stages, and number of membrane modules in each membrane stages) optimisation variables.Together with the highly non-linear total annualised cost (TAC) as the objective function, the optimisation becomes a typical Mixed Integer Non-linear Programming (MINLP) problem, which can be handled by either deterministic optimisation or stochastic optimisation.Due to the high complexity of the designs considered in this work, the built-in Outer Approximation / Equality Relaxation / Augmented Penalty (OAERAP) MINLP optimiser in gPROMS may not easily converge into a good local optima (Chia et al., 2021).Also, this optimiser is not able to perform multi-objective optimisation.Stochastic optimisation is therefore used in this work for both single objective optimisation and multi-objective optimisation.Different stochastic optimisation methods (e.g., Genetic Algorithm (GA), particle swarm optimisation (PSO) and simulated annealing (SA)) have shown great success for optimisation of chemical processes due to their ability to escape poor local optima (owing to the randomness involved), robustness to different structures, no need for providing initial values for each optimisation variable, and capability of performing multi-objective optimisations (Waltermann and Skiborowski, 2017;Yang et al., 2019;Duanmu et al., 2022a).
As gPROMS has no built-in stochastic optimisation method available for solving MINLP problems, we will consider stochastic optimisation using Genetic Algorithm (GA) for single objective optimisation, and the fast and elitist nondominated sorting Genetic Algorithm (NSGA-II) developed by Deb et al. (2002) for multi-objective optimisation.The optimisation is developed in C#, and is connected with gPROMS via the Foreign Process Interface (FPI, a communication protocol used by gPROMS coded in C++) and gO:Run (Process Systems Enterprise 2022), a gPROMS execution-only engine to start simulations in gPROMS externally).Details of the optimisation methods applied are discussed in the respective subsections.

Single objective optimisation
In this work, the single objective optimisation is carried out using Genetic Algorithm (GA) on a workstation with a dual Intel Xeon Gold 6226R CPU (32 cores and 64 processors in total) with a clock speed of 2.90 GHz.Moreover, the total memory capacity is 192 GB with speed of 3200 MHz.To enhance the speed of the optimisation, 40 processors are used for parallel computing.In our experience, the optimisation can then be approximately 20-30 times faster than singlecore optimisation depending on the complexity of designs.A timeout function (20 s in this work) is applied to avoid any "infinitely long" simulations as a few GA chromosomes with poor values of optimisation variables may lead to slow simulations (i.e., hard and/or slow to converge and long time for reinitilisation).Once this happens, this chromosome will be considered as an infeasible simulation and the corresponding thread will be terminated and restarted for the next simulation.As an additional feature, a dynamic bound function is also introduced in the optimisation to avoid unrealistic solutions.For example, for a distillation column, the feed location should always be smaller than the total number of stages of the column.Once this rule is violated during optimisation, the feed location will be regenerated to a random integer within the bounds between its lower bound and current value of total number of stages.
The settings of the GA can be found in Chia et al. (2021).The population size is chosen as five times the number of optimisation variables in each design.The stopping criteria is either when the GA reaches a maximum of 200 generations (i.e., MaxGeneration = 200) or when the objective function stays constant within a given tolerance for 20 consecutive generations (i.e., MaxStallGeneration = 20), depending on whichever comes first.The fitness and constraint tolerances are set as 10 −4 .Elitism is applied in the GA, where the top 10% of the chromosomes are allowed to participate in the next generation directly without going through crossover and mutations.The top 50% of the chromosomes (including the elites) are chosen as the parents for the next generation, where two chromosomes (out of the top 50%) are randomly selected and discrete crossover is applied to produce one child.The selection and crossover procedures are repeated until the remaining 90% of the population is filled up (as the elites took up 10% of the whole population).The mutation probability is initially set as 50% when the simulations are offspec (i.e., at least one constraint is not met), and is decreased to 10% when at least one simulation is on-spec (i.e., all constraints are met).Moreover, to avoid guessing the R value for the penalty of fitness, the penalty function proposed by Deb (2000) is used, where the penalised fitness equals the summation of the worst solution (i.e., the largest fitness for a minimisation task) and the total constraint violations.
To ensure a good optimal design, the very first few optimisations will be performed with wide ranges of each optimisation variables (e.g., the range of the total number of stages in a distillation column is from 1 to 60 in the case studies) to obtain a preliminary optimisation.Then the optimisation will be carried out with a narrower optimisation variable bounds to improve the quality of optimal designs.

Multi-objective optimisation
The multi-objective optimisation method considered in this work is the method proposed by Deb et al. (2002).NSGA-II is in essence a sorting approach for determining the Pareto front in a fast way (i.e., with fewer calculations).Key functions such as selection, crossover, and mutation need to be decided by the user.In this work, the binary tournament selection (Back et al., 2000), discrete crossover (Umbarkar and Sheth, 2015), and uniform mutation (Soni and Kumar, 2014) are applied.The NSGA-II with the selected functions was tested with a few examples presented by Deb et al. (2002), and clear and accurate Pareto fronts indicated a good performance of the applied NSGA-II with the selected functions (not shown).
In this work, as the main focus is the economic performances of each design, the two objective functions chosen to illustrate the multi-objective optimisation are the total capital cost and the total operating cost.Thereby, the relationship between these two costs, and the relationships between the key optimisation variables and costs, can be reflected.In this work, both objectives are related to costs, however, it is also possible to, say, consider cost versus environmental impact in the form of energy consumption.It should be noted that the parallel computing, timeout function and dynamic bound introduced in Section 2.4.1 are also applied for NSGA-II for better performance.
Due to the randomness involved in both GA and NSGA-II, it is necessary to repeat the optimisation of each design a few times (at least five times is used in this work) to get a good local-optimal design for comparison.Especially for NSGA-II, the Pareto front of a single optimisation may not be clear enough, therefore, Pareto front results from repeated optimisations (five repeats in this work) are combined and re-processed to find the final Pareto front in the following case studies.

Case studies
The case studies will consider the performance of the hybrid dividing wall column (H-DWC) against those of the other structures (see Fig. 2), namely the distillation-pervaporation (D-P) structure, the pervaporation-distillation (P-D) structure, and the distillation-pervaporation-distillation (D-P-D) structure.The separation task considered is the separation of the minimum boiling azeotropic ethyl acetate/ethanol mixture, with a feed composition of 0.2∕0.8molmol −1 and a flow rate of 200 kmol h −1 to obtain at least 99 mol% of ethyl acetate and ethanol in their respective product streams.The feed is assumed to be provided as a saturated liquid at 1 bar.The physical properties of the liquid can be described by the UNIQUAC model, while the vapour phase is assumed to be an ideal gas.Due to the lack of openly available membrane information for this mixture, it is assumed in the calculations that the membrane characteristics reported by Tsuyumoto et al. (1997) for the separation of an ethanol/water mixture are applicable also for this mixture (see Table A1 for validation results), and that the ethyl acetate will mainly leave from the retentate, while ethanol will permeate through the membrane and mainly leave from the permeate.The ethanol/water mixture is not considered because beyond the ethanol/water azeotropic point, the separation of ethanol from water using a distillation column is still very difficult as the bubble points of ethanol and water are still very close to each other, thus making the H-DWC structure unsuitable.
The design of the hybrid processes was discussed in Section 2.1 (membrane network) and Section 2.2 (integrating membrane with distillation column).The H-DWC is considered as the energetically equivalent 2-column Petluk arrangement, with a prefractionator and a main column.For the costing, the column diameter of the H-DWC is recalculated by considering the cross-sectional area as the summation of the cross-sectional areas of the prefractionator and the main column (Duanmu et al., 2022b).The height of the H-DWC is chosen as the highest column height of the prefractionator and the main column.
For the maximum allowable membrane temperature, a temperature of 70 ∘ C is considered as the maximum temperature the membrane can withstand.The pressure of the feed stream into the first membrane stage is specified as 5 bar where the retentate stream will remain in liquid phase even with heating and pressure drop (as explained in Section 2.2), and the membrane permeate side pressure is chosen to be 400 Pa (Tsuyumoto et al., 1997).Further details of the membrane system can be found in Appendix A.
For the optimisation, it is assumed that there is no pressure drop in the distillation column so the columns are operated at 1 bar throughout.The membrane network will be located at the distillate end of the column since the azeotropic system is a minimum boiling azeotrope.Because the distillate temperature will always be higher than 70 ∘ C (as the bubble point of the ethyl acetate/ethanol azeotrope is around 71 ∘ C), an additional distillate cooler is required in addition to the total condenser considered in the column model, and the first membrane stage feed heater is therefore removed (i.e., by fixing its binary variable to zero) as previously discussed.The permeate consists of almost only ethanol, hence the permeate stream can for some structures be combined with the bottom stream from the distillation column and collected as a single ethanol product stream.The initialisation and optimisation procedures used (together with the CPU specifications) were previously outlined in Sections 2.3 and 2.4, respectively.
In this work, the objective function for single-objective optimisation is considered to be the total annualised cost (TAC) of the structure, and the objective functions for multiobjective optimisation are the capital cost (CAPEX) and operating cost (OPEX).The TAC is the summation of the annualised CAPEX and OPEX with the annual operating hours assumed to be 8400 h y −1 , eight years of plant life, and a membrane lifetime of two years (Duanmu et al., 2022b).The CAPEX includes the cost of the column shell, trays, condenser, reboiler, membranes, heaters/coolers, and pumps.The distillation column sizing and costing equations can be found in Duanmu et al. (2022b), while the costing equations for the membrane and other pieces of equipment (e.g., pumps, heaters, coolers) can be found in Appendix B. The OPEX includes the heating, cooling and electricity costs.The type of heating or cooling utility is chosen automatically based on the reboiler/heater or condenser/cooler temperatures, and the utility costs can be found in Turton et al. (2012).Details of the cost information can also be found in Appendix B.

Single-objective optimisation
Table 1 shows the key variables for the optimal design of each structure.For the H-DWC, column C1 is the prefractionator and column C2 is the main column of the Petlyuk arrangement as discussed above.The results show that the D-P, P-D, and H-DWC structures require a similar total number of distillation stages, 24, 21, and 25 stages (recall that H-DWC has only one shell and the main column is column C2), respectively, which is much smaller than the D-P-D structure (49 stages in total) as it contains two distillation columns.In this work, it is assumed that the tray spacing is fixed and constant, thus the number of stages on both sides of the wall in the H-DWC has to be the same.It should be noted that the number of stages on either side of the wall can be different if the tray spacing is assumed to be different, or if the column internals installed on either side of the walls are different.
Looking at the designs where the membrane is located after a distillation column (i.e., all designs except the P-D structure), for the total reboiler and condenser duties (summation of reboiler/condenser duties of both columns C1 and C2), it can be seen that the D-P structure requires less total reboiler and condenser duties (2107 kW for D-P vs 2842 kW for D-P-D and 2704 kW for H-DWC for reboiler duty, and 2108 kW kW for D-P vs 2779 kW for D-P-D and 2620 kW for H-DWC for condenser duty).For the P-D structure, the total reboiler and condenser duties (989 kW for reboiler duty and 826 kW for condenser duty) are significantly lower than the others as the feed flow rate (i.e., the retentate stream in P-D structure) into the column is smaller than for the other designs (112 kmol h −1 for P-D structure and 200 kmol h −1 for the other designs).However, the energy is not saved but is instead now required by the membrane system.Moving on to the energy required by the D-P-D and H-DWC structures, it is not surprising that the H-DWC does not save much energy when compared to the D-P-D structure (heating duties in the D-P-D and H-DWC structures are 3305 kW and 3238 kW, respectively; and cooling duties are 3375 kW and 3313 kW, respectively).This is consistent with previous findings that a standard DWC with the wall extended to one end, as considered in this work, will have a similar energy requirement as its conventional counterpart (Kaibel, 2014).
Focusing on the membrane network, it is noted that the P-D structure has the lowest column reboiler duty, but the heaters in its membrane network consumed much more energy when compared to the other structures (e.g., about 3.6 times higher than for the D-P-D structure) due to the fact that the flow rate into the membrane network is much higher in the P-D structure (e.g., 272 kmol h −1 for the P-D structure and 115 kmol h −1 for the D-P-D structure), as most of the ethanol is separated from the system by the membrane network in the P-D structure.
It can also be seen from the results that the D-P-D and H-DWC structures require fewer membrane stages and modules, thus smaller membrane areas (up to 60% reduction), when compared to the other structures because, unlike the membrane network in the P-D and D-P structures where the network is used for final purification of products to meet the product specifications, the membrane networks in the D-P-D and H-DWC structures are only used to help cross the azeotropic point.Upon closer inspection, the D-P-D structure requires slightly less membrane area than the H-DWC structure (522 m 2 for D-P-D vs 600 m 2 for H-DWC).
For the economic comparison, the detailed dimension and cost equations can be found in our previous work (Duanmu et al., 2022b).A stacked bar graph (Fig. 3) and individual donut charts (Fig. 4) are provided to illustrate the economic performances of each individual design.Overall, the D-P-D structure has the lowest total annualised cost (TAC) of $ 2.52 M y −1 , followed by the H-DWC structure (about 2% higher, at $ 2.56 M y −1 ).Although the D-P structure has a lower operating cost ($ 1.62 M y −1 ) than the D-P-D ($ 1.75 M y −1 ) and H-DWC ($ 1.77 M y −1 ) structures, the larger membrane area required by the D-P structure leads to a significant increase in capital cost, which eventually results in a 12% increase of TAC compared with the D-P-D structure ($ 2.81 M y −1 for D-P vs $ 2.52 M y −1 for D-P-D).Not surprisingly, the P-D structure has the highest TAC (at $ 3.66 M y −1 , 45% higher compared to the D-P-D structure) as the entire feed stream plus the recycled material goes through the membrane, thus demanding not only a large membrane system but also a much higher permeate cooling duty due to the need for condensing the low-pressure vapour (400 Pa) in the permeate using expensive refrigerant.It should be noted that although the membrane system in the P-D structure handles much more feed material compared to the other designs, the overall membrane area required is similar to that of the D-P structure, which makes sense as the task of the membrane system in the P-D structure is to cross the azeotropic point but not to separate the mixture completely (i.e., to achieve 99 mol% of each product).However, all the ethanol leaves the system as the product from the membrane system in the P-D structure, which requires a much higher permeate cooling duty than the D-P structure.Moreover, the results clearly show that for each design, the biggest contribution to the TAC is still the operating cost (ranging from about 58% for the D-P structure to 69% for the D-P-D and H-DWC structures).The heating cost contributes the most, not only in terms of the total operating cost, but also in the TAC (up to about 56%).For the capital cost, the D-P-D and H-DWC structures have a lower cost (up to 31% reduction compared to the most expensive P-D structure) due to smaller demand for the expensive membrane (e.g., the D-P-D structure requires about 60% less membrane area compared to the P-D structure).It should be noted that for the D-P-D and H-DWC structures, the capital cost of the distillation column (summation of the column shell and trays, condenser, and reboiler costs) is almost doubled compared to those of the other two designs due to the additional distillation column, i.e. column section, introduced in these two structures.
A further study considering the membrane cost being reduced to 50% of the base case was also performed to consider the impact of the membrane cost, as the membrane has been found to be one of the key factors affecting the costs.The comparison was made just between the D-P-D and H-DWC structures as these were identified as the most promising structures.Both designs were re-optimised with the new membrane cost and the results (detailed designs and results are not shown) show that the H-DWC structure (TAC as $ 2.25 M y −1 ) is now marginally cheaper than the D-P-D structure (TAC as $ 2.26 M y −1 ), which indicates that the D-P-D and H-DWC structures have close economic performances also for different membrane prices.

Multi-objective optimisation
In the multi-objective optimisation, the bounds of the optimisation variables for each structure are narrower (e.g., the bounds of the membrane stages in the H-DWC structure is changed from 1-10 to 4-8) based on their corresponding optimal designs found in the single-objective optimisation (Section 4).This is because it is unnecessary to have a large search space to obtain a clear Pareto front and to establish the relationships between the key variables given the singleobjective insight.Moreover, with wide boundaries, the simulation time is significantly longer.Even with parallel computing (e.g., 50 processors used in this work), it may take a day to obtain one single multi-objective optimisation with wide bounds compared to for example the corresponding optimisation time of about four hours for H-DWC using the narrower bounds.The boundaries of key optimisation variables can be seen from the x-and y-axes from the correlograms (Figs. 6,8,10 and 12).As described in Section 2.4.2, the multi-objective optimisation is performed using NSGA-II (Deb et al., 2002).The population size and maximum generations of each optimisation task is set as 150 and 350, respectively.Before considering the results, a few terms should be defined: .
• Infeasible simulations: Simulations which cannot be converged due to either calculation failure or timeout during optimisation (i.e., unsuccessful simulations).These solutions are not shown on the scatter plots.
• Feasible simulations: Successful simulations, i.e. simula- tions that have not failed, which may either be within or outside the constraints.These are shown on the scatter plots as dark/blue (on-spec) or light/grey (off-spec) points.
• On-spec simulations: Successful simulations which are within the constraints (e.g., achieved the product purities).These are shown on the scatter plots as dark/blue points.
• Off-spec simulation: Successful simulations which are outside the constraints (e.g., did not meet the product purities).These are shown in the scatter plots as light/grey points.
• Pareto front: A set of solutions where the objective func- tion value of a solution cannot be further improved without compensating the objective function value of the other solution (Ngatchou et al., 2005).
• Accumulated Pareto front: A collection of Pareto fronts from all the repeated optimisations (each optimisation is repeated five times in this work).These are shown on scatter plots as larger light/orange points.
• Final Pareto front: The "true" Pareto front is a set of the non-dominated solutions from the accumulated Pareto front.These are shown on the scatter plots as larger dark/ green points.
The results from the multi-objective optimisation are shown in Figs. 5, 7, 9 and 11 for the distillation-pervaporation (D-P), pervaporation-distillation (P-D), distillation-pervaporation-distillation (D-P-D), and hybrid dividing wall column (H-DWC) structures, respectively.The corresponding correlograms (Figs. 6,8,10 and 12) show the correlations between the key continuous optimisation variables for each structure.It should be noted that all Pareto fronts (orange and green points in the scatter plot, e.g., on Fig. 5a), are gathered from all the repeated optimisations (each optimisation is repeated five times in this work).The final Pareto front (green points) is the set of non-dominated solutions by re-analysing all the five Pareto fronts.The feasible points (light/grey and dark/blue points in the scatter plots) are also gathered from repeated optimisations.Due to the large data size of the optimisation (45,000 total simulations for each optimisation (on-or off-spec), i.e., 225,000 simulations for five repeated optimisations), plotting all points together may lead to a very large figure size with minor improvement in the visualisation of the results.Therefore, between one to five set(s) of optimisation results from D-P, P-D, and D-P-D are selected which correspond to 15,000-20,000 on-spec simulations in total (e.g., for D-P a random combination of two out of five sets of the optimisation results sufficed for visualisation purposes as that combination gave about 15,489 on-spec simulations).It should be noted that the selection of optimisation results is randomly selected (i.e., the optimisation results are randomly picked without considering the optimal values).The way of selecting which results to include in the plots will not affect the findings and conclusions as the NSGA-II optimisation method can ensure a good searching space.For H-DWC, there are only about 7000 onspec simulations even when all five sets of optimisation sets are considered because this structure is difficult to converge.However, this does not affect the findings and conclusions, as 7000 points are still large enough to have a clear plot.Since the points in the Pareto fronts are significantly fewer, the Pareto fronts from all five sets of optimisation results are therefore plotted.Similarly, the density plot for the on-spec simulations (e.g., Fig. 5b) is processed according to the points in its corresponding scatter plot.
For each structure, the total membrane area and one additional key variable are chosen when plotting scatter plots for CAPEX (e.g., Fig. 5c) and OPEX (e.g., Fig. 5e) to show their (statistical/probability) correlation.The R 2 values are also shown in the figures to indicate the strength of the correlations.The correlations for CAPEX are defined below (same definitions apply to the correlations for OPEX): .
• Positive correlation with CAPEX: As the value of the vari- able considered (i.e., the x-axis) increases, it is more likely for the possible designs (on-spec/off-spec) to be found at a larger CAPEX.
• Negative correlation with CAPEX: As the value of the variable considered (i.e., the x-axis) increases, it is more likely for the possible designs (on-spec/off-spec) to be found at a smaller CAPEX.
• No clear correlation with CAPEX: As the value of the variable considered (i.e., the x-axis) changes (i.e., increases or decreases), the range of CAPEX for the possible designs (on-spec/off-spec) barely changes.
In the correlograms (e.g., Fig. 6), the density plots on the diagonals show the distribution of the variables that are on their corresponding x-axis (i.e., the y-axis is the number of occurrence while the x-axis is the value of the corresponding variable).The scatter plots (i.e., plots that are not on the diagonal) show the correlation between the pair of variables on its corresponding axes.In the following, each hybrid structure will be considered in turn.

Distillation-Pervaporation (D-P) structure
Starting with the D-P structure, Fig. 5a shows all the feasible simulations obtained from the optimisation, as well as the accumulated and final Pareto fronts.The overall shape is diamond-like, which may be caused by the bounds of the optimisation variables used in the optimisation task.By looking at the detailed results for the off-spec simulations (light/grey points) near the bottom-left corner (not shown), the designs of the corresponding distillation columns have no clear trend but for the membrane networks, the number of membrane stages is most likely to be at the lower bound value (6 membrane stages).However, for the on-spec simulations (dark/blue points), the designs with lower CAPEX and OPEX are most likely to have membrane stages ranging between 8 and 10 stages, and this finding can be explained by looking at Fig. 6.The third column (x-axis with the total membrane area) shows that to achieve the on-spec simulations in the D-P structure, a minimum membrane area of about 1200 m 2 is required.Therefore, a larger membrane size (i.e., more membrane stages) is more commonly seen for the on-spec simulations.The density plot of on-spec simulations (Fig. 5b) shows that most of the designs will result in CAPEX between $ 11 M to $ 13 M and OPEX between $ 1.70 M y −1 to $ 1.95 M y −1 .
The scatter plot for CAPEX vs total membrane area (Fig. 5c) shows clearly that the total membrane area has a strong correlation with CAPEX (i.e., with the increase in membrane area, the CAPEX range of the possible designs is also increasing), which indicates that the membrane design is the most important variable affecting the CAPEX.This makes sense as from the single-objective optimisation (Fig. 4 in Section 3), the membrane contributes to most of the CAPEX for each structure.From Fig. 5d, there is a very weak but positive correlation between distillate flow rate and CAPEX (solid line), which further illustrates that the total membrane area plays the most important role in CAPEX.It is interesting to note that the scatter plot between total membrane area and OPEX (Fig. 5e) has a similar shape to that of the scatter plot between CAPEX and OPEX (Fig. 5a), which makes sense as the total membrane area has a strong correlation with CAPEX, therefore the x-axis in Fig. 5a (which is the CAPEX) can be closely represented by the total membrane area, making Fig. 5a and e look similar.
Moving on to the distillate, Fig. 5f shows a moderate positive correlation with OPEX for both off-spec and on-spec simulations (dashed and solid lines, respectively), indicating that an increase in the distillate flow rate will most likely increase the OPEX.This is because, as the distillate flow rate (D) increases, regardless of the value of the reflux ratio (RR), the vapour flow rate in the column will increase as V = D × (1 + RR).An increase in the vapour flow rate means that there is a higher demand for reboiler duty.Also, as the distillate flow rate increases, the membrane network needs to handle higher throughput, thus the operating costs related to the membrane network will also increase.It is, however, interesting to note that the OPEX bottom line of the on-spec simulation (i.e., the bottom edge of the dark/blue points) in Fig. 5f barely changes.From Fig. 6, it can be seen that for the on-spec simulations, as the distillate flow rate increases, the minimum required reflux ratio clearly decreases.Taking the extreme points for the distillate (i.e., at the highest and lowest distillate flow rate) and the corresponding reflux ratios, it can be calculated that the vapour flow rate (recall that V = D × (1 + RR)) ranges from about 100 × (1 + 0.8) = 180 kmol h −1 to 80 × (1 + 1.8) = 224 kmol h −1 , which is significantly lower when considering the range for the whole search space (including off-spec simulations) which is about 80 × (1 + 0.8) = 144 kmol h −1 to 100 × (1 + 2.0) = 300 kmol h −1 .Therefore, the bottom edge of the OPEX for the on-spec simulations in Fig. 5f almost levels off due to the cancelling effect between the distillate flow rate and the reflux ratio, which in turn leads to an almost negligible increase in the vapour flow rate, thus the reboiler duty remains relatively constant.Moreover, Fig. 4a showed that the heating costs (mainly from the reboiler duty) are the main contributor to OPEX (about 70%), so a relatively constant reboiler duty will implicitly mean a minor change in OPEX as the distillate flow rate increases.
From the density plots in Fig. 6 (diagonal plots), it can be seen that for the three variables considered, in general, the distribution of the on-spec and off-spec simulations superimpose, except when the total membrane area is below around 1200 m 2 .This independent area (i.e., when total membrane area is less than 1200 m 2 ) indicates the range of total membrane area where the simulations are always offspec, i.e., do not meet the specifications.Conversely, the area where the on-spec and off-spec simulations overlap means that the change in the variable has no definite impact on achieving the design specifications of the simulations.The cutting point for the membrane area is expected, as for this structure, the membrane network is entirely responsible for the purification of the lighter component (ethyl acetate) until it achieves the product specification.Thus the membrane area required will be larger than those required by the D-P-D and H-DWC structures, where the membrane is used only to overcome the azeotropic point.
Although there is no clear indication that the simulations will become on-spec or off-spec beyond or upon a certain value for the distillate flow rate or the reflux ratio, it is clear that as their values increase, the number of on-spec simulations increases, meaning that there is a higher chance for an on-spec simulation at higher distillate and reflux ratio when compared to lower distillate and reflux ratio.For the total membrane area, the distributions for the off-spec and on-spec simulations are normal (slightly skewed) with a mean of around 1200 m 2 (which corresponds to the cutting point of off-spec and on-spec simulations) and 1800 m 2 , respectively.

Pervaporation-Distillation (P-D) structure
Moving to the P-D structure, the overall shape of all feasible simulations points (shown in Fig. 7a) looks like an eagle beak, and the final Pareto front (green points) is within the search space, indicating a well-defined search range of the optimisation variables.A unique point of this scatter plot is that there is an unfilled gap (i.e., infeasible simulations) at the top-right corner, which is not seen in the scatter plots of the other structures.The same gap can also be seen in Fig. 8 for the pair plots of distillate flow rate and total membrane area (top-right and bottom-left plots), where the gap is formed when the distillate flow rate is greater than about 125 kmol h −1 and total membrane area is about 2500 m 2 .This may indicate that at large distillate flow rate and membrane area, the simulations become infeasible (as infeasible simulations are not plotted on the graphs).However, theoretically speaking, an increase in distillate flow rate (thus an increase in the recycle flow rate back to the membrane network) will require a larger total membrane area to process the higher throughput, so the top-right corner should contain some feasible simulations.Upon further investigation, it was found that the infeasible simulations are caused by initialisation failure because when another set of initial values near the gap (i.e., around 125 kmol h −1 of distillate and 2500 m 2 of total membrane area) is used, the simulations run successfully.This is a known challenge for hybrid processes as the integration of the membrane network greatly increases the complexity of the whole structure, making the initialisation of the structure a difficult task.Therefore, for cases where the Pareto front is affected by initialisation failure, the simulations may need to be repeated with different sets of initial values.Nevertheless, since the Pareto front (i.e., lower CAPEX and OPEX) in this work is on the opposite side of the gap (i.e., higher CAPEX and OPEX), the infeasible simulations due to failed initialisation do not affect the findings.Fig. 7b shows that most of the on-spec simulations will fall between a narrower range of CAPEX between $ 12 M to $ 16 M and a broader range of OPEX between $ 3 M y −1 to $ 3.5 M y −1 .From Fig. 7c and d, it can be seen that, once again, the total membrane area has a strong (i.e., R 2 = 0.8884) positive correlation with CAPEX, while the distillate flow rate has no clear (i.e., R 2 = 0.0033) correlation with CAPEX.For the correlation with OPEX, Fig. 7e and f show that the total membrane area has no clear correlation with OPEX while reflux ratio has a moderate positive correlation with OPEX, respectively.The finding where Fig. 7e is similar to Fig. 7a, and the less steep OPEX bottom line for the on-spec design in Fig. 7f, are once again observed.However, the reason behind this is different from the explanation for this phenomenon in the D-P structure.In this (P-D) structure, from the results with single-objective optimisation (Fig. 4b), it can be seen that the energy consumption in the membrane system contributes the most to the OPEX.For a small reflux ratio, the energy consumption in the distillation column would be small, however, to achieve the design specification at the outlet of the column, the purity of the lightest component in the retentate stream entering the distillation column will need to be fairly pure, thus leading to a higher energy requirement in the membrane system.For a large reflux ratio, the condition is reversed.Therefore, the OPEX is balanced with the change in the reflux ratio.
Moving on to the correlogram shown in Fig. 8, it can be seen that most of the designs, regardless of the values of the optimisation variables considered, have a higher chance of yielding an on-spec design (i.e., the distribution of the onspec simulations for each variable almost always overshadows the off-spec simulations).The distributions of the off-spec and on-spec designs for the distillate flow rate (top left plot) and reflux ratio (middle plot) completely overlap, indicating that the simulation can be either off-spec or onspec for the whole range of distillate flow rate and reflux ratio investigated.For the reflux ratio, as its value increases, there is a higher chance for the simulations to be on-spec.It can also be observed that there is a clear cut point for the total membrane area, where in order to obtain an on-spec design, the total membrane area has to be greater than 1000 m 2 (bottom right plot).Moreover, the total membrane area for the on-spec simulations is almost normally distributed around 2000 m 2 .

Distillation-Pervaporation-Distillation (D-P-D) structure
The overall shape (see Fig. 9a) of the scatter plot for all feasible D-P-D simulations is droplet-like.The density plot of on-spec simulations (Fig. 9b) shows that most of the designs will result in CAPEX between $ 9 M to $ 12 M and OPEX between $ 2 M y −1 to $ 2.6 M y −1 , which is spread more evenly compared with the other structures.
Considering the scatter plots for CAPEX (Fig. 9c and d), the membrane design (i.e., total membrane area) still plays the most important role in CAPEX, as expected.The next most important variable, the distillate flow rate in column C1, has no clear correlation with CAPEX (i.e., very small R 2 value).For OPEX (Fig. 9e and f), the membrane design shows no clear correlations with OPEX.The distillate flow rate shows a moderate positive correlation with OPEX in general.For onspec simulations, the horizontal line formed by designs with minimum OPEX can be explained the same way as for the D-P structure in Section 5.1 because it also shows a negative correlation with reflux ratio in column C1 as shown by the clear division line for their pair plots in Fig. 10.Moreover, the distillate flow rate and reflux ratio in Column C1 have a stronger effect on OPEX than was seen previously (Section 3), where for the D-P-D structure, the reboiler and condenser duties of column C1 were much higher than the corresponding duties of column C2 for single-objective optimisation.
The correlogram (Fig. 10) shows the relationships between all continuous variables for the D-P-D structure, including distillate flow rate and molar reflux ratio in both columns and the total membrane area.For the distillate flow rate and reflux ratio of column C1, it can be seen that the distribution is heavily skewed towards the right, indicating that there is a higher chance for the simulation to be on-spec when their values are higher.Moreover, it is found that there is a clear minimum requirement for distillate flow rate in column C1 (about 80 kmol h −1 ) as, with too low a distillate flow rate in the first column, some of the light components will remain in the bottom stream leading to an impure product of the heavy component, i.e., an off-spec simulation.Since the azeotropic point is known from the vapour-liquid equilibrium diagram, together with the feed information, the minimum requirement of the distillate can be roughly calculated (not shown).
The distillate and reflux ratio distribution in column C2 is relatively uniform compared to the other variables, with a slightly higher chance of obtaining an on-spec simulation at a larger distillate flow rate and reflux ratio (of column C2) than when their values are lower.The total membrane area has a normal distribution at around 1250 m 2 for the on-spec simulations.A minimum of about 500 m 2 total membrane area is required to achieve an on-spec simulation, which makes sense as a very small membrane area may not separate enough heavy components from the mixture to fulfil the mass balance of the heavy component (recall that the heavy component can only leave the system with the bottom stream from column C1 or from the membrane).Other than the relationship explained for the distillate flow rate and reflux ratio in column C1, there are no clear relationships between any other two variables.Although there may be a similar trend for distillate flow rate and molar reflux ratio in column C2, the search space is not large enough to show the complete relation (but the current search space suffices to find the Pareto front).

Hybrid Dividing Wall Column (H-DWC) structure
Fig. 11a shows the scatter plot of the distribution of all the feasible simulations in terms of CAPEX and OPEX for the H-DWC structure.It can be seen that the overall shape of the distribution is also droplet-like.As for the other structures, the CAPEX is strongly correlated with the total membrane area, thus the same finding as in Section 5.3, where there will always be a design which can yield a low OPEX ($ 1.7 M y −1 ) regardless of the membrane design, can also be found for this structure.Fig. 11b shows that the CAPEX and OPEX are concentrated in the range of $ 7 M to $ 10 M and $ 1.9 M y −1 to $ 2.4 M y −1 , respectively.
Fig. 11c shows that the total membrane area has a strong (i.e., R 2 = 0.9356) positive correlation with CAPEX, and Fig. 11d shows that there is a weak (i.e., R 2 = 0.0841) negative correlation between the bottom flow rate of column C1 and CAPEX.There is no clear correlation between total membrane area and OPEX, and similar shapes and distribution of the simulations in Fig. 11a and e were discussed in Section 5.1.
The reflux ratio of column C2 shows a moderately positive correlation with OPEX (see Fig. 11f) and, unlike the relatively flat bottom line (i.e., the edge of the dark/blue points) of the OPEX in the other designs, the bottom line of the OPEX for H-DWC increases linearly after about 1.5 mol mol −1 , however, it is relatively flat before that.The reason for this is not straightforward, even with the help of the correlogram (Fig. 12).Although the variable is named as the reflux ratio of column C2, it is actually the reflux ratio for the whole column (recall that H-DWC is a single column with one condenser and two reboilers but is modelled as a Petlyuk structure with two columns).Thus, the variables for the H-DWC structure may have even stronger interactions than the D-P-D structure.It should be noted that the R 2 value for the reflux ratio of column C2 is much higher than the R 2 values of the variable that has a stronger correlation with OPEX in the other structures (i.e., the variables presented in Figs.5f, 7f and 9f).
From Fig. 12 it can be seen that some of the optimisation variables are correlated.As the bottom flow rate of column C1 increases (first column plots), the minimum sideflow of column C2 (i.e., the liquid thermal coupling stream flow rate back into column C1, see Fig. 2d) required for the simulation to become on-spec increases linearly (as shown by the welldefined line between the off-spec and on-spec simulations).This is expected as the feed flow rate into column C1 is fixed, and an increase in the bottom flow rate means that the other inlet stream to column C1, which is the side stream from column C2 (i.e., the liquid thermal coupling stream), will increase for the mass balance to be closed.It should also be noted that the sideflow from column C2 into column C1 serves as the reflux flow of column C1.Thus it is not surprising to see that the bottom flow rate in column C1 and the minimum sideflow from column C2 have a strong correlation, as seen in the distillate flow rate and reflux ratio pairs for the other structures.
As the bottom flow rate of column C1 increases, in general, the total membrane area required for an on-spec simulation to happen decreases (Fig. 12, bottom left plot).The bottom of column C1 removes the heavy component (ethanol) from the system, so a larger bottom flow rate in column C1 means less membrane area is required to remove the heavy component.The same reason can be applied to explain the decrease in the minimum total membrane area required as the sideflow of column C2 increases (i.e., an increase in sideflow of column C2 increases the bottom flow rate of column C1, thus decreasing the total membrane area required).Other than that, the bottom flow rate of column C1, sideflow of column C2, and distillate flow rate of column C2 are correlated to a certain degree, and this is because they are directly involved in the overall mass balance, which must be conserved.
It is not surprising that most of the variables in the H-DWC structure are correlated, as the H-DWC itself is a highly integrated design, and it is therefore expected that there are many interactions between the optimisation variables.It can also be seen that there is not a definite value for any of the variables for the simulations to be off-spec or on-spec.However, the distribution of the on-spec simulations for the bottom flow rate of column C1 is skewed towards smaller values.In comparison, the sideflow and distillate flow rate of column C2 are skewed towards larger values, meaning that there is a higher chance to obtain an on-spec simulation at a small C1 bottom flow rate and large C2 sideflow and distillate flow rate compared to large C1 bottom flow rate and small C2 sideflow and distillate flow rate.This finding makes sense as for a small bottom flow rate of C1, through mass balances, a higher vapour flow from C1 to C2 is expected and a larger sideflow of C2 is required to establish the vapour-liquid equilibrium.Also, a small bottom flow rate of C1 means more heavy products need to leave the system from the membrane system leading to a larger distillate flow rate of C2.As for the total membrane area, the on-spec simulations are normally distributed around 800 m 2 .Finally, there is a minimum required total membrane area of 500 m 2 to meet the specifications.

Overall comparison
In previous subsections, the various hybrid structures were considered and discussed individually.Next, their performance will be compared based on their respective Pareto fronts and modifications.A few fitting equations, such as second to fourth order polynomial, exponential, and logarithmic equations are tested.It was found that the second order polynomial fits all the Pareto fronts with a generally good R 2 value (R 2 > 0.85) and the fittings also describe the shape of the scatter plot well, thus the second order polynomial is chosen to describe the relationship between OPEX and CAPEX for the four different structures.Although higher order polynomials may show a higher R 2 value in some of the cases, higher order polynomials are not chosen as they show oscillations in the fitting which is not what the scatter plot shows.Fig. 13a shows that the OPEX is negatively correlated with the CAPEX for all structures.In other words, as the CAPEX goes up, then as expected the OPEX goes down, although to a varying degree depending on the structure.For the D-P structure, the optimal designs (all simulations on the Pareto front are treated equally) have close CAPEX and OPEX ranges.The other structures show the opposite behaviour, where the optimal designs have a larger range of CAPEX and OPEX, meaning that to reduce either OPEX or CAPEX slightly, a fairly large increase is required in the other cost.The D-P-D and H-DWC structures have similar Pareto fronts, which makes sense as they have similar operating principles.The D-P-D and H-DWC structures can yield a design with relatively low OPEX and CAPEX, while the P-D structure is the least economically attractive.This finding is also reflected in the single objective optimisation with total annualised cost (TAC) as the objective function in Section 3.

Impact of plant life
Since TAC is defined as the summation of OPEX and annualised CAPEX, the Pareto front of each structure can be transferred into a TAC plot with the x-axis as the plant life and the y-axis as the TAC.For each plant life, the Pareto front can be re-generated with CAPEX divided by the plant life while the OPEX remains fixed.It should be noted that the range of CAPEX considered in each structure depends on their respective minimum and maximum CAPEX as indicated in Fig. 13a (e.g., the range of CAPEX considered for P-D structure is from about $ 9.5M to $ 12M).Then, the minimum TAC for the specific plant life can be found by summing the x and y values for every point and choosing the minimum value to represent the minimum TAC.
The TAC plot shown in Fig. 13b indicates that the D-P-D and H-DWC structures always have similar TAC, consistently lower than the P-D and D-P structures.The D-P structure initially (plant life less or equal than two years) has similar TAC to that of the P-D structure, but both are significantly more expensive than the other two designs (D-P-D and H-DWC).With increasing plant life, the contribution of CAPEX to TAC is reduced, which leads to a close TAC for all structures, particularly for the D-P, D-P-D and H-DWC structures.

Discussion
To sum up, considering the findings from both single-objective and multi-objective optimisation studies, the D-P-D and H-DWC structures, whose membrane systems are used only to help cross the composition of the azeotropic mixture to the other side of the azeotropic point, is economically superior to the D-P and P-D structures.This is because using the membrane to only help cross the azeotropic point can significantly reduce the need for the membrane, which is beneficial unless the membrane cost is very low.From the multi-objective optimisation study, it can be deduced that the membrane system is the main contributor to the CAPEX in all the structures, no matter the design.
In terms of energy consumption, the D-P-D and H-DWC structures do not show a saving in energy consumption when compared to the D-P structure.In the multi-objective optimisations, the Pareto fronts (see Fig. 13a) of the D-P-D and H-DWC structures have a similar range of operating costs, both higher than that of the D-P structure, which indicates that the saving in TACs for the D-P-D and H-DWC structures are mainly from reduced capital expenditure.This finding is also reflected in Fig. 13b, where the differences in TAC between the D-P, D-P-D, and H-DWC structures become smaller with increasing plant life.In terms of TAC, the D-P-D and H-DWC structures are preferred as they have similar TACs which are always lower than the other two structures, even with increasing plant life.Comparing the D-P-D and H-DWC structures, the H-DWC structure may be preferred as it only requires one column and one condenser.Therefore, space saving may lead to further TAC saving as space is not considered in the TAC in this work.However, the H-DWC usually requires a taller distillation column, which may be a potential limitation.It should be noted that the above suggestions are only made based on the cost, and factors such as controllability or safety should also be considered for a more comprehensive comparison.
It should be noted that the D-P-D and H-DWC structures considered in this work are limited to specific azeotropic mixtures where on both side of the azeotropic point, distillation could be used for separation and for minimum boiling azeotropes.For a maximum boiling azeotropic mixture, the dividing wall in the H-DWC structure would be extended from the middle to the top, and other designs will be significantly affected as well.
The optimal designs are also sensitive to the feed composition (e.g., for a high molar composition (e.g., 0.9 mol mol −1 ) of ethyl acetate, a single membrane system may be preferred) and to membrane properties (e.g., which component will permeate the most).It is therefore difficult to make a general conclusion for different cases.However, the case studies show that the D-P-D and H-DWC structures can in fact reduce the needs of required membrane area significantly, at the expense of higher energy requirement in the distillation system, when compared to the D-P and P-D structures.It is advised to always design, optimise and compare different hybrid distillation-pervaporation processes before making any decisions.

Conclusion
This work compared three commonly seen hybrid distillation-pervaporation structures for the separation of azeotropic systems: distillation followed by pervaporation (D-P), pervaporation followed by distillation (P-D), and distillation followed by pervaporation and then by distillation (D-P-D).This study also considers a hybrid dividing wall column (H-DWC), which is integrated from the D-P-D structure with the dividing wall extending to the bottom of the column, i.e. combining two columns into one shell.
A single-objective optimisation to minimise the total annualised cost (TAC), and a multi-objective optimisation to minimise both capital cost (CAPEX) and operating cost (OPEX), were performed considering a binary minimum boiling azeotropic mixture as the case study.The single-objective optimisation results showed that the D-P-D structure has the least TAC followed closely by H-DWC (2% higher), then D-P (12% higher) and finally P-D structure (45% higher).The lower TAC for the D-P-D and H-DWC structures is mainly due to the smaller membrane system required when compared to the other two structures, as these structures use the pervaporation unit just to cross the azeotropic point and not for product purification.
The multi-objective optimisation results indicated that, for this case study, the membrane system always contributes the most to the CAPEX, regardless of the structure.Furthermore, the relationship between TAC and plant life indicated that the D-P-D and H-DWC structures consistently have similar TAC, and the difference between these two Instead of modelling the membrane as a distributed model, a lumped model where a membrane module is divided into smaller membrane fragments is used to simulate the mass and energy distribution across the membrane.This lumped model approach is widely reported in the literature and is reported to be sufficient for the modelling of a membrane model (Luyben, 2009;Li et al., 2019;Meng et al., 2020;Wu et al., 2020).Through validation with experimental (Tsuyumoto et al., 1997) and simulation (Marriott, 2001) results, it was found that dividing a membrane module into nine membrane fragments, i.e., N frag.= 9, is sufficient to describe the membrane (see model validation results in Table A1).
The membrane used in this work are taken from Marriott and Sorensen (2003a) and Tsuyumoto et al. (1997).The details of the membrane are shown in Table A2.It should be noted that in Tsuyumoto et al. (1997), the membrane area of 6 m 2 is the effective membrane area at about 75% efficiency, however, Marriott and Sorensen (2003a) took the 6 m 2 as the 100% effective membrane area (which is also used in this work).Therefore, instead of directly taking the fibre radius from Tsuyumoto et al. (1997), Marriott and Sorensen (2003a) (and this work) recalculated the fibre radius from 6 m 2 and 3800 fibres: (A1) In the following equations, the subscripts i denotes the component i (or specifically w for water and e for ethanol), feed denotes the feed side, ret denotes the retentate side, perm denotes the permeate side, and mem denotes the membrane layer, and the term N comp denotes the number of components.

A.1.Membrane fluxes
One of the most important equations in a membrane model is the flux equation.Each membrane has its own flux equation, and this work reports only the flux equation used in the case study for the separation of an ethanol/water mixture, which is obtained from Tsuyumoto et al. (1997) for a polyacrylonitrile ultrafiltration hollow-fibre membrane PAN-B5.The necessary information for the membrane can be found in Table A2.A solution-diffusion approach is used to formulate the equation for the flux of water through the membrane, J W (g m −2 h −1 ), which is given by: = + where the term D w0 is the diffusion coefficient of water at infinite dilution, K cw is the sorption coefficient, k dw is a numerical constant, δ is the membrane thickness, γ is the activity coefficient, x is the molar composition, P is the pressure, and P sat is the saturated vapour pressure.γ and P sat are obtained from Multiflash (KBC Advanced Technologies, 2015).
The equation for the flux of ethanol through the membrane, J e (g m −2 h −1 ), is described with a simple equation: where the subscripts feed denotes the feed side, perm denotes the permeate side, and e denotes ethanol, and the terms L p is the permeability constant which is membrane-dependent, ω is the mass fraction, and P is the pressure.Tsuyumoto et al. (1997) claimed that an average value, L p = 5 × 10 −3 g m −2 h −1 torr −1 , could be used for the membrane used in this work.
The terms (D w0 K cw ) and D K k ( 2) in the water flux equation (Eq.(A2)) are given by the equations below (Tsuyumoto et al., 1997) where T feed (K) is the temperature at the feed side.

A.1.1. Modifications to membrane fluxes
The mixture considered in the literature is ethanol-water mixture, but the mixture considered in this work is the ethyl acetate-ethanol mixture.However, due to lack of publicly available data on membranes for separation of ethyl acetateethanol mixtures, the membrane flux equations described in Tsuyumoto et al. (1997) for ethanol-water separation are modified for the separation of ethyl acetate-ethanol mixture according to the boiling points of the components as follows: • The ethanol flux through membrane in this work uses the water flux through membrane in Tsuyumoto et al. (1997) (Eq.( A2)).
• The ethyl acetate flux through membrane in this work uses the ethanol flux through membrane in Tsuyumoto et al. (1997) (Eq.( A3)).
This assumption is reasonable because the aim of this work is to study the feasibility and performances of different hybrid structures relative to each other, and not to generate accurate results for the specific membrane using this mixture, and assumption also used by Barakat (2006).

A.2.Molar balance equations
The component molar balances at the retentate and permeate sides are given by: where x is the molar composition, F is the molar flow rate, A fibre is the surface area of the fibre, J is the flux, and M is the molar holdup.The surface area of the fibre is calculated by: where A mem is the membrane area.Then, under steady-state conditions: perm perm i fibre i , (A10)

A.3.Pressure drop equations
In this work, the permeate side pressure is maintained at 400 Pa (Tsuyumoto et al., 1997), and it is assumed that the pressure drop across the membrane at the permeate side is negligible (Assumption 4).For the retentate side, which is at the fibre side (Assumption 2), one of the most commonly used equations to calculate the pressure drop is the Hagen-Poiseuille equation for laminar flow (Pan, 1986;Lipski and Coˇté, 1990;Marriott, 2001;Kookos, 2002;Katoh et al., 2011;Li et al., 2019): where the subscript fibre denotes the fibre side, the terms ΔP is the pressure change/drop, μ is the dynamic viscosity of the liquid, L mem is the length of the membrane module, V is the vo- lumetric flow rate, and r is the radius.The validity of the Hagen-Poiseuille equation can be examined with the Reynolds number, Re, where if Re < 2100, then the flow in the fibre can be considered to be laminar (Lipski and Coˇté, 1990).

A.4.Energy balance equations
There are many reports in the literature of energy balance equations (Marriott, 2001;Hafrat et al., 2016;Lee et al., 2016;Meng et al., 2020;Babaie and Nasr Esfahany, 2020), and these questions are also used in this work, where: where h is the specific enthalpy, F is the molar flow rate, and M is the molar holdup.Under steady-state conditions, Eq. (A12) can be simplified to: The specific enthalpy is obtained from Multiflash (KBC Advanced Technologies, 2015), where it is a function of the temperature, pressure, and composition of the retentate side, h = f(T, P, x).

B. Costing equations
This section presents the equations used to calculate the costs of the units (other than for the distillation column, which can be found in Duanmu et al. (2022b)) used in this work.The equations for the membrane are taken from González and Ortiz (2002), while the equations for the other units are obtained from Sinnott and Towler (2020).The parameters used in calculating capital, operating, and total annualised costs are shown in Table B1.
In general, the capital cost (CAPEX) and operating cost (OPEX) are used to calculate the total annualised cost (TAC) with the following equation: where the CAPEX and OPEX take into account all the units and utilities, respectively.The plant life and annual operating hours are assumed to be eight years and 8400 h y −1 , respectively, in the case studies.

B.1.Membrane
The costs of the pervaporation membrane are taken from González and Ortiz (2002), which are based on the prices in 2007.Therefore, appropriate scaling using the Chemical Engineering Plant Cost Index (CEPCI) should be applied.The capital cost (CAPEX) of the membrane can be calculated by (González and Ortiz, 2002): where f Lang is the Lang factor (see Table B1), Price mem = $ 1063 m −2 is the price per area of the membrane, and A tot,mem is the total membrane area required.The membrane replacement cost, which is calculated as part of the operating cost (OPEX), can be calculated from (González and Ortiz, 2002): where Price repl. is the membrane replacement cost per area per year in 2002 (Price repl.= $ 200 m −2 y −1 ) taken from González and Ortiz (2002) and t mem is the membrane lifetime which is assumed to be two years (González and Ortiz, 2002).

B.2.Heaters and coolers
The heaters and coolers (including the membrane network heaters) are considered as U-tube shell and tube heat exchangers.The CAPEX of the heat exchanger can then be calculated (Sinnott and Towler, 2020): ) where f Lang and f m are the Lang factor and material factor, respectively (see Table B1), and A HEX (m 2 ) is the heat exchanger require calculated by (Luyben, 2013): where Q HEX is the heating/cooling duty, and U and ΔT are the heat transfer coefficient and typical temperature difference, respectively (see Table B1).
The operating cost for the heat exchanger (heaters or coolers) can be calculated by:

(B6)
where Price util. is the price of the (heating or cooling) utility used, and the type of the utility is decided automatically by Table B1 -Values and references of the parameters used for the calculation of capital, operating, and total annualised costs.
the optimiser depending on the outlet temperatures of the heat exchangers.

B.3.Pumps
The cost equation for the pump is taken from Sinnott and Towler (2020): where f m is the material factor (see Table B1) and V L s ( ) 1 is the inlet volumetric flow rate.
The operating cost of the pump can be calculated by: where Price elec. is the price of the electricity (see Table B1) and P pump (kW) is the power required by the pump.

Fig. 1 -
Fig. 1 -Schematic of the membrane network.The membrane stages (MS) are depicted with the larger boxes, while the membrane modules (MM) are the smaller boxes with a diagonal line enclosed within a membrane stage, and HEX(n) is the feed heater for membrane stage n.There are N ms membrane stages in a membrane network, and there are N mm,n membrane modules in the n th membrane stage.

Fig. 4 -
Fig. 4 -Donut chart showing the portions of annualised capital cost (CAPEX) and operating cost (OEPX), together with the elements that make up the CAPEX and OPEX, in the total annualised cost.

Fig. 5 -
Fig. 5 -Distillation-Pervaporation (D-P) structure: (a,b) distribution of simulations of capital cost (CAPEX) vs. operating cost (OPEX), (c-f) correlation between total membrane area and another main optimisation variable in terms of CAPEX and OPEX for all feasible (successful simulation) points.(Note that on-spec means that all constraints are met).

Fig. 7 -
Fig. 7 -Pervaporation-Distillation (P-D) structure: (a,b) distribution of simulations of capital cost (CAPEX) vs. operating cost (OPEX), (c-f) correlation between total membrane area and another main optimisation variable in terms of CAPEX and OPEX for all feasible (successful simulation) points.(Note that on-spec means that all constraints are met).

Fig. 9 -
Fig. 9 -Distillation-Pervaporation-Distillation (D-P-D) structure: (a,b) distribution of simulations of capital cost (CAPEX) vs. operating cost (OPEX), (c-f) correlation between total membrane area and another main optimisation variable in terms of CAPEX and OPEX for all feasible (successful simulation) points.(Note that on-spec means that all constraints are met).

Fig. 11 -
Fig. 11 -Hybrid dividing wall column (H-DWC) structure: (a,b) distribution of simulations of capital cost (CAPEX) vs. operating cost (OPEX), (c-f) correlation between total membrane area and another main optimisation variable in terms of CAPEX and OPEX for all feasible (successful simulation) points.(Note that on-spec means that all constraints are met).

Fig. 13 -
Fig. 13 -Graphs showing (a) the final Pareto fronts obtained by re-ranking the accumulated Pareto fronts for each structure, and (b) the effect of plant life on the total annualised cost (calculated from the fitted equations).