Distillation operations consume a significant amount of energy. The U.S. Energy Information Administration predicts an increase in energy consumption by U.S. industries from 26 quadrillion British thermal units (Btu) to 36 quadrillion Btu between 2019 and 2050 ( U. S. Energy Information Administration, 2020 ). The bulk chemical and refining industries are the most energy intensive, representing respectively 29% and 18% of total U.S. industrial energy consumption ( U. S. Energy Information Administration, 2020 ). In both sectors, separation processes represent approximately 60% of the energy usage, and about 95% of these separation processes are distillation operations ( U.S. Department of Energy, 2005 ). This significant energy demand has made distillation the focus of many studies aiming to improve its energy efficiency. Among these studies, the division of internal column streams has been shown to be an efficient mechanism for energy improvement. This technique is known as distillation with parallel streams, and it includes the parastillation, with division of the vapor phase; metastillation, with liquid division, and the best-known dividing-wall column (DWC), with liquid and vapor divisions.
The DWC integrates two conventional columns of the Petlyuk system into a unique column shell. This configuration, in comparison to the conventional one, reduces the number of condensers and reboilers, and thereby decreases the capital and operational costs by up to 30% and 40%, respectively ( Kiss and Suszwalak, 2012 ). Reports point to 100 DWCs operating industrially ( Luyben, 2013 ). However, DWC configurations have been used for specific processes, where two conventional columns are replaced by one DWC, mainly in the context of ternary mixture separations. Single distillation columns, as the ones used in binary separations, can be replaced by parastillation and metastillation columns, with energy and capital savings. These are less-known techniques that consider the vapor (parastillation) or the liquid (metastillation) internal flow divisions, into two or more parallel streams.
Two-liquid divisions, in metastillation, can decrease the stage area up to 30%, in comparison to a corresponding conventional distillation column ( Mizsey et al., 1993 ). Additionally, metastillation columns may present Murphree efficiencies greater than those reported for conventional columns ( Gouvêa, 1999 ). On the other hand, the two-vapor divisions, in parastillation columns, lead to reductions in column height of up to 50%, when compared with conventional distillation columns ( Biasi, 2016 ;Gouvêa et al., 20 0 0 ;Heucke, 1987 ;Meirelles et al., 2017 ;Mizsey et al., 1993 ). Furthermore, replacing a conventional distillation column by a parastil-lation column can lead to a decrease in energy consumption of about 30% ( Biasi, 2016 ;Canfield and Jenkins, 1986 ;Gouvêa, 1999 ;Meszaros and Fonyo, 1990 ;Moraes, 2006 ).
The construction of lab-scale parastillation columns were reported by Belincanta et al. (20 06 ,20 05 ). The authors investigated the separation of the ethanol-water mixture and did not report operational issues related to the phase division. Even with the important advantages presented by para-and metastillation, there is only one report about the industrial installation of three parastillation columns in England, ( Canfield and Jenkins, 1986 ) and none about metastillation. The report by Canfield and Jenkins (1986) mentions that one of the parastillation columns presented a reduction in the reflux ratio of 41%, over conventional distillation. However, the authors do not provide further information, stating proprietary issues. They concluded that industrial parastillation columns were completely satisfactory, improving separation and reducing energy consumption over conventional distillation ( Canfield and Jenkins, 1986 ). There is no further report about the industrial use of para-and metastillation columns. This can be attributed to the lack of scientific and technological information about these columns in comparison with conventional distillation equipment, and the specific difficulties of simulating paraand metastillation processes. Most of the few academic works on para-and metastillation columns have investigated simplified processes, but did not consider potential industrial cases. Except for the works of Meszaros and Fonyo (1990) and Mizsey et al. (1993) , that consider hydrocarbon multicomponent mixtures, all other works consider simple binary separations. Furthermore, there is no work concerning the control structure for para-and metastillation columns. Some information about the hydrodynamics of parastillation columns is provided in Belincanta et al. (20 06 ,20 05 ). Given the lack of information on these new types of columns, greater industrial acceptance requires more robust research on the subject, since current distillation processes require equipment with optimized performance that is designed on the basis of wellestablished methodologies.
This paper presents a generic tool for simulation of any type of distillation equipment, including conventional columns, paraand metastillation columns and the three different DWC configurations. Up to now, Aspen Plus® and other commercial simulators do not cover para-and metastillation columns. Although studies have shown that it is possible to simulate the DWC equipment using commercial simulators ( Kiss, 2013 ), this alternative is not straightforward to execute. A DWC is represented, in the simulators, by the integration of more than one column. However, the resulting system of equations is not solved as a fully coupled system using a Newton-type method, but instead by a nonlinear Gauss-Seideltype iteration which alternatingly solves the equations of one column given values for the other. This solution approach has poorer convergence properties than those based on the fully coupled system.
In this work, we present a new unified model that allows the simulation of distillation columns with any integer number of phase divisions. Our model is based on the so-called MESH equations ( Seader et al., 2011a ) and avoids the simplifying assumptions made in short-cut methods for columns with divided streams in the vapor or liquid phases ( Heucke, 1987 ;Meirelles et al., 2017 ). The proposed model uses a unique set of MESH equations for all columns, including parastillation, metastillation, dividing-wall column (DWC), and conventional distillation. By small adjustments of the indices of the MESH equations, it is also possible to adapt the model for simulation of other DWC configurations, as the bottom and top-DWC. The model is simply adjusted, to a specific column, by specifying few parameters and, therefore, can be used to analyze capital or operational cost savings that may be obtained by using columns with phase divisions, in comparison to conventional distillation columns. Newton's method is used to solve the model equations and does not suffer from the convergence problems observed with previous approaches. Previous simulations ( Biasi, 2016 ;Canfield, 1984 ;Canfield and Jenkins, 1986 ;Gouvêa, 1999 ;Meszaros and Fonyo, 1990 ;Mizsey et al., 1993 ) required a tailored algorithm for each type of column, which is time consuming and therefore costly. Using a unique set of equations, our approach reduces the implementation time and facilitates the comparison of different column configurations. This procedure allows systematic analyses and comparisons of such different configurations and, therefore, supports the design of new, energy efficient and capital cost reducing columns.
We generalize the mathematical model of the MESH equations by introducing two new variables: θ and β. These two variables represent the number of phase divisions in the liquid ( θ ) and vapor ( β) phases. It enables the simulation of columns with multiple phase divisions, using a unique set of equations. The suggested mathematical model was implemented and validated. Conventional distillation and DWC were validated comparing our results with those from Aspen Plus® commercial simulator and from Luyben (2013) . The validation for metastillation and parastillation was done by comparing with results from the literature ( Gouvêa, 1999 ;Gouvêa et al., 20 0 0 ). We also provide some details about the convergence and some tips used to decrease the computational time required by the simulation.
A case study illustrates the use of para-and metastillation columns applied to the production of hydrous bioethanol ( Batista et al., 2012 ), which is a case of industrial interest. Results demonstrate that para-and metastillation techniques lead to savings in equipment and/or operational costs. It is possible to reduce operational costs by about 19% through the replacement of conventional distillation by parastillation. Additionally, the height of parastillation columns may be reduced by 33%, when compared to the height of conventional columns. These reductions can decrease the total annual cost of distillation processes in 15%. These results, associated with the unified model for simulating columns with multiple phase divisions, may increase the industrial interest in para-and metastillation columns. In addition, the unified model simplifies the simulation of the DWC equipment and allows a faster comparison of the performance of different column configurations.
The methodologies for simulating distillation columns can be classified as short-cut methods or rigorous simulation procedures. The short-cut methodologies are based on simplified models, for example, models that assume constant molar overflow along the column and/or constant relative volatility. One well known method is the MacCabe-Thiele approach. Meirelles et al. (2017) presented a general approach for the calculation of parastillation and metastillation columns based on the adaptation of the original McCabe-Thiele methodology. The short-cut approaches are an important tool for preliminary-design studies ( Seader et al., 2011a ). They can be useful for a fast but rough calculation of the benefits of the alternative column configurations in comparison to conventional distillation processes. However, these methods are based on several simplifications, they do not adequately consider the simultaneous heat and mass transfer occurring inside the distillation apparatus, and are not suitable for simulation of multicomponent processes. In these cases, a rigorous method is necessary, such as the approach developed by Naphtali and Sandholm (1971) .
The Naphtali and Sandholm (1971) methodology considers multicomponent mixtures, non-ideal solutions and Murphree efficiencies, in multiple feed and multiple side-stream columns. This methodology is based on the MESH equations, which consist of mass balances (M); equilibrium relationships (E); sums of the molar fractions of the components in the liquid and vapor phases (S); and enthalpy balances (E). The set of nonlinear MESH equations can be solved by any iterative method suitable for non-linear systems. Newton's method was suggested by Naphtali and Sandholm (1971) and implemented by Fredenslund et al. (1977) . The Naphtali and Sandholm (1971) methodology was originally proposed for the simulation of conventional distillation columns, but it was subsequently extended to parastillation by Meszaros and Fonyo (1990) and by Mizsey et al. (1993) and to para-and metastillation by Gouvêa (1999) and by Biasi (2016) . All these works have implemented separate algorithms for each type of column, and they have been limited to only two partitions of the selected phase.
We formulate a new set of MESH equations suitable for the simulation of distillation columns with any integer number of divisions of the vapor or/and liquid internal streams. Such a model is suitable for the simulation of conventional distillation, parastillation, metastillation and dividing-wall (DWC) columns. The main novelty of our work is the possibility of simulating a variety of internal configurations of distillation columns using a single model, and the numerical implementation of this model that can be adjusted to different columns simply by setting input parameters.
Our model for columns with multiple phase divisions is an extension of the method proposed by Naphtali and Sandholm (1971) , which was originally developed for conventional distillation columns. The Naphtali and Sandholm (1971) method combines the S equations (sum of the molar fractions of components in liquid or vapor phases equal to one) with other MESH equations to eliminate 2 N variables. It reduces the system to a total of ( 2 C + 1 ) • N equations, where N and C are the total number of stages and components, respectively. In this approach the total flow rates and component molar fractions are combined to generate equations based on component flow rate variables ( Seader et al., 2011a ). In brief, the variables to be solved for are: N • C component vapor flow rates, N • C component liquid flow rates and N temperatures. These equations can be solved by any nonlinear solver.
The main novelty of our model is the incorporation of the number of phase divisions, into the MESH equations. This allows the simulation of different parallel stream columns, without further model modification. Previous works have developed a new set of equations for each process with parallel streams, without considering columns with more than two phase divisions ( Biasi, 2016 ;Gouvêa et al., 20 0 0 ;Meszaros and Fonyo, 1990 ;Mizsey et al., 1993 ). Traditionally, the MESH equations for conventional distillation considers that a generic stage n receives the vapor and liquid flows from stages n -1 and n + 1 , respectively, when counting the stages from bottom to top ( Naphtali and Sandholm, 1971 ). Our extension considers that a stage n receives liquid from stage n + θ and vapor from n -β, where θ and β are the numbers of liquid and vapor phase divisions, respectively.
The assumption of θ liquid phases and β vapor phases allows the simulation of different configurations. For exam ple, a conventional distillation column with one liquid and one vapor stream can be simulated considering β = θ = 1 . This specification leads to the original model proposed by Naphtali and Sandholm (1971) , in which a conventional generic stage n receives vapor from stage n -1 , and liquid from stage n + 1 . In parastillation, the vapor flow is divided into two or more ( β ≥ 2) ascending parallel streams, which contact a single downward stream of liquid ( θ = 1 ). This means that a generic stage n receives liquid from n + 1 , but the inlet vapor is from n -β. Analogously, in the case of metastillation, the division occurs in the liquid stream. Two or more downward liquid flows ( θ ≥ 2) contact a single upward vapor stream ( β = 1 ). For a DWC with two column divisions in the entire apparatus, there are two liquid flows ( θ = 2 ) and two vapor flows ( β = 2 ), dividing the column into two parts by a wall. In a DWC the phases from one side of the column do not contact the phases from the other side of the wall. Prior to the simulations, one must enter the internal number of phase divisions ( θ and β). This procedure will automatically select the column configuration.
In our model, the inlet liquid comes from θ stages above, and the vapor from β stages below, counting the stages from bottom to top. For all configurations, the columns have a total number of N stages, with stage n = 1 being the reboiler and stage n = N being the condenser. Fig. 1 illustrates the stages numbering scheme for four different columns configurations: conventional distillation, metastillation (with θ = 2 ), parastillation (with β = 2 ), and DWCs with whole division of the column. For all the cases, a sequential numbering scheme has been applied, counting from the bottom to the top. This sequential scheme is often used to represent conventional distillation columns, as seen in Seader et al. (2011a) . Some authors also consider counting from top to bottom, as often used in Aspen Plus®. The counting orientation do not affect the solutions and the algorithm convergence.
For conventional distillation ( Fig. 1 ), the reboiler is the first stage ( n = 1 ); the second stage ( n = 2 ) is the first tray of the column; the second tray is the third stage ( n = 3 ), and so on until it reaches the top of the column, where the condenser is the last stage ( n = N). Metastillation columns have the same numbering scheme as the conventional distillation, in which the stages are numbered consecutively from bottom to top. That is, stage ( n -1 ) is located just below the stage n , while stage ( n + 1 ) is immediately above stage n . This counting scheme is valid for metastillation columns with any number of phase divisions, although Fig. 1 only illustrates a column with two liquid divisions.
For parastillation with two vapor divisions and for DWC equipment (with whole division of the column), one needs to consider stages from different sides of the column as consecutive stages. It means that the stages from one side have even numbers, while the stages on the other side have odd numbers. With this procedure, the first tray on the right side of these columns is assigned stage n = 2 , while the second tray of this same column side is stage n = 4 . On the left side, the first tray is assigned stage n = 3 and the tray immediately above this one is stage n = 5 . Overall, the right side contains stages n = 2 , 4 , . . . , n -2 , n, n + 2 , . . . and the left side contains stages n = 3 , 5 , . . . , n -1 , n + 1 , n + 3 , . . . . Parastillation columns with more than two divisions ( β ≥ 2) are not illustrated in Fig. 1 . However, the counting scheme for these columns is analogous to the column with two divisions ( β = 2 ). For example, for three divisions ( β = 3 ), the first column partition (side) contains stages n = 2 , 5 , 8 . . . , n, n + 3 , . . . . The second partition contains stages n = 3 , 6 , 9 , . . . , n + 1 , n + 4 , . . . . The final third column partition contains stages n = 4 , 7 , 10 , . . . , n + 2 , n + 5 , . . . .
The assumption of θ liquid and β vapor phases and the counting scheme illustrated in Fig. 1 permits us to represent a generic stage n . This stage represents any arbitrary configuration receiving liquid from n + θ and vapor from n -β, as shown in Fig. 2 . As in the MESH equations, the total liquid and vapor molar flows are denoted by L n and V n , respectively. The component molar flows are denoted by l n, i (for liquid phase) and v n, i (for vapor phase), where n is the stage and i the component. The enthalpies of the liquid, vapor and feed flows are H L n , H V n and H L/V f n , respectively. The molar flow of each component i of the external feed is represented by f n, i , and the liquid and vapor side-streams by S L n and S V n , respectively. The separation efficiency at stage n for component i is represented by the Murphree efficiency E MV , n, i . The temperature at stage n is denoted by T n . The MESH equations are still dependent on the reflux ratio r n and the phase-equilibrium ratio or k -value for each component i, k n, i . The k-values usually depend on stage pressure, temperature and liquid/vapor molar fractions.
The MESH equations are summarized in Table 1 . In addition to these equations, we need to add equations that express the total liquid and vapor flows (
v n, j ), the phaseequilibrium ratios ( k n, i ) and the enthalpies ( H L n and H V n ) , in terms of the component molar flows ( l n, i and v n, i ) and temperatures ( T n ). These equations are standard and similar to those found in reference books ( Kister, 1992 ;Seader et al., 2011b ;Treybal, 1981 ) and Discrepancies functions for distillation columns with parallel streams with any number of liquid ( θ) or vapor ( β) phase division, valid for all stages (n from 1 to N) and components (i from 1 to C).
Component mass balance ( j = 1), for components used in commercial simulators as Aspen Plus® ( Aspentech, 2001 ).
For each case study one must select the best equilibrium relations and enthalpy equations to represent the simulated mixture. For more details about the phase-equilibrium ratios k n, i , see Seader et al. ( Seader et al., 2011c , pp. 35-84, Eq. 2.19) and Kister (1992 , pp. 3-11, Eq. 1.1, 1.11). For details about enthalpy, see Treybal ( Treybal, 1981 , pp. 357-363) and Aspentech ( Aspentech, 2001 , Eq. 34, 39). We solve the set of MESH equations for the molar flows l n, i and v n, i and temperatures T n . As in Fredenslund et al. (1977) and Seader et al. (2011a) , for a given set of variables, the residuals of the equations in Table 1 are called discrepancies or errors and are denoted by Dj ( n, i ), where j represents the type of function, n the stage of the column, and i refers to the component of the mixture of C components.
In the MESH equations of Table 1 the total number of stages, number of phase divisions, the feed conditions (temperature, pressure, component flow and position), the side streams conditions (flow and position), the Murphree efficiencies, the stage pressure, the reflux ratio, and the distillate flow are specified parameters. Specifying the reflux ratio and the distillate flow, it is possible to avoid the specification of the heat flow added to the reboiler ( QR 1 ) and removed from the condenser ( QC N ). This is done by replacing the enthalpy balances from stages n = 1 and n = N by global mass balances. In Table 1 , these mass balances represent columns with partial condenser, with distillate flow as vapor. For a total condenser, discrepancy equations D 3(1) and D 3(N) must be replaced by the correspondent global mass balances. With the solution of the MESH equations, the heat flows can be calculated as indicated by Eqs. ( 1) and (2) .
The equations shown in Table 1 are suitable for simulation of DWCs by dividing the whole column into two parts. This means that there is no contact between the liquid or vapor phase from different sides of the column, except in the condenser and in the reboiler. This procedure does not represent the Petlyuk configuration ( Luyben, 2013 ), in which the phase division occurs just in the middle of the column and the bottom and top stages are composed by conventional distillation trays ( Fig. 3 ). Furthermore, the conventional DWC can be adapted to bottom-and top-DWCs ( Fig. 3 ). The bottom-DWC considers the division of the column at the bottom into two parts, while the top is composed by conventional distillation trays. This type of column has two reboilers, with each one receiving the liquid flow from one side of the column. There is no mixture of the outlet liquids from each side of the column. Analogously, in the case of top-DWC, the phase division occurs only at the top of the column. This column has two condensers that collect the vapor flows from the two top column divisions.
For simulations of DWC columns containing conventional distillation stages on top and/or on its base, small changes in the indexes of the MESH equations must be considered. Additionally, the extra condensers and reboilers, present in top-and bottom-DWCs, respectively, require extra equations. For these DWCs, we adopted a different counting system for the stages, than the one represented in Fig. 1 . This system considers that different sides of the column may contain different numbers of stages. Fig. 3 shows the new counting scheme. Additionally, the new set of equations considers the internal phase division in proportions different than 50%. The left side of the column receives rv % of the vapor phase and rl % For bottom-DWC, the number of stages in the bottom section is equal to zero and the first stage of the column is na with index 2. For top-DWC the number of stages at the top is null, so the last stage of the column is ( nc -1 ) = ( N -1 ) . For some numerical example, let us consider a traditional DWC with 7 trays in the bottom section, 10 trays in the top section, 15 trays in the right side, and 20 trays in the left side. The parameters, of Eqs. ( 3
and nc = ( 17 + 20 ) = 37 . Note, that for this column the nd parameter has no physical meaning once the stage 1 represents one of the reboilers ( Fig. 3 ). This bottom-DWC has a total of 48 stages, considering 45 trays plus two reboilers and one condenser.
The DWC columns can be simulated considering the main MESH equations, represented by Eqs. ( 7) -( 9) . These equations are valid for all components (1 < i < C ) and for all stages ( 0 < n < N + 1 ). Eqs. ( 7) -( 9) depends of the indexes for the liquid ( Index liquid ) and vapor ( Index vapor ) phases, which represents the inlet flows, and of the multiplication factors ( MF liquid for the liquid phase and MF vapor for the vapor phase), which represent the percentage of inlet flow that each stage receives. The indexes and multiplication factors of these equations are represented in Table 2 . In Eq. ( 8) , QR n and QC n represent the heat added to the reboilers and removed from the condensers, respectively. The main reboilers are in stages n = 0 and n = 1 , and the main condensers in stages n = N and n
The heat added to the reboilers and removed from the condensers can be specified, or one can specify other properties, as reflux and boil-up ratio and/or bottom and top flows. In this case, the energy equations ( Eq. ( 8) ) must be replaced by global or local mass balances containing the new specified properties.
a N/A = Does not apply, means that there is no equation for this stage.
b Null index, means that ( l ( Inde x liquid ) ,i ) or ( v ( Inde x v apor ) ,i ) are equal to zero. c Equations with two induces means that there are two liquid ( l
Alternative discrepancy equations to the energy balances.
Replacement for stages 0 and 1 Replacement for stages N and N + 1
is based on the equations presented by Seader et al. (2011a) , and contains some replacement suggestions. These equations consider partial condensers, with distillate product as vapor.
The mass balance (M), enthalpy balance (H) and equilibrium relations (E) correspond to a total of ( 2 C + 1 ) • N equations and ( 2 C + 1 ) • N variables ( Tables 1 and2 ). The variables are the molar liquid and vapor flows ( l n, i and v n, i ) and the temperatures ( T n ). The system of nonlinear equations can be solved by any method suitable for nonlinear relations. For conventional distillation, Newton's method was suggested by Naphtali and Sandholm (1971) and described in detail by Seader et al. (2011a) and Fredenslund et al. (1977) . We implemented a slightly modified Newton's method using MatLab®. We opted for a particular way of implementing Newton's method to have a better control throughout the calculation process. Specifically, we organize iterates to ensure that variables are always within physically significant ranges. For example, the calculation of component flows was organized in such a way as to generate always non-negative values, as recommended by Fredenslund et al. (1977) . See Section 3.3 for additional details.
Newton's method uses the linearization of the equations, and each iteration of the method requires the solution of a linear system with the equations given by the Jacobian matrix. The Jacobian matrix is composed of the derivatives of all discrepancy functions ( Tables 1 and2 ) with respect to all variables ( l n, i , v n, i , T n ). We have used finite difference approximations of the Jacobians. For problems with one and two phase divisions, we also used analytic Jacobians, for comparison purposes. For these cases, we compared the numerical results obtained with analytical Jacobians against those with finite difference approximations and both gave essentially the same convergence behavior and computed solutions.
The Jacobian matrix is sparse since each discrepancy function depends only on a few variables. For conventional distillation, the Jacobian matrix has a simple sparsity pattern. Non-zero entries are located at the three neighboring diagonals ( Seader et al., 2011a ), since in conventional distillation stage n depends only on the stages immediately above and below it. In the case of columns with different β's and θ 's the sparsity pattern for the Jacobian is a bit more complicated. The discrepancy functions for stages n from ( β + 2 ) to ( N -( θ + 1 ) ) are dependent only on output variables of stages ( n + θ ), ( n ), and ( n -β). For stages from 2 to ( β + 1 ), the discrepancy functions are also dependent on variables of stage n = 1 , and for stages from ( N -θ ) to ( N -1 ) they are dependent on n = N.
As mentioned before, Newton's method is used to approximately solve the discrepancy equations Dj ( n, i ), which are represented in Table 1 and Eqs. ( 7) -( 9) . Newton's method is terminated when the sum of the squares of the discrepancy functions is small:
where ε is the residual tolerance. Fredenslund et al. (1977) adopted an ε value of 0.1. However, we observed that for ε = 10 -4 there is still a substantial difference between the latest iterates. Therefore, we suggest an ε value of 10 -6 or smaller. Note, however, that to reduce the residual from 0.1 to 10 -6 usually requires only a few iterations, due to the quadratic convergence of Newton's method close to a solution of the MESH equations, when the Jacobian matrix is nonsingular.
The implementation of the MESH equation, as the design of Jacobian matrices, are reported in detail by Seader et al. (2011a, pp. 393-400 ). The authors illustrated a flowchart representing the algorithm of Newton-Raphson method applied to conventional distillation columns ( Seader et al., 2011a , p. 397). The main differences between this flowchart and the procedure here presented is: (i) the replacement of the traditional MESH equations by the ones presented in Tables 1 and2 , and (ii) the introduction of the new parameters β and θ , used to represent the number of phase divisions. The β and θ parameters must be specified in the beginning of the simulations, together with other process information, such as the number of stages, the reflux ratio, the feed conditions, etc.
Different types of distillation processes can be activated automatically through the β and θ parameter assignment. For example, the conventional column can be simulated specifying β = θ = 1 . With this specification, the equations of Table 1 resume to the MESH equations reported to conventional distillation columns ( Naphtali and Sandholm, 1971 ). To select the parastillation option, one must specify θ = 1 and β as any integer number bigger or equal to 2 ( β ≥ 2). For metastillation, one must specify β = 1 and θ ≥ 2. The equations, of Table 1 , also represent a DWC column whole divided into two parts. To simulate this column, the parameters are β = θ = 2 . However, Table 1 does not cover the DWC corresponding to the Petlyuk configuration, and the bottom and top-DWCs processes. To simulate these configurations new discrepancy equations must be implemented, considering the equations presented in Section 2.1.2 . These adaptations introduce more four variables ( na, nb, nc, nd ), that controls the location of the division wall in the column. These variables can be assigned by specifying the number of stages in each column section (as seen in Eqs. ( 3) -( 6)).
In order to validate our proposed model, we compared the results generated by our model with the literature and with results generated using Aspen Plus®. The first comparison considered conventional distillation columns, which can be simulated by the present model assigning the β and θ parameters equal to the unit ( β = θ = 1 ). The initial iterate to start Newton's method were calculated based on the procedure described by Fredenslund et al. (1977) . The simulations were carried out on an Intel Core i7 processor running at 3.6 GHz.
Newton's method applied to 150 equations (30 conventional distillation stages, reflux ratio of 5.04) converged in 7 iterations, considering a methanol (50 mol%) + ethanol (50 mol%) mixture, the Peng-Robinson thermodynamic package, and 100% of Murphree efficiency. It required a run time of approximately 5 s. Applied to 1500 equations (300 conventional distillation stages and reflux ratio of 5.04) Newton's method converged in 15 iterations, and required a run time of 70 s. Note that, the construction of conventional distillation columns with 300 stages is impractical. However, this example was used to verify the method's convergence under extreme cases. It is important to highlight that columns with phase divisions require more stages than conventional ones, without a necessary increase in the column height, justifying this kind of analysis. For example, a parastillation column, with six vapor divisions ( β = 6 ) and 300 stages, has the same height as a conventional column with 53 stages. Conventional columns with more than 40 stages are common in ethanol distilleries plants ( Batista et al., 2013( Batista et al., , 2012 ) ) and in other processes.
In addition to the two-component conventional distillation ( β = θ = 1 ), we also considered a multicomponent distillation for validation purposes. Specifically, the feed was composed of water (518.04 kmol/h), ethanol (14.36 kmol/h), isoamyl alcohol (1.62.10 -2 kmol/h), ethyl acetate (8.73.10 -4 kmol/h), and glycerol (4.74.10 -2 kmol/h). For this example, the NRTL model was adopted to determine the non-ideality of the liquid phase. The vapor phase and the enthalpies were considered ideal. The total number of stages, including the condenser and reboiler, was 36. The reflux ratio was set at 4.5, the distillate flow at 17 kmol/h. The side-stream of 0.3 kmol/h was located three stages above the feed stage. The mass fraction of ethanol in the distillate was 93%.
Considering the multicomponent distillation, Newton's method solved the 396 equations (36 conventional distillation stages and five components) in 115 iterations. The larger number of iterations is due to the isoamyl alcohol, which accumulates in the middle of the column hampering the convergence. Without this component, Newton's method converges in 27 iterations. The results obtained with the conventional distillation were used as the initial guess for the simulation of parastillation and metastillation columns. Keeping constant the number of equations ( 396), Newton's method required 75 and 79 iterations to solve the equations for meta-and parastillation columns, respectively.
The conventional distillation results were compared with results generated by Aspen Plus®. For all components and stages the molar fraction and temperature profiles computed using our new approach agreed with those computed by Aspen Plus®. The condenser and reboiler duties were also in accordance (the difference from the two simulations was within 2%). The correlation coefficient was 0.99 for ethanol, water and glycerol, and 0.98 for isoamyl alcohol and ethyl acetate. The differences may be associated to the fact that even though both methods used the NRTL thermodynamic package, the enthalpy calculations are different for both simulators. This work considers ideal enthalpy, while Aspen Plus® also considers excess enthalpy. Nevertheless, the observed differences were very small, occurring mainly in stages where the concentration of the minor components tended to zero.
There are few works about meta-and parastillation with two divisions in the literature, and no work concerning rigorous distillation with more than two divisions. The methanol + water distillation was simulated in parastillation columns by Gouvêa et al. (20 0 0) , and in metastillation columns by Gouvêa (1999) . The columns had an inlet feed flow of 100 kmol/h, composed by 60 mol% of methanol and 40 mol% of water, and a distillate flow of 40 kmol/h. The reflux ratio was 2.0 for parastillation and conventional distillation and 1.1 for metastillation. The authors used the UNIQUAC model to determine the activity coefficient, in this work we used NRTL. In both works vapor phase and enthalpies were considered ideal. Parastillation columns with two vapor divisions can be simulated considering θ = 1 and β = 2 , in the model here presented ( Section 2.1.1 ). For metastillation with two liquid divisions the assign parameters are θ = 2 and β = 1 .
The molar fraction profiles generated by our methodology were in great agreement with those reported by Gouvêa et al. (20 0 0) and Gouvêa (1999) . The correlation coefficient was 0.999 for both processes. The reboiler and condenser temperatures were equal to those reported by Gouvêa et al. (20 0 0) and Gouvêa (1999) . In the work of Gouvêa et al. (20 0 0) , the heat flows in the reboiler and condenser were inverted; this was confirmed by simulations in Aspen Plus® for corresponding conventional columns. Considering this inversion, the heat flows obtained with our simulations were equal to those reported by Gouvêa et al. (20 0 0) for parastillation. Note that the MESH equations of Table 2 , with the correct substitution of the number of phase divisions ( β or θ equal to 2), are identical to the equations reported by Gouvêa et al. (20 0 0) and Gouvêa (1999) .
Finally, we validated the DWC simulations by comparing our results to those reported by Luyben (2013, pp. 355-369) for a conventional DWC with division just in the middle part of the column. For this example, the conventional DWC was simulated considering the adaptations in the MESH equations, suggested in the Section 2.1.2 . The column was fed with 3600 kmol/h (at 358 K) of a mixture containing: 30 mol% of benzene, 30 mol% of toluene, and 40 mol% of o-xylene. Product purities were 99 mol%, in the three flows of the column (bottom, side and top product). Chao-Seader physical properties were used by Luyben (2013) . Our simulations used the Peng-Robinson thermodynamic property package. For a fair comparison of the results, the configurations proposed by Luyben (2013) were also simulated in Aspen Plus® considering the Peng-Robinson thermodynamic property package. The columns were simulated in Aspen Plus® considering two configurations: the combination of four Radfrac columns, and one MultiFrac column with the Petlyuk configuration.
For the DWC validation, we used a column with 12 and 8 stages in the bottom and top sections, respectively, and 24 stages in each side of the middle section, plus one condenser and one reboiler. This configuration was accessed in the algorithm by assigning the parameters of Eqs. ( 3) -( 6) . Considering the number of stages in each section of the DWC, the input parameters were: nd = ( 12 + 1 ) = 13 , na = ( 13 + 1 ) = 14 , nb = ( 14 + 24 ) = 38 , and nc = ( 38 + 24 ) = 62 . The results generated by the new unified model were slightly different from those we obtained with Aspen Plus® using the column configurations reported by Luyben (2013) and the Peng-Robinson thermodynamic package. However, the molar fraction and temperature profiles were very similar in our Mat-Lab® and Aspen Plus® simulations. The differences are associated with the adopted considerations, and with the fact that the Aspen Plus® simulation, based on the work of Luyben (2013) , considers the combination of more than one column for simulating the DWC equipment, while in our new unified model we consider one unique column, as the real equipment should be. The differences between our MatLab® and Aspen Plus® results with each other and with Luyben's results were of the same magnitude as those between the two different Aspen Plus® procedures adopted by Luyben (2013) : the combination of four Radfrac columns or the Petlyuk setting. This means that the slight differences observed in the results should be attributed to the different configurations used for representing DWC equipment and to the differences in the thermodynamic packages used. In this way, we concluded that the new unified model was also validated for simulating dividing-wall columns.
For all configurations mentioned above, Newton's method converged. However, columns with larger number of stages and/or divisions required more iterations. This may be due to inadequate initial guesses. Additionally, the convergence is slower for problems with a larger number of variables. Since columns with parallel streams can accommodate more stages per equipment height, these columns usually have more stages than conventional columns. This means that having appropriate procedures to generate initial estimates for the simulation of columns with parallel streams is very important. Procedures to generate initial estimates of the variables in the case of conventional distillation columns are also useful for simulating alternative columns. These procedures include the linear interpolation indicated by Fredenslund et al. (1977) and the Bubble-Point method reported by Seader et al. (2011a) . The short-cut method presented by Meirelles et al. (2017) is also a good alternative for estimating initial values in the case of binary mixtures. This method can be used to estimate the internal flows, while the temperatures are assumed to be the bubble temperature of each stage. However, the converged traditional distillation ( β = θ = 1 ) seems to be a better initial guess for the alternative columns.
In the case of distillation with multiple internal streams we observed that a small number of internal divisions has no effect on the convergence of Newton's method. However, the equations become more difficult to solve when the number of partitions becomes larger. For columns with more than four divisions a careful selection of initial values was important for the practical convergence of Newton's method. We used the results from rigorous conventional distillation simulations as initial estimates for the alternative columns. The total number of stages in the conventional column should be equal to the total number of stages in the alternative columns, such that the total number of molar liquid and vapor flows and temperatures is the same in both columns. The mapping between stage numbers follows the one indicated in Fig. 1 . Other specifications, such as reflux ratio and distillate flow, are the same for both columns. The conventional distillation temperature profile can be directly used in the alternative columns as initial estimative. The component molar flow must be previously divided by the number of phase divisions. More details of this procedure are given in Section 3.4 .
For some column configurations Newton's method converges more slowly than for others. In these cases it may be necessary to use a sequence of simulations using configurations that are easier to converge to generate starting values for simulating a more difficult configuration. For example, simulations with larger reflux ratios exhibit better convergence properties. Therefore, one can perform a sequence of simulations, each one using a smaller reflux ratio than the previous one, and use the converged solution from the previous calculation as the initial estimate for the next. We also observed that simulations with smaller values of Murphree efficiency have better convergence properties, and a series of simulations with sequential Murphree efficiency values can also be used. Furthermore, simulations of columns with fewer components or fewer stages have better convergence properties than those for columns with larger number of components or stages.
Again, the corresponding results can be used as starting point for a next simulation, provided values for the additional variables are incorporated in the initial estimate file. For instance, in the case of a higher number of stages the replication of the values of the variables near the condenser and/or reboiler generate a good set of estimate values. In addition, it is important that variables like temperatures and molar flows obey given physical restrictions. For instance, component molar flows and temperatures must remain non-negative ( Fredenslund et al., 1977 ;Naphtali and Sandholm, 1971 ). In this case, if a computed molar flow or temperature is negative, it is reset to a small positive number.
For conventional columns and for para-and metastillation with two divisions we always obtained the same solution, even when Newton's method was started with different initial iterates. Specifically, we simulated different configurations, with θ and β less than or equal to 3, considering binary and multicomponent mixtures. Independent of the initial iterate, Newton's method always converged to the same solution for these columns. However, for columns with more than four streams, we observed that different procedures to generate initial iterates for Newton's method may lead to different solutions, given the same column specification. This is due to the high degree of non-linearity of the equations. Some of the computed solutions may have no physical meaning. However, it is possible that a system has more than one physical solution, i.e. multiple steady states ( Chavez et al., 1986 ). In this case, the processes may have more than one stationary point. This means that two equal-size columns operating under the same conditions -same feed conditions, same reflux ratio and same distillate flow -may operate at different temperature and composition profiles. Columns with multiple stationary points have already been observed in practice ( Gaubert et al., 2001 ;Güttinger et al., 1997 ).
To illustrate the effect of different initial estimates in the system solution, we simulated parastillation columns with 2-8 phase divisions. The saturated liquid feed was composed by methanol (50 mol%) and ethanol (50 mol%) at a flow rate of 450 kmol/h. The vapor distillate flow was 225 kmol/h. The columns height was fixed in 11.22 m and the number of stages were adjusted, for each number of phase divisions. The Murphree efficiency and the column pressure were 100% and 1 atm, respectively, for all stages. The Peng-Robinson thermodynamic package was adopted, for equilibrium and enthalpy calculus. The columns were configurated to obtain a top and bottom product with 95 mol% of methanol and ethanol, respectively.
Considering the above column specifications, the temperature in the condenser and in the reboiler should be around 65 °C (dew point) and 77 °C (bubble point), respectively. Using these temperatures to generate the initial estimate for the temperature profile based on a linear interpolation, as suggested by Fredenslund et al. (1977) , we obtained results consistent with conventional distillation, for all the divisions. However, using a linear interpolation between 80 and 100 °C as an initial iterate for Newton's method -an example of a poorly guessed initial profile estimate -the computed solution is not physically meaningful for 5 and 7 divisions, and the system did not converge for 4, 6 and 8 divisions. For 5 and 7 divisions, some stage temperatures were around 110 °C, which is higher than the reboiler temperature ( ∼77 °C). These solutions also presented stages with vapor flow tending to zero. A complete set of results are reported in the Supplementary Material, for a parastillation column with 7 vapor divisions.
Poor initial estimates seem to slow down the convergence of conventional distillations, but do not affect its final solution. Once the temperature profile is similar in conventional distillation and in parastillation columns, converged profiles of conventional distillation can be used as initial estimate for parastillation columns. We simulated the parastillation columns using as initial estimate the temperature profiles obtained from the simulation of the conventional distillation. For all cases, the results presented physical meaning and were equivalent to those obtained using linear interpolation between 65 °C and 77 °C as in the initial estimate. The use of converged profiles from conventional distillation as initial estimative for parastillation and metastillation is highly encouraged, especially in cases where one does not have a good initial estimate of the temperature profile. It is important to highlight that this example considers a simple and ideal distillations process. However, industrial cases often involve separation systems that are highly non-ideal. For these systems, a good initial estimative is even more important.
We illustrate the use of our proposed model via the simulation of parastillation and metastillation of alcoholic mixtures. This case study illustrates the use of the mathematical model presented here, but it also expands understanding of the paraand metastillation processes. The first example considers the binary ethanol/water mixture. Although we selected this system as a generic case study, it should be noted that this mixture is a simplified version of the system involved in the production of hydrous ethanol, which can be used in flexible fuel or ethanol-powered cars. In a second example, a multicomponent mixture is considered. This second approach is based on a real industrial case.
We first evaluated the minimum number of stages and minimum reflux ratio by the rigorous procedure described by Kister (1992) and Treybal (1981) . The simulated columns were maintained at a constant pressure of 1 atm and the Murphree efficiency was considered 70% for all stages, except for the reboiler, where the Murphree efficiency was 100%. The feed flow, at the temperature of 94 °C, was composed of 518.04 kmol water/h (9332.65 kg water/h) and 14.35 kmol ethanol/h (661.50 kg ethanol/h). The top ethanol concentration was fixed at 0.9300 (kg ethanol/kg total) at a rate of 17.00 kmol total/h. A partial condenser with distillate as vapor was considered. The thermodynamic properties were calculated by the NRTL model, with the parameters adjusted by Batista et al. (2012) .
In the case of metastillation columns, the feed was introduced at θ consecutive trays. Each tray received 1/ θ of the total feed flow, where θ is the number of liquid phase partitions. This procedure distributes the effects of the liquid feed to all internal flows, once the introduction of the feed in just one tray is equivalent to the addition of the entire feed amount in just one liquid internal flow. This division avoids the imbalance that can be induced by the external feed stream ( Biasi, 2016 ). The same reasoning can be extended to vapor feed streams in the case of parastillation columns, but in the present case study only liquid feeds were tested. Fig. 4 represents the minimum number of stages and minimum reflux ratio. The minimum number of stages to obtain an ethanolic concentration of 93.00 mass% in distillate was 15 stages in conventional distillation, 20 for metastillation with two phase divisions, and 23 and 33 for parastillation with two and three vapor divisions, respectively. The required number of stages to perform a specific separation, fixing the reflux ratio, increases with the increase of the number of phase partitions. However, the minimum reflux ratio, of 2.1, was the same for all processes. This result is in accordance with Meirelles et al. (2017) , who demonstrated that the minimum reflux ratio does not depend on the number of phase divisions.
Parastillation columns require more stages than conventional distillation to complete the same separation, although it is possible to allocate more stages per column height, as reported by Jenkins (1985) for two phase divisions. In fact, when one keeps the same tray spacing, it is possible to allocate approximately β additional parastillation stages per column height than in conventional distillation. The total number of stages (including the reboiler and the condenser) in parastillation columns ( N P/β S ) with the same height as the conventional ones can be calculated using the following Eq. ( 11) ,
where β is the number of vapor divisions and N D S is the total number of stages of the corresponding conventional distillation (including the condenser and the reboiler). A similar equation was reported by Meirelles et al. (2017) for parastillation columns with a distillation tray in the feed stage. In the present case, all trays are parastillation stages.
Fig. 5 shows the relation between the number of stages required for conventional distillation and for parastillation (with β= 2) to obtain the same overall separation with ethanol concentration in the distillate equal to 93.00 mass%. The solid line The number of parastillation stages corresponding to the number of distillation stages for columns with the same height were calculated by Eq. ( 11) .
b The column diameter corresponds to the diameter of the enriching section, slightly larger than the corresponding stripping section diameter.
represents parastillation columns with the same height of the corresponding conventional distillation ones. On the other hand, fixing the same reflux ratios of the corresponding conventional distillation, parastillation columns require lower numbers of stages than that indicated by the prior continuous line ( Fig. 5 ) and these columns also have a lower height than conventional ones. Canfield (1984) , evaluating a binary system under total reflux, concluded that parastillation columns requires 33% more stages than conventional ones. Similarly, Gouvêa (1999) , also working with binary systems, but with partial reflux ratios, observed an increase of 55% in the number of parastillation stages in comparison to conventional columns. These results indicate that it is always possible to obtain the same degree of separation by using a parastillation column with a smaller height than the conventional one, without any increase of the reflux ratio. Another possibility is to use a parastillation column with the same height as the conventional one, but operating with a lower reflux ratio. The first option leads to a reduction in the capital costs, while the second one reduces the operational costs. A combination of both alternatives is also possible.
The conclusions, obtained with the simplified binary distillation, can be extended to the multicomponent industrial plants for bioethanol distillation. This second case study is based on the column configurations studied by Batista et al. (2012) . The columns were maintained at a constant pressure of 1 atm and 70% Murphree efficiency was considered for all stages, except for the reboiler, where the efficiency was 100%. The feed flow of 532.46 kmol/h (9997.14 kg/h), at a temperature of 94 °C, was composed of 518.04 kmol water/h (9329.92 kg water/h), 14.36 kmol ethanol/h (661.37 kg ethanol/h), 1.62.10 -2 kmol isoamyl alcohol/h (1.42 kg isoamyl alcohol/h), 8.73.10 -4 kmol ethyl acetate/h (7.69.10 -2 kg ethyl acetate/h), and 4.72.10 -2 kmol glycerol/h (4.35 kg glycerol/h). The top ethanol concentration was fixed in 0.9300, in mass fraction, at a rate of 17.00 kmol total/h ( ∼707 kg/h). A partial condenser with distillate as a vapor stream was considered. The liquid side stream, rich in isoamyl alcohol, was set at 0.30 kmol/h ( ∼6 kg/h), located three stages above the feed tray. The thermodynamic properties were calculated by the NRTL model with the parameters adjusted by Batista et al. (2012) .
The investigation included simplified estimations associated with the column construction, involving equipment dimensioning and costs. The column diameter was calculated based on the methodology described in Treybal (1981) . The capital cost included heat exchangers, column shell and trays values and was calculated based on equations reported by Douglas (1988) and Kiss (2013) , considering a M&S index of 1468.6. It was adopted a heat-transfer coefficient of 0.852 and 0.568 kW/(K .m 2 ) for the condenser and reboiler, respectively. The column height was calculated considering the tray spacing of 0.6 m, with a 10% of tolerance. For all columns, the material of construction was stainless steel (Fc = 3.67). The energy cost considered a steam price of $8.22/GJ. The total annual cost (TAC) considers the sum of the capital and operational costs, normalizing the capital costs by a payback period of 3 years and the operational costs by 8760 operating hours per year. The results of this investigation are presented in Table 4 .
A conventional distillation column with 36 stages requires a reflux ratio of 4.5 to achieve a distillate with at least 93 mass% of ethanol. The same operation conducted in a parastillation column, with the same number of stages and same reflux ratio, produces a distillate with ethanol concentration below the specification ( < 93 mass%). For the same reflux ratio, 49 parastillation stages are required to achieve the same degree of separation obtained in the conventional distillation with 36 stages. This parastillation column presents 36% more stages than the conventional one, although it has a height 33% smaller than the conventional column, contributing to a reduction of the capital costs and TAC of 16% and 9%, respectively. The capital costs are reduced using columns with smaller heights compared to the conventional ones, while the operational costs can be reduced by fixing the column height and decreasing the reflux ratio. According to Eq. ( 11) , a distillation column with 36 stages is comparable in height to a parastillation column ( β= 2) with 69 stages. This parastillation column still represents a reduction of capital costs due to the smaller column diameter caused by the lower reflux ratio. Furthermore, the alternative parastillation column (with β= 2) allowed a reduction of the reflux ratio equal to 24%, corresponding to a decrease of 19% in operational costs, when compared with the conventional distillation column ( Table 4 ). The total annual cost (TAC) can be decreased by 15% by substitution of the conventional distillation column with a parastillation column of the same height.
In this work, we propose a generic mathematical model for simulation of columns with parallel streams, including parastillation, metastillation and divided-wall columns. The initial challenge was the formulation of a unique set of MESH equations able to describe, with no further implementation, those different distillation processes. The model was implemented and validated considering different processes reported in the literature. Results showed good convergence and were in great agreement with previous literature results. The proposed approach displays several advantages compared to conventional methods. For instance: (i) the use of the traditional MESH equations avoids the simplifying assumptions made in short-cut methods; (ii) the possibility of simulating a comprehensive set of processes with a unique mathematical model, that is simply adjusted to the chosen column configuration by specifying a few parameters; (iii) fast simulation and convergence for different processes; (iv) the possibility of simulating column configurations not covered by conventional simulators; and (v) the possibility of investigating and comparing different column configurations in order to support the design of new, energy efficient and capital cost reducing columns.
Additionally, we applied the proposed model to the distillation of alcoholic mixtures of industrial interests. The case study discussed in the article indicates that alternative columns, especially the parastillation equipment, allows substantial reductions in the reflux ratio or in the equipment cost by decreasing its height, when compared with conventional distillation.