Design of flexible heat exchanger network for multi-period operation

Design of flexible heat exchanger network for multi-period operation

Chemical Engineering Science 61 (2006) 7730 – 7753 www.elsevier.com/locate/ces Design of flexible heat exchanger network for multi-period operation W...

385KB Sizes 17 Downloads 136 Views

Chemical Engineering Science 61 (2006) 7730 – 7753 www.elsevier.com/locate/ces

Design of flexible heat exchanger network for multi-period operation W. Verheyen, N. Zhang ∗ Centre For Process Integration, CEAS, The University of Manchester, PO Box 88, Manchester M60 1QD, UK Received 19 April 2006; received in revised form 9 August 2006; accepted 13 August 2006 Available online 30 August 2006

Abstract Heat exchanger networks (HENs) increase heat recovery from industrial processes by matching hot and cold streams to exchange heat and reducing utility consumption. The design of HENs is a very complex task which generally involves mixed-integer non-linear programming (MINLP). This work evaluates and critically compares existing HEN design methods. It then presents a systematic methodology in the design of HENs under multiple periods of operation. The model presented in this work is a superstructure-based MINLP model which minimises the total annualised cost containing heat exchanger area cost and utility costs. The model is based on the superstructure by Yee and Grossmann [1990. Simultaneous optimisation models for heat integration—II, heat exchanger network synthesis. Computer & Chemical Engineering 14(10), 1165–1184], which was later formulated for multiple periods by Aaltola [2002. Simultaneous synthesis of flexible heat exchanger network. Applied Thermal Engineering 22, 907–918]. It includes a multi-period simultaneous MINLP model to design the HEN structure, and an NLP model to improve the solution and allow for non-isothermal mixing. Modifications to Aaltola’s model include the use of maximum area per period in the area cost calculation of the MINLP objective function, and the removal of slack variables and weighed parameters from the existing NLP improvement model. The new model has been applied to one industrial case study, demonstrating that the new combined MINLP–NLP model can obtain better solutions by not relying on the average area assumption in the MINLP stage. 䉷 2006 Elsevier Ltd. All rights reserved. Keywords: Heat exchanger network; Multi-period; Mathematical modelling; Design

1. Introduction to HEN design The world’s increasing energy prices and limited fossil fuel reserves push industries towards better utilisation of energy and increased energy savings. Heat exchanger networks (HENs) form an important part in the overall chemical processes since they utilise available heat from processes in the most economic way. The design of the HEN links the process flowsheet with the utility system and generally involves a large fraction of both the overall plant capital cost but, most important, of the overall plant running costs in terms of energy requirements. Therefore, a proper design for heat integration is one of the key factors to a profitable industry.

∗ Corresponding author. Tel.: +44 161 3064384.

E-mail address: [email protected] 0009-2509/$ - see front matter 䉷 2006 Elsevier Ltd. All rights reserved. doi:10.1016/j.ces.2006.08.043

The HEN synthesis problem has been formulated for the first time by Masso and Rudd (1969). Since then, the synthesis and design of HENs has received considerable attention in lots of different aspects. The synthesis and design of HENs has been very well studied during the past few decades. The goal is to recover energy from the process in the most economical way by matching hot process streams to cold process streams. Any energy still needed is then supplied by hot and cold utilities. The main focus of heat exchanger design strategies is in the design of HENs for fixed parameters. However, it is possible that there are significant changes in the environment of a plant. There are two major groups of possible changes of parameters. First, uncertainties are deviations of a certain parameter around one nominal value. These may be caused by outside environmental conditions or by poor control of this parameter. Second, periodical changes in operating conditions can be either seasonal changes or intentional changes to the system operating conditions. The flexible HENs that can remain operable and

W. Verheyen, N. Zhang / Chemical Engineering Science 61 (2006) 7730 – 7753

optimum under these variations then become “resilient” HENs (uncertainty parameters) or “multi-period” HENs (periodical changes). This paper presents a renewed summary and comparison between current heat exchanger design methods. The main objective is to develop a systematic methodology for the design of multi-period HENs, based on a suitable single period model, which is both accurate with reasonable solving times. The developed model should be efficient for small-scale problems and should be capable of providing good solutions to larger scale problems as well. 2. Literature review Since the first heat exchanger design methodologies were published a few decades ago, several groups of design strategies evolved. The earliest methods relied on a combination of thermodynamic principles and heuristic methods. Later more systematic approaches were developed, from pinch analysis to mathematical programming to stochastic methods. This chapter reviews and compares the different HEN design methodologies, including the design of flexible and multi-period HENs. The last section contains a summary of the strengths and weaknesses of each group of methods. 2.1. Pinch analysis and thermodynamics 2.1.1. Pinch analysis concept The first systematic methods in the HEN design were introduced by Hohmann (1971) and later recalled and refined by Linnhoff and Flower (1978a,b). Both recognised that it is possible to identify targets that predict the maximum amount of energy that can be recovered, as well as the minimum number of units of the HEN. This method was named the ‘pinch design’ method and allows targets to be set which are independent of the structure of the network and are only based on stream data. Pinch analysis is now globally known as a thermodynamic analysis tool which is based on the first and second law of thermodynamics. Pinch analysis sets design targets in terms of minimum utility consumption, minimum number of units, and minimum area and cost for a given minimum temperature difference. Using pinch analysis, the problem is reduced to two smaller regions, and HEN can be designed above the pinch and below the pinch separately. 2.1.2. Heat exchanger design method For graphical purposes, the grid diagram is used for the presentation of the HEN. Hot streams and cold streams are represented as straight lines going from left to right and from right to left, respectively. After the pinch has been located by the problem table algorithm (either manually or in the form of a linear optimisation program), one can then separate the grid diagram into two separate regions, allowing the design of two smaller scale networks. Pinch design of each subsystem starts at the location of the pinch, where there are more limiting constraints. According to pinch analysis, a HEN can be designed

7731

by the following rules: • No heat exchanger is allowed to exchange heat across the pinch. • No side of any heat exchanger can have a temperature difference below DTmin . • Cold utilities are only placed below the pinch, hot utilities above the pinch. Any violation of these rules results in increased utility consumption. A more detailed description of the HEN design procedure under the pinch analysis can be found in the literature (Smith, 2005). However, the ultimate objective is to find a HEN which does not use minimum hot and cold utilities, but which has minimum total annualised cost (TAC). As the level of energy recovery increases (DTmin decreases) and the use of utilities is minimised, the temperature differences, which are the driving forces to the heat exchange, are decreased as well. This results in a sharp increase in heat exchanger area. Also a higher level of energy recovery might feature a greater number of units. There is a clear trade-off between number of matches, area cost (capital costs) and utility cost (operating costs). To minimise capital costs, Linnhoff (1979) established a targeting method to obtain the minimum number of matches for a given set of streams and minimum approach temperature. The method relies on detecting any heat exchange loops or utility paths and minimising the number of exchangers by removing all loops and paths. If a pinch exists, this procedure can be applied to the two separate sections independently. Pinch analysis is a powerful tool for screening of existing processes for retrofit as well as providing targets of possible heat recovery in grassroots designs. The main strength of pinch analysis is that it can provide targets as well as provide graphical interpretation and visualisation of the problem, which makes it easier to understand. Many sequential approaches rely on the pinch analysis in different stages for obtaining minimum hot and cold utility demands and networks which feature the minimum number of units. However, it does not take into account heat transfer coefficients and heat exchanger areas properly and therefore might be misleading in some cases. Also, the composite curves and grand composite curve are assumed to be straight lines, meaning that the heat capacity flow rate needs to be constant over the entire temperature range. 2.2. Sequential approaches When there first was the need to design HENs in mathematical programming, the optimisation techniques available were not adequate to solve the optimisation problems, and the performance of the existing computers did not allow the solving of more complex models. Pinch theory provided the right background for a number of sequential design procedures based on problem decomposition. The sequential design consists of the three following sub-problems (Biegler et al., 1997): • minimum utility costs; • minimum number of units; • minimum investment cost of networks.

7732

W. Verheyen, N. Zhang / Chemical Engineering Science 61 (2006) 7730 – 7753

Papoulias and Grossmann (1983) proposed a transhipment model to predict the minimum utility consumption using linear programming (LP) and the minimum number of units using mixed-integer linear programming (MILP). The sources in the model represent the hot streams and hot utilities. The destinations represent the cold streams and cold utilities. Each temperature interval (warehouse) accepts heat from the hot streams, supplies heat to the cold streams, and passes on a heat residual to the following temperature interval. This heat residual is the excess heat which could not be utilised by the cold streams in that temperature interval. Since the heat flows from the sources and to the destinations in each temperature interval are fixed, the only variables are the heat residuals which are passed on into the next interval. Therefore the proposed optimisation model is linear and minimises the use of hot and cold utilities while keeping all heat residuals positive. There is a strong analogy between the transhipment model and the grand composite curve in the pinch design method, which is also based on heat residuals. Other transportation models have been developed for these purposes (Cerda et al., 1983). For design of networks which feature the minimum number of units, the same two models are most commonly used. This design is formulated as a MILP, in which the integers represent the existence of a match between two streams. In 1986, Floudas et al. developed a model which can minimise TAC in a two-stage procedure. First the MILP model minimises the number of heat exchangers. Then this fixed heat exchanger structure is used within a non-linear programming (NLP) model to obtain the minimum TAC of the network. Since an MILP program can have multiple solutions, it is important to obtain all possible mixed-integer solutions and perform NLP programming on all obtained solutions. Even though this decomposition reduces the size of the problem significantly, there is no rigorous analysis of the trade-off between energy, number of units and area, and therefore sequential methods can lead to sub-optimal solutions. For example, HENs which could be optimum in terms of TAC, but feature more than the minimum number of heat exchangers, are excluded from the search space. This is a major drawback for sequential approaches. The most recent sequential approaches were proposed by Zhu (1995) and Zhu et al. (1995a,b,c) by decomposing the problem into enthalpy intervals. Later Zhu (1997) proposed an automated synthesis approach, consisting of an MILP model to select the matches and an mixed-integer non-linear programming (MINLP) model for determining the cost-optimal network. 2.3. Simultaneous approaches Simultaneous optimisation techniques consider all the tradeoffs mentioned previously into one rigorous model, which in most cases is an MINLP formulation. Simultaneous approaches have shown to be superior to sequential approaches in most cases. The problems that arise with simultaneous approaches are the complexity of the model and the difficulty in attaining the optimal solution. Some assumptions need to be made to

make the problem become solvable, which are: • which superstructure to use; • assumptions to reduce the model size (number of equations, variables); • assumptions to reduce the model complexity (non-linear, non-convex). Ciric and Floudas (1991) combined the transhipment model of Papoulias and Grossmann (1983) for the match selection with the NLP model to determine the minimum investment cost network (Floudas et al., 1986) into one MINLP formulation for a specified minimum temperature difference. This method consists of one single stage optimisation in which all variables are optimised simultaneously. The hyperstructure used by the authors is shown in Fig. 1. While this hyperstructure embeds all different alternative structures, it involves more non-linear heat balance constraints and is therefore more difficult to solve. At the same time, a stage-wise simplified superstructure formulation has been developed by Yee et al. (1990) and extended by Yee and Grossmann (1990). This method does not rely on any division into temperature or enthalpy intervals and features linear constraints from the following assumptions: • • • •

isothermal mixing, no split stream flowing through more than one exchanger, utilities at the end of the superstructure and no stream bypass.

The superstructure presented by Yee et al. (1990) is shown in Fig. 2. At each stage in the superstructure, each hot stream is split into the number of cold streams and one possible heat exchanger is placed at every branch. As shown, utilities are located at the end of the superstructure. The number of stages can be decided by the designer and is normally chosen to be equal to the minimum of the number of hot and cold streams. After increasing the number of stages a few times, there will be no more improvement in the objective value. Even though the set of constraints is linear, the objective function remains non-linear, non-convex, and of large size. The assumptions to formulate this last superstructure reduce the problem to a manageable size, but with the loss of some possible HEN configurations. Floudas (1995) has observed that a number of structures cannot be obtained by the superstructure of Yee et al. (1990). These are shown in Fig. 3. As can be seen in Fig. 3, the structures excluded are the following: • splits with two or more exchangers in series on one branch; • bypass streams which are only feasible for non-isothermal mixing; • any combination of the first two structures. The main motivation for using a simplified superstructure is to end up with a mathematical model that features only linear constraints while all non-linearities appear only inside the

W. Verheyen, N. Zhang / Chemical Engineering Science 61 (2006) 7730 – 7753

1

4

H1-C1

5

6

8

C1

2

9

H2-C1

7733

7

12

10

11

13 18

3

14

H3-C1

17

15

16

Fig. 1. Hyperstructure of Ciric and Floudas (1991).

Stage 1

Stage 2

Stage k

1,1,1

1,1,2

1,1,k

1,2,1

1,2,2

1,2,k

Steam

H1

Steam

H2

.. .

.. .

.. .

1,n,1

1,n,2

1,n,k

2,1,1

2,1,2

2,1,k

2,2,1

2,2,2

2,2,k

.. .

.. .

.. .

2,n,1

2,n,2

2,n,k

.. .

.. .

.. .

m,1,1

m,1,2

m,1,k

m,2,1

m,2,2

m,2,k

C1

CW

C2

CW

Steam

Hm

.. .

.. .

.. .

m,n,1

m,n,2

m,n,k

Cn

CW

Fig. 2. Simplified stage-wise superstructure of Yee et al. (1990).

objective function. Even though some HEN structures cannot be generated by the model, Yee and Grossmann (1990) illustrated that good HEN structures can be obtained. Due to the relative ease of obtaining solutions while still incorporating most of the possible configurations, the simplified superstructure could prove very effective for the design of multi-period HENs, and this superstructure can be used and extended to the multi-period model. Bjork and Westerlund (2002) evaluated the assumption of isothermal mixing in their recent publication, and proposed a methodology based on the same superstructure without the assumption of isothermal mixing. Removing this assumption greatly increases the number of variables and introduces nonlinear and non-convex equations as constraints. Therefore, a

global optimisation procedure was suggested which creates and solves convexified sub-problems. This global optimisation method makes it possible to reach the global optimum under the assumptions of constant heat capacity flow rates and heat transfer coefficients. The limitation of this method, however, is when the problem is of large size (such as multiple periods), it might be impossible to solve it globally. 2.4. Multi-period and flexible HEN design 2.4.1. Flexibility criteria and flexible HEN design While good design methodologies exist for the design of HENs with fixed parameters, changes in operating parameters might cause deviation from optimality or even cause control

7734

W. Verheyen, N. Zhang / Chemical Engineering Science 61 (2006) 7730 – 7753

θ2 Feasible region R

Ψ(d,θ) = 0

θ2N

θ1N

θ1

Fig. 4. Feasible region of operation.

Fig. 3. Structures excluded from simplified HEN superstructure.

and safety problems. Therefore, a good HEN should be optimal for a set of nominal conditions as well as provide the flexibility to handle a range of operating conditions. Changes in operating conditions might originate from uncertainty in parameters during operation or from different operating conditions, which then refers to multi-period operation. This section reviews the literatures on the design of flexible HENs under variable operating conditions, both on uncertainty as well as multi-period approaches. The literature presented is limited to grassroots design of HENs. Marselle et al. (1982) firstly developed the concept of resilient HENs which can tolerate uncertainties in temperatures and flow rates. This method, however, needs a manual combination of a series of optimal designs under different worst-case scenarios, making the application of this method to large-size problems practically impossible. To measure the degree of flexibility of a HEN, Saboo et al. (1985) proposed the resilience index (RI). The RI can be used to compare between different options on a quantitave basis. Swaney and Grossmann (1985a,b) introduced a flexibility index, which represents a maximum deviation of an uncertain parameter while still remaining within the feasible region. The flexibility index provides a basis for different designs to be

compared as well as providing the necessary information about which critical points in the feasible region that limit the design. This concept is illustrated in Fig. 4. The point in the middle represents the nominal operating conditions for parameters. The surrounding constraints mark the region R for which there is a feasible solution. The rectangles inside show the maximum deviation of the two parameters for which the network remains feasible. Tantimuratha et al. (2001) proposed a two-stage design procedure, where the first stage consists of area targeting using area target models developed by Briones and Kokossis (1999). For a selected set of matches from the first stage, the network is optimised with an iterative procedure to minimise annual costs. This method, however, could still lead to sub-optimal solutions because it relies on decomposition. 2.4.2. Multi-period HEN design For multi-period operations, the different operating modes resemble the corner points of the feasible region. Floudas and Grossmann (1986) introduced a multi-period MILP model based on the transhipment model of Papoulias and Grossmann (1983). This model for minimum utilities and minimum number of matches uses pinch points at each operating period. The next stage is a reformulated NLP model that develops the network configuration, presented by Floudas and Grossmann (1987). Here again, the decomposition of the problem into different stages significantly reduces the size of the problem, but does not take into account the trade-off between area, number of units, and energy rigorously. The model proposed by Papalexandri and Pistikopoulos (1994a,b) simultaneously explores alternatives in an MINLP problem. This of course limits the size of the problem to relatively small-scale problems. Aaltola (2002) proposed a model which simultaneously optimises the multi-period MINLP problem for minimum costs and flexibility, without relying on sequential decomposition. This model, based on the stage-wise HEN superstructure

W. Verheyen, N. Zhang / Chemical Engineering Science 61 (2006) 7730 – 7753

representation of Yee et al. (1990), has been successfully applied to industrial HEN problems. Since the design of multiperiod heat exchangers involves the use of the same heat exchangers for different heat loads, this will generally involve bypassing part of a stream or even an entire stream from the heat exchanger unit. The formulation of Aaltola (2002) allows the elimination of bypass modelling, which would introduce non-linear constraints into the formulation. Chen and Hung (2004) proposed a three-step approach for designing flexible multi-period HEN, which is based on the stage-wise HEN superstructure representation of Yee et al. (1990) and the mathematical formulation of Aaltola (2002). The authors introduced the maximum area consideration in the objective function, and decomposed the problem into three main iterative steps: simultaneous HEN synthesis, flexibility analysis, and removal of infeasible networks. The number of iterations required in this approach is of concern for industrial problems.

MINLP Model

Solve Relaxed NLP Problem

All Binaries Integral?

Solve MILP Master Problem

Solve NLP Subproblem

2.5. Selection of deterministic optimisation algorithms Most sequential approaches offer another advantage, which is that at each stage it is possible or even guaranteed that the global optimum can be obtained. Therefore, the application of existing MILP and NLP solvers is sufficient to solve the HEN design problem. Simultaneous approaches, however, are generally formulated as MINLP problems, which have an extensive number of local optima that can be readily obtained. While introducing the HEN design problem based on their superstructure, Yee and Grossmann (1990) discussed the complexity of the model and proposed to solve MINLP formulation by using the combined penalty function and outer-approximation method, firstly published by Viswanathan and Grossmann (1990). This solver has been incorporated under the name DICOPT as the standard GAMS (GAMS Development Corporation, 2001) solver for MINLP problems. This is the reason why this solver has been used for solving the MINLP problems formulated in this work. In their paper, Viswanathan and Grossmann give a description and comparison between different existing MINLP solver methods before proposing their combined method. The newly proposed method was based on the outerapproximation/equality-relaxation (OA/ER) algorithm, which relies solving an NLP sub-problem and MIP master problem sequentially in a loop structure. As described by its authors, the newly developed solver includes a new MILP master problem which incorporates a penalty function method that penalises any violations to the non-linear constraints. This new solver does not need an initialisation of the integer variables for the master MILP problem, because the MILP problem starts with the solution of a relaxed NLP problem. A summary of the solver algorithm is shown in Fig. 5. As has been discussed by the authors, this solver method has shown to often lead to the global optimum. Also, the chance of reaching the global optimum is increased if a good initial solution can be obtained by the relaxed NLP problem. Finally, one should keep in mind that this solving strategy cannot guarantee

7735

A Higher NLP Soln Value?

Optimal Solution obtained

Fig. 5. Combined penalty-function/outer-approximation method.

that the global optimum is reached, but it can provide satisfying results in terms of minimising the total annual cost of a HEN. 2.6. Summary Until now, there exist two multi-period HEN design models. The first one was formulated by Floudas and Grossmann (1986) and is based on the sequential approach. Recently, the first simultaneous multi-period HEN design approach was formulated by Aaltola (2002). The sequential multi-period model fails to take into account the trade-offs resulting from different objective functions, and could lead to sub-optimal problems. However, it is less demanding in terms of computer solving time. The simultaneous multi-period model takes into account those trade-offs in a rigorous way, but relies on some simplifications in order to linearise the set of constraints. The simultaneous model proposed in the following section tries to take away some of those assumptions in order to reach better solutions. 3. Improved formulation for multi-period HEN design Based on the information presented in the previous section, the simultaneous design approach has been chosen because its rigorous framework can obtain optimal solutions without any major decomposition. The simplified superstructure formulated by Yee et al. (1990) has been extended to multi-period HEN de-

7736

W. Verheyen, N. Zhang / Chemical Engineering Science 61 (2006) 7730 – 7753 Stage 1

Stage 2

1,1,1

1,1,2

1,1,k

1,2,1

1,2,2

1,2,k

Stage k

Steam

H1

.. .

.. .

.. .

1,n,1

1,n,2

1,n,k

2,1,1

2,1,2

2,1,k

2,2,1

2,2,2

2,2,k

C1

CW Steam

H2

.. .

.. .

.. .

2,n,1

2,n,2

2,n,k

.. .

.. .

.. .

m,1,1

m,1,2

m,1,k

m,2,1

m,2,2

m,2,k

C2

CW

Steam

Hm

.. .

.. .

.. .

m,n,1

m,n,2

m,n,k

Cn

CW

Fig. 6. Simplified stage-wise superstructure by Yee et al. (1990).

sign by Aaltola (2002). In this section, Aaltola’s model has been presented, as well as a set of modifications made to this model to obtain more accurate solutions. These modifications make the model more detailed at the expense of introducing nonlinear functions into the set of constraints. The overall model consists of two stages: an MINLP simultaneous model which optimises the HEN structure by assuming iso-thermal mixing, and an NLP improvement model that takes away this assumption. In contrast to the model published by Aaltola (2002), the proposed model includes the maximum area cost into the objective function at both stages. 3.1. Simultaneous MINLP model 3.1.1. Superstructure The simultaneous models found from the literature all rely on some form of superstructure. The formulation of this superstructure is critical to the optimisation problem; simple superstructures are robust and easy to solve, at the expense of removing some feasible structures from the solution space. On the other hand, complex and rigorous superstructures which include all possible HEN structures in the solution space might be very difficult to solve and may cause very high computational times. The MINLP model used in this work relies on the stagewise simplified HEN superstructure developed by Yee et al. (1990) and extended to multiple periods by Aaltola (2002). This simplified superstructure consists of a number of stages, in which different matches between hot streams and cold streams are allowed. For each stage, each hot stream is split into the equivalent number of cold streams and is directed to a heat exchanger which matches this hot stream to each cold stream. At the end of the stage, a mixer combines all split streams back into one single hot stream, which will be the feed into the following stage. At each stage, the inlet and outlet temperatures of hot and

Hot stream i

t(i,k+1,p)

t(i,k,p)

H1 H2

C1

C2 Cold stream j

t(j,k,p) Stage k Period p

q(i,j,k,p)

t(j,k+1,p)

Fig. 7. Reduced version of superstructure for one stage and four streams.

cold streams are treated as variables, while hot utilities and cold utilities are situated at the end of the HEN superstructure. From this formulation, the superstructure does not depend on any division or partitioning into temperature intervals or enthalpy intervals, and the pinch point is not required. Fig. 6 shows again the simplified superstructure, which is easily extended to multiple periods. A reduced version incorporating two hot streams and two cold streams has been shown in Fig. 7 to identify all necessary variables. This structure consists of: • four heat exchangers per stage plus one utility per stream; • one temperature variable per stage per stream; • one binary variable assigned per heat exchanger. The assumption of isothermal mixing at the end of each stage greatly simplifies the formulation of the model. It specifies that the outlet temperature of each heat exchanger in each branch

W. Verheyen, N. Zhang / Chemical Engineering Science 61 (2006) 7730 – 7753

is the same and therefore also equals to the exit temperature of the stream for that stage. It removes a large portion of the temperature variables and also removes a set of non-linear mixing equations from the constraints. Moreover, in the way there is no need to define the flow rates in each branch as extra variables, and they can be calculated from the solution by energy balances across that heat exchanger. As a result of this assumption, the feasible region is bounded by a set of linear constraints, which improves the optimisation procedure and making the model more robust and easier to solve. The assumption of isothermal mixing has been investigated in detail by Bjork and Westerlund (2002). They also suggested methods of removing the isothermal mixing assumption while convexifying the non-convex terms in the model. In this way, it is possible to obtain a global optimum which satisfies the set of convex constraints. The main limitations of the assumption are the following (Floudas, 1995): 1. Resulting HEN structure may feature more exchanger units than required if a stream splitting takes place. 2. HEN structures that are only feasible with non-isothermal mixing may be excluded. 3. HEN structures in which a split stream goes through several exchangers in series are neglected, as was already shown in Fig. 1. This simplified superstructure is suitable for multi-period HEN design because it can provide a rigorous model without decomposition, which is manageable in size, and can provide good results. To overcome the assumption of isothermal mixing, Aaltola (2002) suggested an NLP improvement model to allow non-isothermal mixing. 3.1.2. Area calculations The heat transfer area for each individual heat exchanger can be calculated from the following information: 1. the heat load of the particular match; 2. the inlet and outlet temperatures of both hot and cold streams; 3. the overall heat transfer coefficient for the given pair of streams A(i, j, k, p) =

q(i, j, k, p) , LMTD(i, j, k, p) · U (i, j )

(1)

where i is the set of hot streams, j the set of cold streams, k the stage number, and p the period of operation. The log mean temperature difference (LMTD) is calculated as follows: LMTD(i, j, k, p) [ti(i, k, p)−tj(j, k, p)]−[ti(i, k+1, p)−tj(j, k+1, p)] = . ln((ti(i, k, p)−tj(j, k, p))/(ti(i, k+1, p)−tj(j, k+1, p))) (2) The variables are illustrated in Fig. 8.

ti(i,k,p)

St age k

7737

ti(i,k+1,p)

Heat Exchanger tj(j,k,p)

tj(j,k+1,p)

Fig. 8. Illustration of variables to calculate LMTD.

When the temperature differences on both sides of the heat exchanger are equal, the calculation of LMTD might cause numerical difficulties as a result of division by zero. For this reason, a number of different approximations are generally used, which are: 1. average LMTD, 2. Paterson approximation, 3. Chen approximation. The average LMTD just takes the average of both temperature differences: LMTD(i, j, k, p) [ti(i, k, p)−tj(j, k, p)]+[ti(i, k+1, p)−tj(j, k+1, p)] = . 2 (3) This approximation is very rough and will predict large errors when the temperature differences at both sides of the heat exchanger are not of similar size. The error is a large underestimation of the true heat transfer area. The approximation published by Paterson (1984) is stated as follows: LMTD(i, j, k, p) = 23 · [(ti(i, k, p) − tj(j, k, p)) · (ti(i, k + 1, p) − tj(j, k + 1, p))]0.5 +

1 6

· [ti(i, k, p) − tj(j, k, p)]

+ [ti(i, k + 1, p) − tj(j, k + 1, p)].

(4)

Using Paterson for calculating LMTD generally tends to slightly under-estimate the required heat transfer area. The approximation published by Chen (1987) is stated as follows: LMTD(i, j, k, p)  = (ti(i, k, p)−tj(j, k, p)) · (ti(i, k+1, p)−tj(j, k+1, p))  (ti(i, k, p)−tj(j, k, p))+(ti(i, k+1, p)−tj(j, k+1, p)) 1/3 . · 2 (5) Chen approximation for calculating LMTD generally tends to slightly over-estimate the actual heat transfer area. The benefit of using this equation is that it will predict a zero value for

7738

W. Verheyen, N. Zhang / Chemical Engineering Science 61 (2006) 7730 – 7753

To determine the temperatures at each stage, stage-wise heat balances are necessary, which define the total amount of heat rejected or received by a stream in any stage:

Comparison of Different LMTD approximations 120

LMTD (C)

100

[ti(i, k, p) − ti(i, k + 1, p)] · F i(i, p) =

80

Log LMTD average Paterson Chen

60 40



q(i, j, k, p),

j ∈CP

k ∈ ST, i ∈ HP, p ∈ PR,

(8) 

[tj(j, k, p) − tj(j, k + 1, p)] · Fj (j, p) =

q(i, j, k, p),

i∈H P

20

k ∈ ST, j ∈ CP, p ∈ PR.

0 0

At both sides of the superstructure, the variable temperatures need to be assigned by the inlet and outlet temperatures of all streams. The inlet temperatures of the hot streams will be at temperature location k = 1, and the inlet temperatures of the cold streams will be at temperature location k = NOK + 1:

50 100 150 200 250 outlet temperature difference (C) Fig. 9. Comparison of LMTD approximations.

LMTD if either side of the heat exchanger has a zero temperature difference. Fig. 9 summarises the accuracy of the different LMTD approximations while varying the outlet temperature difference of the heat exchanger. It confirms that both the Paterson approximation and the Chen approximation perform well, while Paterson slightly over-estimates the true LMTD and Chen slightly under-estimates the true LMTD. In this work, Paterson approximation has been used to calculate LMTDs. This formulation tends to be slightly easier to solve under more extreme values of the variables. 3.1.3. Multi-period MINLP model This multi-period MINLP model minimises the TAC. The costs included are capital costs for heat exchanger area and unit and operating costs for utility consumption. The match between two streams is denoted by a binary variable. This binary variable takes the value of 1 if a heat exchanger exists in any of the operating periods, and takes the value of 0 if no heat exchanger exists in any of the operating periods. The continuous variables to be optimised are stream temperatures for each stage and heat loads for each match in each period. The constraints consist of a series of equality constraints representing the heat balances at different levels, and a number of inequality constraints for assessing feasibility and defining logical operations. The overall stream heat balances make sure that each stream receives the amount of heating or cooling that is required for each period: [Tiin (i, p) − Tiout (i, p)] · F i(i, p)   = q(i, j, k, p) + qcu (i, p),

i ∈ HP, p ∈ PR,

k∈ST j ∈CP

(6) [Tjin (j, p) − Tjout (j, p)] · Fj (j, p)   = q(i, j, k, p) + qhu (j, p),

(9)

j ∈ CP, p ∈ PR.

k∈ST i∈H P

(7)

Tiin (i, p) = ti(i, 1, p),

i ∈ HP, p ∈ PR,

Tjin (j, p) = tj(j, NOK + 1, p),

j ∈ CP, p ∈ PR.

(10) (11)

To make sure that there is a monotonic decrease or increase in temperature from the inlet to the outlet temperature, the following constraints are used. Also the temperatures at the outlet of the superstructure are limited by the stream outlet temperatures. If no utility exists from one given stream, then the outlet temperature of the superstructure equals to the stream outlet temperature. If utilities exist for the given stream, then the superstructure outlet temperature can be either greater than the stream outlet temperature (cold utility) or less than the stream outlet temperature (hot utility): ti(i, k, p) ti(i, k+1, p),

k∈ST , i∈HP, p∈PR,

(12)

tj(j, k, p) tj(j, k+1, p),

k∈ST, j ∈CP, p∈PR,

(13)

Tiout (i, p) ti(i, NOK + 1, p), Tjout (j, p)tj(j, 1, p),

i ∈ HP, p ∈ PR,

j ∈ CP, p ∈ PR.

(14) (15)

Similar to the energy balances for the main set of heat exchangers in the superstructure body, there are energy balances for the utility matches, forming the relationship between the outlet temperatures of the superstructure and the stream outlet temperatures: [ti(i, NOK + 1, p) − Tiout (i, p)] · F i(i, p) = qcu (i, p), i ∈ HP, p ∈ PR,

(16)

[Tjout (j, p) − tj(j, 1, p)] · Fj(j, p) = qhu (j, p), j ∈ CP, p ∈ PR.

(17)

The following logical constraints form the relationship between the existence of a match, its heat load, and the upper bound to this heat load. When a match exists, heat load q is bounded by its upper bound, and when no match exists, q is forced

W. Verheyen, N. Zhang / Chemical Engineering Science 61 (2006) 7730 – 7753

7739

• weighed cold utility costs; • weighed hot utility costs.

to zero: q(i, j, k, p) − Qup · z(i, j, k)0, i ∈ HP, j ∈ CP, k ∈ ST, p ∈ PR,

(18)

Finally, the objective function is defined as

qcu (i, p) − Qup · zcu (i)0,

i ∈ H P , p ∈ PR,

(19)

min TAC ⎡

j ∈ CP, p ∈ PR,

(20)

= AF · ⎣

qhu (j, p) − Qup · zhu (j )0,

z(i, j, k), zcu (i), zcu (j ) ∈ {0, 1}.

dt(i, j, k, p)ti(i, k, p) − tj(j, k, p) + DTup · (1 − z(i, j, k)), i ∈ HP, j ∈ CP, k ∈ ST, p ∈ PR,

(22)

dt(i, j, k + 1, p) ti(i, k + 1, p) − tj(j, k + 1, p) + DTup · (1 − z(i, j, k)), i ∈ HP, j ∈ CP, k ∈ ST, p ∈ PR,

(23)

dtcu (i, p)ti(i, NOK+1, p)−Tcu , out+DTup · (1−zcu (i)), i ∈ HP, p ∈ PR,

(24)

dthu (j, p)Thu , out − tj(j, 1, p) + DTup · (1 − zhu (i)), j ∈ CP, p ∈ PR.

(25)

The temperature difference on each side of each heat exchanger is then limited by dt(i, j, p)DTmin .

(26)

Next, the total hot utility available can be limited by the following constraint:  qhu (j, p)HUup , j ∈ CP, p ∈ PR. (27) j ∈CP

In the model proposed by Aaltola (2002), to avoid modelling of bypasses in the MINLP model, his objective function considers the area of one match to be the mean value of areas in different periods, and therefore there is an under-estimation of the total area cost and the TAC. He reasoned that the introduction of maximum area into the formulation would introduce non-linearities into the set of constraints or non-linearities with discontinuous derivatives into the objective function. Therefore, the objective function in that paper contains the following elements: • unit costs for all heat exchangers including utility exchangers; • mean area costs for all matches (i, j, k); • mean area costs for cold utility matches; • mean area costs for hot utility matches;

Cf · z(i, j, k)

i∈H P j ∈CP k∈ST

(21)

To ensure feasible driving forces for the heat exchangers, temperature differences for each side of the heat exchanger are calculated and activated by binary variables if the match exists. When no match exists between a pair of streams at a given stage, there is an upper bound for temperature difference active. For the hot and cold utilities, similar equations exist on the variable temperature side:

  

+

 

Cf · zcu (i) +

i∈H P CU

 

⎤ Cf · zhu (j )⎦

j ∈CP H U



   1 · C N OP p∈P R i∈H P j ∈CP k∈ST  B q(i, j, k, p) · LMTD(i, j, k, p) · U (i, j )  B   1 qcu (i, p) +AF · C· · N OP LMTDcu (i, p) · U (i, cu) p∈P R i∈H P  B  1  qhu (j, p) +AF · C· · NOP LMTDhu (j, p) · U (hu, j ) + AF ·

p∈P R

j ∈CP

 DOP(p)  + Ccu · qcu (i, p) · NOP p∈P R

i∈H P

p∈P R

j ∈CP

 DOP(p)  + · Chu · qhu (j, p). NOP

(28)

3.2. Improvements to the MINLP model In the above model, the heat balances are similar to those in the original single period model and extended to multiple periods. Also the temperature feasibility constraints and logical constraints are quite easily extended from single period to multiple periods. The accuracy of the model hereby relies on the two major assumptions made to formulate this model: • the assumption of iso-thermal mixing and; • the average heat exchanger area in the objective function. Furthermore, it should be mentioned that the assumption of average area in the MINLP model can lead to sub-optimal solutions. The following results are copied from the report of one of the optimisations performed with the average area MINLP model (before the NLP improvement stage). Fig. 10 shows that there are two matches between streams H1 and C1, occurring at stages 1 and 2 of the superstructure. They also show that the first match exists at period 1 only, with an area of 19, 000 m2 , and the second match exists for periods 2 and 3, with a heat transfer area of 16, 600 m2 to be installed. The calculation performed by the objective function related to this area cost is as follows: TAC = · · · + AF · C · [ 13 · (16, 566.13 + 14103.026) +

1 3

· 18, 899.69],

7740

----

W. Verheyen, N. Zhang / Chemical Engineering Science 61 (2006) 7730 – 7753

has been used for LMTD calculations:

VARIABLE A.L area (in m2)

q(i, j, k, p) , LMTD(i, j, k, p) · U (i, j ) i ∈ HP, j ∈ CP, k ∈ ST, p ∈ PR,

A(i, j, k)  INDEX 1

= H1

P1

P2

P3

Ahu (j ) Acu (i) 

C1.1 C1.2

16566.130

14103.026

18899.690

Fig. 10. Results from average area MINLP stage.

while the actual cost of the installed area would be equal to TAC = · · · + AF · C · (16, 566.13 + 18, 899.69). Even though the NLP improvement model has not been applied to the solution yet, it would be generally not logical and not optimal to have different heat exchangers operating in different periods, as shown in this case. Also, these solutions will provide worse starting points to the NLP improvement model, fixing the heat exchanger units. A more thorough spreadsheet calculation of the TAC of the entire network shows an error between the true area cost and the average area cost of approximately 2.5 million ¥. It is the intention of this work to investigate the assumption of average area in the MINLP stage and propose an MINLP model which can predict HEN featuring minimum TAC cost with the calculation of maximum areas. Chen and Hung (2004) considered the issue of maximum area in their formulation by using the discontinuous function of MAX in the objective, which results in more expensive computation. The calculation of the maximum area, however, can also be achieved by allowing more non-linear equations into the constraints body. In order to account for the maximum area per period of each heat exchanger, the following changes have been made to the existing model: New variables are added into the formulation: A(i, j, k) is the maximum area of match i, j at stage k of the superstructure, Ahu (j ) is the maximum area of the hot utility matches of stream j , Acu (i) is the maximum area of the cold utility matches of stream I . The new equations added to the body of constraints limit the new variable A to be greater than or equal to each of the previously defined areas. The minimisation of the objective function will push the areas towards its minimum possible value, which is the maximum area per period. Another advantage is that, by using inequality constraints, it is easier for the solver to find a feasible solution. In these equations, Paterson approximation

qhu (j, p) , LMTDhu (j, p) · U (hu, j )

qcu (i, p) , LMTDcu (i, p) · U (i, cu)

(29) j ∈CP, p∈PR,

(30)

i ∈ HP, p ∈ PR. (31)

Finally, the objective function for the new MINLP model will be defined as the sum of all capital costs (heat exchanger unit costs and maximum heat exchanger area costs) and operating costs (hot utility costs and cold utility costs). Eq. (32) represents the new and improved objective function, including the maximum areas instead of the average areas per period min TAC ⎡

  

= AF · ⎣

Cf · z(i, j, k)

i∈H P j ∈CP k∈ST

+

 

Cf · zcu (i) +

i∈H P CU

  

+ AF ·

 



C · A(i, j, k)B

C · Ahu (j )B + AF ·

j ∈CP

Cf · zhu (j )⎦

j ∈CP H U

i∈H P j ∈CP k∈ST

+ AF ·





C · Acu (i)B

i∈H P

 DOP(p)  + Ccu · qcu (i, p) · NOP p∈P R

i∈H P

p∈P R

j ∈CP

 DOP(p)  + · Chu · qhu (j, p). NOP

(32)

3.3. Modified NLP model 3.3.1. Existing NLP model Aaltola (2002) proposed an NLP improvement model to improve the solution of the simplified multi-period MINLP model for two reasons: • to take away the assumption of isothermal mixing; • to take into account the maximum area into the objective function. The author decided to introduce the average area into the objective function to avoid any non-linear constraints or discrete elements in the objective function. In this way, his MINLP model becomes fast and robust to solve. The slack variables used in the NLP model represent the decrease of exchanger area with respect to some maximum exchanger area, which can be obtained from the results of the MINLP model. The slack variables are then minimised together with the TAC function inside the objective function.

W. Verheyen, N. Zhang / Chemical Engineering Science 61 (2006) 7730 – 7753

Most constraints of the NLP model can be easily taken and rephrased from the MINLP model. The balances and constraints governing this model are described next: The two overall heat balance equations, (1) and (2), define the heat load required by each stream and are therefore unaltered with respect to the MINLP model: [Tiin (i, p) − Tiout (i, p)] · F i(i, p)   q(i, j, k, p) + qcu (i, p), =

i ∈ HP, p ∈ PR,

k∈ST j ∈CP

(33) [Tjin (j, p) − Tjout (j, p)] · Fj(j, p)   q(i, j, k, p) + qhu (j, p), =

j ∈ CP, p ∈ PR.

k∈ST i∈H P

(34) Similarly, the stage heat balances define the stage temperatures for all streams and are therefore also unaltered:  [ti(i, k, p) − ti(i, k + 1, p)] · Fi(i, p) = q(i, j, k, p), j ∈CP

k ∈ ST, i ∈ HP, p ∈ PR,

(35)

[tj(j, k, p) − tj(j, k + 1, p)] · Fj(j, p)  = q(i, j, k, p), k ∈ ST, j ∈ CP, p ∈ PR.

(36)

The energy balances for the utility matches are independent on any changes of variables inside the superstructure and are therefore unaltered: [ti(i, NOK + 1, p) − Tiout (i, p)] · Fi(i, p) = qcu (i, p), (37)

[Tjout (j, p) − tj(j, 1, p)] · Fj(j, p) = qhu (j, p), j ∈ CP, p ∈ PR.

ti(i, k + 1, p) =

 1 · fh (i, j, k, p) · ths (i, j, k, p), Fi(i, p) j ∈CP

i ∈ HP, k ∈ ST, p ∈ PR, tj(j, k, p) =

(43)

 1 fc (i, j, k, p) · tcs (i, j, k, p), · Fj(j, p) i∈H P

j ∈ CP, k ∈ ST, p ∈ PR.

(44)

If a stream split does not exist, the stage outlet temperatures and branch flow rates are assigned directly: ti(i, k + 1, p) = ths (i, j, k, p), i ∈ HP, j ∈ CP, k ∈ ST, p ∈ PR,

(45)

tj(j, k, p) = tcs (i, j, k, p), i ∈ HP, j ∈ CP, k ∈ ST, p ∈ PR,

(46)

fh (i, j, k, p) = Fi(i, p), i ∈ HP, j ∈ CP, k ∈ ST, p ∈ PR,

i ∈ HP, j ∈ CP, k ∈ ST, p ∈ PR.

(47)

(38)

j ∈CP

(48)

The assignment of stream inlet and outlet temperatures is independent of any changes of variables within the superstructure, and remains unaltered. Temperature feasibility constraints to maintain the monotonic increase or decrease in stage temperature are also unchanged with respect to the MINLP model: Tiin (i, p) = ti(i, 1, p),

i ∈ HP, p ∈ PR,

Tjin (j, p) = tj(j, NOK + 1, p),

When parallel heat exchangers exist, i.e., when there is a split, then the following mass balance equations define the flow rates of the split streams:  fh (i, j, k, p)=Fi(i, p), i∈HP, k∈ST, p∈PR, (39) 

When a stream split exists, the final interval exit temperatures can be calculated by the following mixing equations:

fc (i, j, k, p) = Fj(j, p),

i∈H P

i ∈ HP, p ∈ PR

7741

ti(i, k, p) ti(i, k+1, p), tj(j, k, p) tj(j, k+1, p),

k∈ST, i∈HP, p∈PR, k∈ST, j ∈CP, p∈PR,

Tiout (i, p) ti(i, NOK + 1, p), Tjout (j, p)tj(j, 1, p),

j ∈ CP, p ∈ PR,

i ∈ HP, p ∈ PR,

j ∈ CP, p ∈ PR.

(49) (50) (51) (52) (53) (54)

i∈H P

Where a stream split exists, the following temperature feasibility constraints make sure that the heat exchanged between the streams is positive:

At the stages where a stream split exists, it is necessary to define additional energy balances at the level of the individual heat exchanger in order to define their outlet temperatures:

ti(i, k, p) ths (i, j, k, p),

fc (i, j, k, p)=Fj(j, p),

j ∈CP, k∈ST, p∈PR.

(40)

fh (i, j, k, p) · [ti(i, k, p) − ths (i, j, k, p)] = q(i, j, k, p), i ∈ H P , j ∈ CP, k ∈ ST, p ∈ PR,

(55)

tcs (i, j, k, p)tj(j + 1, k, p), (41)

fc (i, j, k, p) · [tcs (i, j, k, p) − tj(j, k + 1, p)] = q(i, j, k, p) i ∈ HP, j ∈ CP, k ∈ ST, p ∈ PR.

i ∈ HP, j ∈ CP, k ∈ ST, p ∈ PR,

(42)

i ∈ HP, j ∈ CP, k ∈ ST, p ∈ PR.

(56)

The three logical constraints that assign the relationship between the existence of a match and its corresponding heat load

7742

W. Verheyen, N. Zhang / Chemical Engineering Science 61 (2006) 7730 – 7753

reappear also: q(i, j, k, p) − Qup · z1 (i, j, k)0, i ∈ HP, j ∈ CP, k ∈ ST, p ∈ PR,

(57)

qcu (i, p) − Qup · zcu1 (i)0,

(58)

qhu (j, p) − Qup · zhu1 (j )0,

i ∈ HP, p ∈ PR, j ∈ CP, p ∈ PR,

z1 (i, j, k), zcu1 (i), zcu1 (j ) ∈ {0, 1}.

(59) (60)

When a match exists, it is necessary to ensure feasible driving forces for that heat exchanger. When there is no split of a certain hot or cold stream, the driving force will be calculated from the respective stage temperatures and constraints (61), (62) are activated:

In these equations, MAXarea is a parameter which represents the maximum area obtained from the MINLP model results. Sometimes an increase of 1–5% in addition to the MINLP result is favourable. The formulation of the constraints for the slack variables makes sure that the objective function minimises the TAC while taking into account the maximum of the areas over all periods of operation. However, the introduction of constraint (68) contains non-linearities inside the LMTD calculation. Therefore, there is an inconsistency in his pursue for obtaining linear constraints by using average area in the first stage, but allowing non-linear constraints in the second stage. The objective function is defined as follows:  DOP(p)  minobj = Ccu · qcu (i, p) · NOP p∈P R

dt(i, j, k, p)=ti(i, k, p)−tj(j, k, p)+DTup · (1−z1 (i, j, k)), i ∈ HP, j ∈ CP, k ∈ ST, p ∈ PR, (61) dt(i, j, k + 1, p) = ti(i, k + 1, p) − tj(j, k + 1, p) + DTup · (1 − z1 (i, j, k)), i ∈ HP, j ∈ CP, k ∈ ST, p ∈ PR. (62) When there is a stream split of any hot or cold stream, the driving force is calculated from the heat exchanger inlet and exit temperatures, and constraints (63) and (64) become active: dt(i, j, k, p) = ti(i, k, p) − tcs (i, j, k, p) + DTup · (1 − z1 (i, j, k)), i ∈ HP, j ∈ CP, k ∈ ST, p ∈ PR, (63) dt(i, j, k + 1, p) = ths (i, j, k, p) − tj(j, k + 1, p) + DTup · (1 − z1 (i, j, k)), i ∈ HP, j ∈ CP, k ∈ ST, p ∈ PR. (64) Also, the driving forces for the hot utility and cold utility matches can be defined and constrained in a way similar to the MINLP model: dtcu (i, p) ti(i, NOK + 1, p) − Tcu , out + DTup · (1 − zcu1 (i)), i ∈ HP, p ∈ PR, (65) dthu (j, p) Thu , out − tj(j, 1, p) + DTup · (1 − zhu1 (i)), j ∈ CP, p ∈ PR.

(66)

Finally, the minimum approach temperature can be defined by dt(i, j, k, p)Dmin , i ∈ HP, j ∈ CP, k ∈ ST, p ∈ PR.

(67)

According to the NLP formulation by Aaltola (2002), these final constraints (68) and (69) define the slack variables which are related to the decrease in area: q(i, j, k, p) = MAXarea (i, j, k)·LMTD(i, j, k, p)·U (i, j )−s(i, j, k, p), i ∈ HP, j ∈ CP, k ∈ ST, p ∈ PR, (68) smin (i, j, k)s(i, j, k, p), i ∈ HP, j ∈ CP, k ∈ ST, p ∈ PR.

(69)

i∈H P

 DOP(p)  + · Chu · qhu (j, p) NOP p∈P R j ∈CP    − Cw · smin (i, j, k).

(70)

i∈H P j ∈CP k∈ST

In this model, the solution is achieved by an iterative procedure while altering the weighed parameter and re-running the model until the minimum TAC is achieved. As can be noticed, it is not possible to directly read the objective value of TAC from these optimisation results. Since the optimisation program is able to provide the maximum area per heat exchanger as its output, the bypass fractions can be calculated by a straightforward procedure, described in more detail by Aaltola (2002). 3.3.2. Improvements to the NLP model The use of slack variables and a weighed parameter into the NLP formulation are necessary to account for the maximum area per heat exchanger and still obtain a set of linear constraints. When non-linearities are allowed into the set of constraints, it is not necessary any more to introduce these additional variables and this weighed parameter. This section discusses changes which are made to the previous NLP model to provide an overall MINLP–NLP model which could perform better than the existing model. The new variables added into the formulation are the same as in the MINLP model: A(i, j, k) is the maximum area of match (i, j ) at stage k of the superstructure, Ahu (j ) the maximum area of the hot utility matches of stream j , Acu (i) the maximum area of the cold utility matches of stream I . Again, the minimisation of the objective function will automatically push the areas towards its minimum possible value, which is the maximum area per period: q(i, j, k, p) , LMTD(i, j, k, p) · U (i, j ) i ∈ HP, j ∈ CP, k ∈ ST, p ∈ PR,

A(i, j, k)

Ahu (j ) Acu (i) 

qhu (j, p) , LMTDhu (j, p) · U (hu, j )

qcu (i, p) , LMTDcu (i, p) · U (i, cu)

(71) j ∈CP, p∈PR,

(72)

i ∈ HP, p ∈ PR. (73)

W. Verheyen, N. Zhang / Chemical Engineering Science 61 (2006) 7730 – 7753

Finally, the objective function for the new NLP model will be defined as

i∈H P

Next, the multi-period MINLP–NLP models are compared to each other using the multi-period stream data of this case study. Finally, the new proposed model has been evaluated in terms of computer solving time and its capability of handling problems of increased size in terms of number of streams, number of stages, and number of periods.

 DOP(p)  + · Chu · qhu (j, p) NOP p∈P R



+ AF · ⎝

j ∈CP

  

Cf · z1 (i, j, k)

i∈H P j ∈CP k∈ST

+



Cf · zhu1 (j ) +

j ∈CP



+ AF · ⎝



+

j ∈CP

4. A case study

Cf · zcu1 (i)⎠

C · A(i, j, k)B

i∈H P j ∈CP k∈ST





i∈H P

  

C · Ahu (j )B +



which is formulated as • a linear objective function with non-linear constraints, or • a non-linear objective function with linear set of constraints.

 DOP(p)  minobj = Ccu · qcu (i, p) · NOP p∈P R

7743

⎞ C · Acu (i)B ⎠ .

(74)

i∈H P

3.4. Comments Aaltola (2002) described that the main goal in his thesis is to provide a simplified model which can give accurate and good solutions for large-scale industrial processes. To be able to achieve this goal, a set of simplifications are made. Firstly, the simplified superstructure of Yee et al. (1990) has been applied and extended to multiple periods. This superstructure has shown its ability to provide good results in many cases. The assumption of isothermal mixing has been overcome by the introduction of an NLP improvement model. By extending the single period model to multiple periods, other assumptions need to be made in order to preserve the linearity of the set of constraints. This has been done by taking into account the average area into the objective function. In this way, non-linear constraints and nondiscrete terms in the objective function can be avoided. However, this simplification has an effect on the obtained solution as well as on the formulation of the NLP improvement model. In this second model, the author makes use of slack variables to account for the decrease in area. Also, in the objective function, a weighed parameter needs to be used to maximise this area reduction. From Eq. (61) it can be seen that, in fact, there are non-linear constraints present in the NLP model, in the form of LMTD functions. This justifies the reason of changing the objective function in the NLP model, shown in Eq. (67). This model contains no slack variables, no weighed parameters (which mean less user interaction and less time consuming) and has non-linear constraints to the same degree. Also, in the new NLP model, the true objective function is directly minimised, hereby confirming the optimality of the solution. In the following cast study, firstly, a single period case study has been used to compare the results of the MINLP model,

The case study presented in this work is the HEN for the vacuum gas oil (VGO) hydrotreater unit of an oil refinery. The VGO hydrotreating unit takes VGO and treats and upgrades it to more useful products. During this hydrotreating operation, the boiling range of the feedstock is not altered significantly. Normal hydrotreater reaction temperatures are around 300.400 ◦ C and the reactions are catalytic. Moreover, the process of hydrotreating consumes large amounts of hydrogen. In the hydrotreating operation, the feed is first mixed with the hydrogen and heated to the reactor inlet temperature. This can be done by a single furnace but is more efficiently done by reusing waste process heat first and using the furnace to provide the high temperature heat. Although many reactions take place inside the hydrotreating units, the net reaction is of exothermic nature and the reactor outlet temperature tends to be somewhat higher than the inlet temperature. Then, the reactor effluent is cooled down and passed to a series of flash separators to separate and recycle hydrogen for reuse. Then the products are passed through a hydrogen sulphide stripper to remove most hydrogen sulphide present. Finally a series of distillation columns (either simple or complex) separate light hydrocarbons, gasoline, diesel and gas oil. During the duration of the hydrotreating process, the catalyst used in the reactor gradually loses its activity. Therefore, to compensate the loss of reaction rate due to loss of catalyst activity, the reactor inlet temperature is gradually increased. This causes not only disturbances to the HEN in terms of stream inlet and outlet temperatures, but also in terms of reactor outlet compositions, which are the result of changes in operating conditions. This change in composition then affects the flow rates of diesel and gasoline formed, as well as the amount of hydrogen used and recycled. These changes in operating conditions give rise to the need for designing multiple period HENs which are not only feasible for these different conditions, but can provide optimum solutions for all operating conditions. HYSYS simulation software has been used to model the steady-state performance of the VGO hydrotreater unit and subsequent separation units. The operating conditions and some unknown parameters have been obtained from the specification of the distillate products in terms of boiling point curves. From the simplified flow sheet, the necessary hot streams (streams which need to be cooled down) and cold streams (streams which

7744

W. Verheyen, N. Zhang / Chemical Engineering Science 61 (2006) 7730 – 7753

Table 1 Stream identification

Table 3 Extracted stream data for MOR

Stream name

Description of stream in flow sheet

Stream

Inlet temperature (◦ C)

Outlet temperature (◦ C)

H1 H2

Reactor outlet stream to be cooled down Diesel stream from distillation side stripper to be cooled down Gas oil stream from reboiler to be cooled down Reactor inlet stream to be heated to reaction temperature Stream from low pressure separator heated to H2S stripper Stream from H2S stripper heated to distillation inlet temperature Side stripper reboiler

Heat capacity flow rate F (kW/K)

H1 H2 H3 C1 C2 C3 C4

406 160 362 72 62 220 250

60 40 60 365 210 370 290

205.0 198.8 136.4 210.3 141.0 175.4 318.7

H3 C1 C2 C3 C4

Table 4 Extracted stream data for EOR

Table 2 Extracted stream data for SOR Stream

Inlet temperature (◦ C)

Outlet temperature (◦ C)

Heat capacity flow rate F (kW/K)

H1 H2 H3 C1 C2 C3 C4

393 160 354 72 62 220 253

60 40 60 356 210 370 284

201.6 185.1 137.4 209.4 141.6 176.4 294.4

Stream

Inlet temperature (◦ C)

Outlet temperature (◦ C)

Heat capacity flow rate F (kW/K)

H1 H2 H3 C1 C2 C3 C4

420 160 360 72 62 220 249

60 40 60 373 210 370 286

208.5 175.2 134.1 211.1 140.5 174.5 271.2

need to be heated up) have been identified, and the necessary stream data have been collected from the simulation. This simulation have been repeated for three operating periods:

as input to the MINLP–NLP model to allow a robust and reasonably fast solution. There are no data present of any occasional scenarios, which should be feasible, and therefore the HEN design does not involve the feasibility test.

• start-of-run (SOR), • end-of-run (EOR), • mid-of-run (MOR), which uses the average values of SOR and EOR.

4.1. Single period evaluation and comparison

These three periods are assumed to have equal durations. In the model it is possible to specify the duration of each period by adjusting parameter DOP(p). The data extracted from the flowsheet are: • stream inlet and outlet temperatures; • stream specific heats at heat exchanger inlet and outlet; • stream mass flow rates. The identified streams are presented in Table 1. In this case, the distillation column uses steam stripping in the main column and therefore does not need a reboiler at this position. At higher distillation temperatures, steam stripping proves to be more beneficial than the use of a reboiler in many cases. The necessary data obtained from the HYSYS simulation are combined in Tables 2–4 representing periods SOR, MOR, and EOR, respectively. The heat capacity flow rates shown in this table are the average heat capacity flow rates of each stream. It is possible to obtain linear or higher level polynomial functions of Cp with respect to temperature. Constant Cp values have been chosen

This section discusses the comparison of the single period MINLP model proposed by Yee and Grossmann (1990), extended with the NLP model for allowing non-isothermal mixing, in two different formulations: • A non-linear objective function with a set of linear constraints, which will be the basis of the average area multiperiod model proposed by Aaltola (2002). • A linear objective function (in this case the area cost exponent is unity) with non-linear constraints, which will be the bases of the maximum area multi-period model, proposed in this work. In this comparison, both models are run a number of times, with different values of minimum allowable temperature difference (DTmin ) and different values of maximum hot utility available (HUup ). Although in most cases these constraints are artificial and are not active, they still influence the solving path and will lead the program towards a set of different local minima. Then, from this set, the best local minimum is chosen and is presented as the solution to the heat exchanger design problem for this case. The ranges that are investigated are: • DTmin from 10 ◦ C up to 35 ◦ C, • HUup from 50 MW down to 10 MW.

W. Verheyen, N. Zhang / Chemical Engineering Science 61 (2006) 7730 – 7753

7745

Table 5 Results from single period comparison Single period comparison DTmin (◦ C)

Model > HUup

New model Objective function

Model by Yee and Grossmann Objective function

DTmin (◦ C)

Model > HUup

New model Objective function

Model by Yee and Grossmann Objective function

10

50,000 45,000 40,000 35,000 30,000 25,000 20,000 15,000 10,000

6,684,732 6,536,714 7,897,972 7,453,795 7,775,521 6,407,668 7,158,953 INFES INFES

7,444,592 7,446,259 7,444,592 7,452,128 7,773,855 6,407,668 6,997,515 INFES INFES

25

50,000 45,000 40,000 35,000 30,000 25,000 20,000 16,000 10,000

6,359,872 6,304,400 6,303,146 6,306,067 6,303,146 6,363,657 7,098,245 INFES INFES

6,328,205 6,324,872 6,417,735 6,324,872 6,306,067 6,407,668 6,407,668 INFES INFES

15

50,000 45,000 40,000 35,000 30,000 25,000 20,000 15,000 10,000

6,355,915 6,539,951 6,324,872 7,453,795 6,545,980 6,407,668 7,158,953 INFES INFES

6,355,915 6,255,495 6,992,455 7,452,128 7,776,062 6,407,668 7,157,286 INFES INFES

30

50,000 45,000 40,000 35,000 30,000 25,000 20,000 18,000 10,000

6,324,872 6,326,538 6,328,205 6,328,205 6,306,067 6,357,009 INFES INFES INFES

6,253,829 6,324,872 6,326,538 6,324,872 6,324,872 6,406,002 INFES INFES INFES

20

50,000 45,000 40,000 35,000 30,000 25,000 20,000 15,000 10,000

6,324,872 6,417,735 6,329,871 7,452,128 6,326,538 6,407,668 6,982,160 INFES INFES

6,329,871 6,329,871 7,444,592 7,453,795 6,326,538 6,407,668 7,157,286 INFES INFES

35

50,000 45,000 40,000 35,000 30,000 25,000 20,000 18,000 10,000

6,312,304 6,306,067 6,304,400 6,306,067 6,306,067 6,363,657 INFES INFES INFES

6,324,872 6,324,872 6,324,872 6,324,872 6,324,872 6,361,991 INFES INFES INFES

It should be noted that, for larger values of DTmin , the composite curves are horizontally shifted away from each other and this increases the minimum hot and cold utilities that are needed. Therefore, no feasible networks will exist for the lowest values of HUup at increased DTmin . The ranges of HUup and DTmin have been carefully chosen based on a number of observations. The lower limit to DTmin can in theory be reduced to zero, since it will not become an active constraint at these low values. When the driving force (temperature difference) becomes very small, the area increases drastically and the result would be uneconomical. In practice, driving forces lower than 10 ◦ C are not advised for hydrocarbon oil heat exchangers because this will cause excessive fouling of the heat exchanger unit. On the other hand, it is been determined that, if a lower bound on allowable driving force is set to be 40 ◦ C or greater, the values of objective function obtained from the model are becoming less economical. Similarly, for HUup , the smallest value for HUup is set by the feasibility limit and the targets for minimum possible hot utility duty. Furthermore, for values of HUup greater than 50 MW, there has been no solution reported which approaches the best local minimum obtained. At the low values of DTmin , the formulation with linear constraints tends to calculate excessively large derivatives, while the formulation with the linear objective function has no problems handling this parameter. To obtain a valid comparison between two

models, both need to be based on the same cost functions and the same economic parameters. For simplicity, the cost functions and economic parameters used in this work have been adopted from Aaltola (2002) and are stated next. The cost of a heat exchanger has been evaluated as follows: Heat exchanger cost = Cf + C · (Area)B , where Cf is the fixed cost of a heat exchanger = 8333.3 ¥, C is area cost coefficient for a heat exchanger = 641.7 ¥/m(2B) , B is area exponent = 1. Furthermore, one needs to specify the following economic criteria: Ccu is the cost of cold utility = 1.3 ¥/kWy, Chu is cost of hot utility = 115.2 ¥/kWy, annualisation factor used = 0.2. Table 5 shows a list of optimal solutions reached by both formulations (combined MINLP–NLP) for a range of different DTmin and HUup . Table 5 shows that results obtained by the formulation with linear objective function and non-linear constraints are of similar magnitude, compared to the model which adds all nonlinearities into the objective function and keeps a linear set of constraints. Local minima vary from 6.25 million ¥ up to 7.77 million ¥ or more. To compare the best results obtained by both formulations, the DTmin is limited to values of 10 ◦ C or more for physical feasibility. As was discussed previously, this lower bound makes sure that there is no excessive fouling occurring in the

7746

W. Verheyen, N. Zhang / Chemical Engineering Science 61 (2006) 7730 – 7753 Stage 1

Stage 2 4.3 MW

15.9 MW

Table 6 Multi-period simultaneous model comparison

H1 3.8 MW

18.4 MW

Multiperiod comparison

7.1 MW

DTmin

Model-> HUup

New model Objective function

Aaltola’s model Objective function

20

50,000 45,000 40,000 35,000 30,000 25,000 20,000 17,000 10,000

6,439,285 6,418,713 6,435,752 6,491,251 6,435,752 7,625,869 7,295,796 INFES INFES

6,439,481 8,157,217 6,422,375 6,439,481 6,422,375 6,584,265 7,731,043 INFES INFES

25

50,000 45,000 40,000 35,000 30,000 25,000 20,000 18,000 10,000

6,416,403 7,591,377 6,431,203 6,410,355 6,431,203 6,641,023 INFES INFES INFES

6,424,042 6,353,888 6,439,481 6,424,042 6,421,253 6,578,335 INFES INFES INFES

30

50,000 45,000 40,000 35,000 30,000 25,000 20,000 18,000 10,000

6,436,234 6,435,856 6,436,234 6,320,508 6,437,870 6,644,294 INFES INFES INFES

6,441,148 6,439,481 6,364,146 6,439,481 6,422,375 6,584,265 INFES INFES INFES

H2 H3 C1 4.4 MW

46.9 MW C2

14.0 MW

21.0 MW C3 12.4 MW C4

9.1 MW

Fig. 11. Solution for the new model at HUup = 30 MW and DTmin = 25 ◦ C.

Stage 1

Stage 2 11.9MW

H1 22.2 MW H2 7.1 MW H3 C1 13.4 MW

46.1 MW C2 21.0 MW C3

14.1 MW

12.4 MW C4 9.1 MW

Fig. 12. Solution for existing model at HUup = 50 MW and DTmin = 30 ◦ C.

heat exchangers. The best local optimum is reached by the new formulation with non-linear constraints for a DTmin boundary of 25 ◦ C and upper hot utility bound of 30 MW. The best local optimum which is reached from the formulation of Yee and Grossmann (1990) with linear constraints was obtained at values of DTmin and HUup of 30 ◦ C and 50 MW, respectively. The final minimum cost network obtained for a single period by both formulations is shown in Figs. 11 and 12. The first network (Fig. 11) features 11 heat exchangers, with a hot utility consumption of 27.5 MW and a cold utility consumption of 41.4 MW. These values can be compared to the minimum hot and cold utility consumption from pinch analysis at the same DTmin of 25 ◦ C: minimum hot utility consumption: 17.1 MW minimum cold utility consumption: 30.8 MW. This network has a TAC of 6.30 million ¥ and has a reasonable level of energy recovery. The second network (Fig. 12) features nine heat exchangers, with hot utility consumption of 27.5 MW and cold utility consumption of 41.2 MW. These values are also compared to the minimum hot and cold utility consumption from pinch analysis at a DTmin of 30 ◦ C: minimum hot utility consumption: 18.9 MW minimum cold utility consumption: 32.6 MW. The network has a TAC of approximately 6.25 million ¥ and features the same level of energy recovery with a larger minimum approach temperature.

To conclude, it is shown that, for this scale of problem, both formulations predict similar results. For both formulations, the global optimum cannot be guaranteed due to non-convexities in the model, and a set of local optima is produced, based on the adjustment of two control parameters, which are HUup and DTmin . From this set of local optima, it is possible to select the most cost effective and controllable or flexible network structure. 4.2. Multi-period model comparison In this section, the new model which evaluates maximum areas in the MINLP phase is compared with the multi-period simultaneous model, proposed by Aaltola (2002), which takes into account average areas in the MINLP phase. For the MINLP phase, both models use CONOPT as the default NLP solver and CPLEX for the relaxed MIP problem. Table 6 shows the results and comparison of both combined MINLP–NLP models, for minimum temperature differences ranging from 20 up to 30 ◦ C. In this table, the term “INFES” is used when it is impossible to achieve a feasible network featuring those upper limits to hot utility for the given DTmin . These observations are confirmed by using the pinch analysis to show that the minimum hot utility target is violated.

W. Verheyen, N. Zhang / Chemical Engineering Science 61 (2006) 7730 – 7753

Similar to the single period comparison, DTmin has been chosen to vary from 20 up to 30 ◦ C for the following reasons: • Some heat exchangers are able to handle minimum temperature differences which approximate zero. The choice of minimum temperature difference depends mainly on the fluid used, on the application, and on the heat exchanger type. For general hydrocarbon oil heat exchangers, temperature differences below 10 ◦ C on either side of the heat exchanger cause early fouling of the exchanger. Many hydrocarbon oil exchangers prefer a temperature difference of 20 ◦ C or higher. The lowest minimum has been obtained in this range of temperatures. • The maximum value of 30 ◦ C has been chosen because, at higher bounds on temperature difference, there is an increase in annualised cost due to the increase in hot and cold utilities used. The optimal network for the new model is obtained at a temperature difference bound of 30 ◦ C and hot utility bound of 35 MW and has a total annual cost of approximately 6.32 million ¥. The optimal network for Aaltola’s multi-period model is obtained at a temperature difference bound of 25 ◦ C and a hot utility bound of 45 MW and features a total annual cost of more than 6.35 million ¥. Furthermore it can be observed from the above table that, in overall, the average annual network cost obtained from the new model is lower than the average annual network cost obtained by the previous model. 4.3. Model results The optimal results from both models, after the MINLP stage as well as after the NLP stage, are reported in detail in this section. 4.3.1. Results from MINLP stage (Aaltola’s model) Tables 7 and 8 show the results from the MINLP stage of Aaltola’s model at a DTmin of 25 ◦ C and HUup at 45 MW for all three periods, in terms of: • Amax , the maximum area of the heat exchanger for the three periods; • A/Amax , the ratio of used area over maximum area of the heat exchanger; • fh , the flow rate of the hot stream through the heat exchanger; • fc , the flow rate of the cold stream through the heat exchanger; • th (k), the hot side heat exchanger inlet temperature of the hot stream; • th (k + 1), the cold side stage temperature of the hot stream, which also equals to the exchanger outlet temperature of the hot stream (by isothermal mixing); • tc (k), the hot side stage temperature of the cold stream, which also equals to the exchanger outlet temperature of the cold stream;

7747

Table 7 Results from MINLP stage (Aaltola’s model) Match (i, j, k)

1,1,1

1,1,2

1,4,1

3,2,2

3,3,1

p=1 Amax A/Amax fh fc th (k) th (k + 1) tc (k) tc (k + 1) Q

m2 % kW/K kW/K ◦C ◦C ◦C ◦C kW

2333 53.0 100.6 209.4 393.0 302.7 293.6 250.2 9088

7701 99.0 201.6 209.4 302.7 117.6 250.2 72.0 37,313

1623 74.5 101.0 294.4 393.0 302.7 284.0 253.0 9126

4058 100.0 137.4 141.6 263.9 111.4 210.0 62.0 20,951

2440 95.3 137.4 176.4 354.0 263.9 290.2 220.0 12,383

p=2 Amax A/Amax fh fc th (k) th (k + 1) tc (k) tc (k + 1) Q

m2 % kW/K kW/K ◦C ◦C ◦C ◦C kW

2333 49.0 84.5 210.3 406.0 300.2 291.9 249.4 8946

7701 100.0 205.0 210.3 300.2 118.2 249.4 72.0 37,313

1623 100.0 120.5 318.7 406.0 300.2 290.0 250.0 12,750

4058 99.2 136.4 141.0 264.4 111.4 210.0 62.0 20,875

2440 100.0 136.4 175.4 362.0 264.4 295.9 220.0 13,311

p=3 Amax A/Amax fh fc th (k) th (k + 1) tc (k) tc (k + 1) Q

m2 % kW/K kW/K ◦C ◦C ◦C ◦C kW

2333 100.0 130.8 211.1 420.0 290.9 320.3 240.4 16,882

7701 93.4 208.5 211.1 290.9 120.4 240.4 72.0 35,547

1623 78.0 77.7 271.2 420.0 290.9 286.0 249.0 10,033

4058 99.7 134.1 140.5 265.0 109.9 210.0 62.0 20,795

2440 94.5 134.1 174.5 360.0 265.0 293.0 220.0 12,741

Table 8 Utility related results from MINLP stage (Aaltola’s model) H1

H2

H3

C1

C3

Amax

m2

1523

2863

896

1051

975

P1: A/Amax Q1

% kW

93.5 11,603

93.1 22,207

100.0 7067

82.9 13,079

100.0 14,080

P2: A/Amax Q2

% kW

95.9 11,930

100.0 23,855

99.2 7012

100.0 15,371

94.1 12,996

P3: A/Amax Q3

% kW

100.0 12,584

88.1 21,020

95.5 6694

82.5 11,124

96.3 13,429

• tc (k + 1), the cold side heat exchanger inlet temperature of the cold stream; • Q, the heat duty of the heat exchanger. From the MINLP model, the structure of the HEN is fixed and is shown in Fig. 13. The NLP improvement model will then optimise the distributions of the flow rates into the parallel heat exchangers, if there exist any.

7748

W. Verheyen, N. Zhang / Chemical Engineering Science 61 (2006) 7730 – 7753 Stage 1

Stage 2

Cold Utility

Table 9 Results from NLP improvement stage (Aaltola’s model)

H1

Match (i, j, k)

1,1,1

1,1,2

1,4,1

3,2,2

3,3,1

p=1 Amax A/Amax fh fc th (in) th (out) tc (in) tc (out) Q

m2 % kW/K kW/K ◦C ◦C ◦C ◦C kW

2333 100.0 122.2 209.4 393.0 269.9 227.6 299.5 15,052

7701 100.0 201.6 209.4 273.1 111.4 72.0 227.6 32,597

1623 98.4 79.4 294.4 393.0 278.0 253.0 284.0 9126

4058 100.0 137.4 141.6 263.9 111.4 62.0 210.0 20,951

2440 95.3 137.4 176.4 354.0 263.9 220.0 290.2 12,383

p=2 Amax A/Amax fh fc th (in) th (out) tc (in) tc (out) Q

m2 % kW/K kW/K ◦C ◦C ◦C ◦C kW

2333 90.2 84.5 210.3 406.0 262.2 237.2 295.0 12,155

7701 100.0 205.0 210.3 284.5 115.0 72.0 237.2 34,754

1623 100.0 120.5 318.7 406.0 300.2 250.0 290.0 12,750

4058 99.2 136.4 141.0 264.4 111.4 62.0 210.0 20,875

2440 100.0 136.4 175.4 362.0 264.4 220.0 295.9 13,311

p=3 Amax A/Amax fh fc th (in) th (out) tc (in) tc (out) Q

m2 % kW/K kW/K ◦C ◦C ◦C ◦C kW

2333 100.0 139.7 211.1 420.0 296.6 241.5 323.2 17,251

7701 100.0 208.5 211.1 289.1 117.4 72.0 241.5 35,796

1623 95.0 68.7 271.2 420.0 274.0 249.0 286.0 10,033

4058 100.0 134.1 140.5 264.9 109.8 62.0 210.0 20,795

2440 94.8 134.1 174.5 360.0 264.9 220.0 293.1 12,758

H2 H3 C1 1,1,1

1,1,2 C2 3,2,2 C3 3,3,1 C4

Hot Utility

1,4,1

Fig. 13. Final structure after MINLP stage (Aaltola’s model).

The structure features 10 heat exchangers, amongst which five utility heat exchangers and five process-to-process heat exchangers, with only one stream split in the first stage. The areas of each heat exchanger, the stage temperatures, the heat loads and the stream flow rates through each exchanger are reported in Tables 7 and 8. This network configuration is optimal for Aaltola’s MINLP–NLP model. Since it features only one stream splitting, the NLP improvement model will rather improve the objective function in terms of more accurate prediction of network cost (based on true maximum area), rather than improve the objective function to allow non-isothermal mixing. 4.3.2. Results from NLP improvement stage (Aaltola’s model) By using the above structure as input into Aaltola’s NLP improvement model, the following tables report the finalised HEN with optimal distribution of flow rates, featuring minimum TAC. Table 9 contains the maximum areas, the area fractions, the heat duties, the flows and temperatures for each process-toprocess heat exchanger in the superstructure. Next, Table 10 shows the areas, the area fractions and the heat duties of the utility heat exchangers. Previously the assumption of iso-thermal mixing allowed the heat exchanger inlet and outlet temperatures to be represented by the stage temperatures, th (k), th (k + 1), tc (k) and tc (k + 1). When non-isothermal mixing is allowed, it is necessary to assign different temperature variables. In the tables, th (in) and th (out) represent the inlet and outlet temperatures of the hot stream into the heat exchanger units. Similarly, tc (in) and tc (out) represent the inlet and outlet temperatures of the cold stream into the heat exchanger unit. Table 10 shows that there are now seven utility heat exchangers. This can happen if the MINLP stage prints out a solution where one match exists with zero area. It is possible to overcome this by adding the following line into the MINLP stage: z(i, j, k)

 1 · A(i, j, k, p). NOP

(75)

p∈P R

This will force the binary variable representing the match (z) to become zero whenever the average area equals zero (or less

Table 10 Utility related results from NLP improvement stage (Aaltola’s model) H1

H2

H3

C1

C3

Amax

m2

1470

2863

896

1017

975

P1: A/Amax Q1

% kW

89.3 10,356

93.1 22206.6

100.0 7067

79.0 11,832

100.0 14,080

P2: A/Amax Q2

% kW

95.4 11,281

100.0 23,855

99.2 7012

100.0 14,721

94.1 12,996

P3: A/Amax Q3

% kW

100.0 11,965

88.1 21,020

95.3 6676

81.4 10,506

96.2 13,411

than one) as well. If one, however, would like to investigate the effect of allowing this match to exist, it is not necessary to add this constraint into the MINLP formulation. The updated HEN structure showing heat transfer areas for all matches is shown in Fig. 14. The final network obtained by Aaltola’s multi-period model developed by Aaltola (2002) features a minimum TAC of 6.35 million. The multi-period solution is slightly more expensive than the single period solution obtained previously

W. Verheyen, N. Zhang / Chemical Engineering Science 61 (2006) 7730 – 7753 Stage 1

Stage 2 1,470 m2

7749

Table 12 Utility related results from MINLP stage (new model)

H1

H1

H2

H3

C1

C3

C4

Amax

m2

1471

2863

1070

900

1026

194

P1: A/Amax Q1

% kW

90.8 10,602

93.1 22,207

100.0 9100

78.4 10,099

100.0 15,085

79.1 3007

P2: A/Amax Q2

% kW

94.2 11,094

100.0 23,855

98.3 8900

92.6 11,426

96.0 14,296

100.0 3697

P3: A/Amax Q3

% kW

100.0 11,984

88.1 21,020

97.1 8812

100.0 11,661

96.4 14,410

0.0 0

2,863 m2 H2 896 m2 H3 C1 1,017 m2

2,333 m2

7,701 m2 C2 4,058 m2 C3

975 m2

2,440 m2 C4 1,623 m2

Fig. 14. Final network from NLP stage (Aaltola’s model).

Table 11 Results from MINLP stage (new model)

Stage 1

Stage 2

Cold Utility

H1

Match (i, j, k)

1,1,1

1,1,2

1,2,2

1,4,1

3,1,2

3,3,1

p=1 Amax A/Amax fh fc th (k) th (k + 1) tc (k) tc (k + 1) Q

m2 % kW/K kW/K ◦C ◦C ◦C ◦C kW

2390 100.0 146.0 209.4 393.0 283.0 307.8 231.1 16,056

2913 100.0 78.6 84.2 283.0 112.6 231.1 72.0 13,402

3428 100.0 122.9 141.6 283.0 112.6 210.0 62.0 20,951

1331 71.0 55.6 294.4 393.0 283.0 273.8 253.0 6119

4257 100.0 137.4 125.2 271.2 126.2 231.1 72.0 19,924

1900 100.0 137.4 176.4 354.0 271.2 284.5 220.0 11,378

p=2 Amax A/Amax fh fc th (k) th (k + 1) tc (k) tc (k + 1) Q

m2 % kW/K kW/K ◦C ◦C ◦C ◦C kW

2390 100.0 132.8 210.3 406.0 280.6 310.7 231.5 16,657

2913 100.0 79.6 83.1 280.6 114.1 231.5 72.0 13,261

3428 100.0 125.4 141.0 280.6 114.1 210.0 62.0 20,875

1331 100.0 72.2 318.7 406.0 280.6 278.4 250.0 9053

4257 100.0 136.4 127.2 273.9 125.2 231.5 72.0 20,285

1900 100.0 136.4 175.4 362.0 273.9 288.5 220.0 12,011

p=3 Amax A/Amax fh fc th (k) th (k + 1) tc (k) tc (k + 1) Q

m2 % kW/K kW/K ◦C ◦C ◦C ◦C kW

2390 100.0 133.6 211.1 420.0 285.9 317.8 232.9 17,914

2913 100.0 85.0 89.0 285.9 117.5 232.9 72.0 14,319

3428 93.1 123.5 140.5 285.9 117.5 210.0 62.0 20,795

1331 100.0 74.8 271.2 420.0 285.9 286.0 249.0 10,033

4257 100.0 134.1 122.2 272.3 125.7 232.9 72.0 19,659

1900 100.0 134.1 174.5 360.0 272.3 287.4 220.0 11,759

H2 H3 C1 1,1,1

1,1,2

3,1,2 C2 1,2,2 C3

3,3,1 C4 Hot Utility

1,4,1

Fig. 15. Final structure after MINLP stage (new model).

(6.25 million ¥). This seems reasonable because the other two periods (MOR and EOR) require heating to a higher temperature and therefore will also need to be cooled down more than during the first period (SOR). 4.3.3. Results from MINLP stage (new model) Tables 11 and 12 show the results from the MINLP stage of the new and improved model at a DTmin of 30 ◦ C and HUup at

35 MW for all three periods. The obtained structure features a total twelve heat exchanger units, amongst which six processto-process heat exchangers and six utility heat exchangers. Also this structure has a total of two stream splits, of which one in the first stage and one in the second stage. From the new and improved MINLP model, the structure of the HEN is fixed again, and it is shown in Fig. 15. The new NLP improvement model will then optimise the distributions of the flow rates into the parallel heat exchangers, if any stream splits exist. The structure features 12 heat exchangers, amongst which six utility heat exchangers and another six process-to-process heat exchangers, with a total of two stream splits. The areas of each heat exchanger, the stage temperatures, the heat loads and the stream flow rates through each exchanger are reported in Tables 11 and 12. This network configuration is optimal for this new MINLP–NLP model. The NLP improvement model will use this structure as its inputs and will generate a minimal TAC network allowing non-isothermal mixing at all the mixers of the stream splits. 4.3.4. Results from NLP improvement stage (new model) By using the above structure as input into the new NLP improvement model, the following tables report the finalised HEN with optimal distribution of flow rates, featuring minimum

7750

W. Verheyen, N. Zhang / Chemical Engineering Science 61 (2006) 7730 – 7753

Table 13 Results from NLP improvement stage (new model)

Stage 1

Stage 2 1,423 m2

H1

Match (i, j, k)

1,1,1

1,1,2

1,2,2

1,4,1

3,1,2

3,3,1

2,863 m2 H2 998 m2 H3

p=1 Amax A/Amax fh fc th (in) th (out) tc (in) tc (out) Q

m2 % kW/K kW/K ◦C ◦C ◦C ◦C kW

2916 100.0 120.9 209.4 393.0 248.9 214.5 297.7 17,427

2240 100.0 60.3 64.9 262.5 108.2 72.0 215.3 9298

3994 100.0 141.3 141.6 262.5 114.3 62.0 210.0 20,946

1441 100.0 80.7 294.4 393.0 283.0 253.0 283.2 8876

4146 100.0 137.4 144.5 267.5 118.0 72.0 214.2 20,547

2104 100.0 137.4 176.4 354.0 267.5 220.0 287.389 11,889

p=2 Amax A/Amax fh fc th (in) th (out) tc (in) tc (out) Q

m2 % kW/K kW/K ◦C ◦C ◦C ◦C kW

2916 100.0 124.6 210.3 406.0 257.9 220.1 307.9 18,451

2240 100.0 73.2 72.7 267.6 120.6 72.0 220.1 10,768

3994 100.0 131.8 141.0 267.6 109.2 62.0 210.0 20,875

1441 100.0 80.4 318.7 406.0 282.6 250.0 281.1 9919

4146 100.0 136.4 137.6 270.0 120.5 72.0 220.2 20,391

2104 100.0 136.4 175.4 362.0 270.0 220.0 291.6 12,549

p=3 Amax A/Amax fh fc th (in) th (out) tc (in) tc (out) Q

m2 % kW/K kW/K ◦C ◦C ◦C ◦C kW

2916 100.0 137.3 211.1 420.0 270.5 223.0 320.3 20,531

2240 100.0 85.6 81.7 273.4 130.1 72.0 222.0 12,255

3994 100.0 122.9 140.5 273.4 104.2 62.0 210.0 20,795

1441 100.0 71.2 271.2 420.0 279.0 249.0 286.0 10,033

4146 100.0 134.1 129.4 268.4 122.0 72.0 223.7 19,630

2104 100.0 134.1 174.5 360.0 268.4 220.0 290.4 12,281

C1 868 m2

2,916 m2

2,240 m2

4,146 m2 C2

1,000 m2

3,394 m2 C3 2,104 m2 C4

150 m2

1,441 m2

Fig. 16. Final structure after NLP improvement stage (new model).

Table 14 Utility related results from NLP improvement stage (new model) H1

H2

H3

C1

C3

C4

Amax

m2

1423

2863

998

868

1000

150

P1: A/Amax Q1

% kW

93.7 10,584

93.1 22206.6

97.7 7966

95.0 12,209

100.0 14,574

8.7 249

P2: A/Amax Q2

% kW

96.4 10,926

100.0 23,855

100.0 8257

100.0 12,020

95.7 13,758

100.0 2831

P3: A/Amax Q3

% kW

100.0 11,431

88.1 21,020

100.0 8319

100.0 11,137

96.2 13,888

0.0 0

Table 15 Utility consumption of optimal heat exchanger network Period

Hot utility consumption (kW)

Cold utility consumption (kW)

P1 (SOR) P2 (MOR) P3 (EOR)

27,033 28,609 25,025

40,756 43,038 40,770

are used instead as heat exchanger hot and cold stream inlet and outlet temperatures. Tables 13 and 14 show the new values for areas, heat loads, and temperatures in the network structure. The updated HEN structure showing heat transfer areas for all matches is shown in Fig. 16. Fig. 16 shows the final HEN, containing 12 heat exchangers and featuring a total annual cost of 6.32 million ¥. This network features two binary stream splits of stream H1. After careful inspection, it might be worthwhile investigating the effect of allowing the split to cross both stages without mixing between the stages. In this way, there will be only one stream split and it might be possible to combine the two heat exchangers which match streams H1 and C1 in both stages. The total utility consumption per period is shown in Table 15. 4.4. Model performance

TAC. Table 13 contains the maximum areas, the area fractions, the heat duties, the flows and temperatures for each process-to-process heat exchanger in the superstructure. Next, Table 14 shows the areas, the area fractions and the heat duties of the utility heat exchangers. Again, instead of using th (k) and th (k + 1), which were assigned to stage temperatures and were equal to heat exchanger inlet and outlet temperatures by the assumption of iso-thermal mixing, now th (in) and th (out)

The solution to the improved MINLP–NLP program was obtained in 45 s using an AMD Athlon姠 1.66 GHz processor. The solution to Aaltola’s model was obtained by the same processor in also 45 s. The case study used in this work is of relatively small scale, and therefore it is necessary to measure the performance of both models versus problems of larger scale. Table 16 shows the performance of the new model when the number of periods is increased or decreased. The standard case is the case investigated in this work, having a total of three hot streams, four cold streams, three periods of operation and two stages. Increasing the number of periods increases the solving time, the number of variables and the number of equations. However, the number of binary variables is not a function of the number of periods and therefore remains constant.

W. Verheyen, N. Zhang / Chemical Engineering Science 61 (2006) 7730 – 7753 Table 16 Sensitivity of model performance versus number of periods Number of periods

Solve time (s)

Single equations

Single variables

Discrete variables

Single 2 Standard 4 5 6

27 30 42 50 104 157

572 1143 1714 2285 2856 3427

280 485 690 895 1100 1305

43 43 43 43 43 43

Table 17 Sensitivity of model performance versus number of stages Number of stages

Solve time (s)

Single equations

Single variables

Discrete variables

Single stage Standard 3 stages 4 stages

34 42 46 65

1168 1714 2260 2806

465 690 915 1140

31 43 55 67

Table 18 Sensitivity of model performance versus number of streams Number of streams

Solve time (s)

Single equations

Single variables

Discrete variables

Standard Add H4 and C5 Add H5 and C6 Add H6 and C7

42 94 170 335

1714 2710 3928 5368

690 1088 1574 2148

43 69 101 139

Table 17 shows the sensitivity of the model performance for varying number of stages, with two stages as standard. The table shows that increasing the number of stages slightly increases the solving time, while also increasing the number of equations, number of variables, and number of binary variables. A single stage would not be a very efficient way to design a proper HEN, but increasing the number of stages to four or more would increase the complexity of the model without providing any improvement in the model results. Finally, Table 18 shows the sensitivity of the model versus the number of hot and cold streams. After every run, the number of hot streams and cold streams has been incremented by one and the performance of the model has been assessed. As can be seen, the solving time increases substantially after increasing the number of hot and cold streams. Also there is a large increase in number of single equations and variables, with binary variables reaching 140 for a set of 13 streams in three periods and two stages. 4.5. Comments This case study showed the application of the newly developed model to a typical case study, the design of a multi-period HEN for the VGO hydrotreater unit of an oil refinery. First, a

7751

comparison was made between a single period model which allowed non-linear constraints and featured a linear objective function, and the existing single period model which features only linear constraints and a non-linear and non-convex objective function. From this comparison it was shown for both formulations to predict similar results. Then the proposed multi-period model was tested rigorously against Aaltola’s multi-period simultaneous. It is shown that, by accounting for the cost maximum areas in the objective function of the MINLP model, one can reach solutions which are more economically attractive than the results obtained by Aaltola’s multi-period model. It is also shown that, for the scope of this case study, the average TAC is approximately 5% lower than the average TAC obtained from Aaltola’s model. 5. Conclusions In this work, a simultaneous MINLP model has been developed for the design of HENs in multi-period operation. Simultaneous heat exchanger design methods have been chosen to be extended to multiple periods because they can provide a rigorous framework that can obtain good solutions with a moderate degree of complexity. The NLP accounts for the maximum area per period at the expense of allowing non-linear functions in the set of constraints. The assumption of isothermal mixing has been overcome by an NLP improvement model that allows mixing at different temperatures. It has been shown in this work that the superstructure of Yee and Grossmann (1990) can be applied for designing multiperiod HENs. Furthermore, the assumption of average area in the objective function of the MINLP model developed by Aaltola (2002) has been removed, and it has been shown by a case study that better results can be obtained by allowing the maximum area into the objective function of the MINLP model. As the problem size increases, the simultaneous model becomes more complex and more difficult to solve, especially when non-linear terms are present in the set of constraints. Therefore, a sensitivity analysis showed the satisfactory performance of the improved model with respect to an increase in number of periods, number of stages, or number of streams. Notation Abbreviations CC GA GCC HEN LMTD LP MILP MINLP NLP SA TAC

composite curves genetic algorithm grand composite curve heat exchanger network logarithmic mean temperature difference linear programming mixed-integer linear programming mixed-integer non-linear programming non-linear programming simulated annealing total annual costs

7752

TS VGO

W. Verheyen, N. Zhang / Chemical Engineering Science 61 (2006) 7730 – 7753

Tabu search vacuum gas oil

zcu1 (j )

Indices i j k p

hot process stream or hot utility cold process stream or cold utility stage number or temperature interval period of operation

Sets CP HP PR ST

set of a cold process stream j set of a hot process stream i set of a operation period, p=1, . . . , NOP set of a stage in the superstructure, k = 1, . . . , NOK

Parameters AF B C CP Ccu Cf Chu Cw DOP(p) DTup epsi Fi (i, p) Fj (j, p) HUup MAXarea (i, j, k)

NOK NOP Qup SplitC Split H Tiin Tiout Tjin Tjout U

z1 (i, j, k)

zhu1 (i) Variables A(i, j, k) Acu (i) Ahu (j ) dtcu (j, k, p) dthu (i, k, p)

annualisation factor, dimensionless exponent for area cost, dimensionless area cost coefficient for heat exchanger, ¥/unit specific heat capacity, kJ/kg K per unit cost for cold utility, ¥/unit fixed charge for heat exchanger unit, ¥/unit per unit cost for hot utility, ¥/unit weight parameter for decrease of areas, dimensionless duration of period, dimensionless upper bound on temperature difference, ◦C exchanger minimum approach temperature, ◦ C heat capacity flow rate of hot stream i, kW/K heat capacity flow rate of cold stream j , kW/K upper bound on total hot utility available, kW upper bound on heat transfer area for exchanger connecting streams i and j in stage k, m2 number of stages, dimensionless number of periods, dimensionless upper bound on heat exchange capacity, kW existence of split on cold stream j at stage k, dimensionless existence of split on hot stream i at stage k, dimensionless inlet temperature of hot stream, ◦ C outlet temperature of hot stream, ◦ C inlet temperature of cold stream,◦ C outlet temperature of cold stream,◦ C overall heat transfer coefficient,kW/m2 K

existence of match (i, j ) in stage k,dimensionless existence of cold utility match for hot stream i,dimensionless existence of hot utility match for cold stream j ,dimensionless

dt(i, j, k, p) fc (i, j, k, p) fh (i, j, k, p) obj q(i, j, k, p) qcu (i, p) qhu (j, p) s(i, j, k, p) smin (i, j, k) tcs (i, j, k, p) ti(i, k, p) ti(i, k + 1, p) ths (i, j, k, p) tj(j, k, p) tj(j, k + 1, p)

maximum area for match of hot stream i and cold stream j in stage k,m2 maximum area for match of hot stream i and cold utility,m2 maximum area for match of cold stream j and hot utility, m2 temperature difference for match of hot stream i and cold utility in period p, ◦ C temperature difference for match of cold stream j and hot utility in period p, ◦ C temperature difference for match (i, j ) at temperature location k in period p, ◦ C heat capacity flow rate of cold stream fraction related to exchanger i, j, k, kW/K heat capacity flow rate of hot stream fraction related to exchanger i, j, k, kW/K objective function, total annualised cost, dimensionless heat exchanged between hot streami and cold stream j in period p, kW heat exchanged between hot stream i and cold utility in period p, kW heat flow exchanged between cold stream j and hot utility in period p, kW slack variable, dimensionless minimum of s(i, j, k, p) for match i, j, k, dimensionless temperature of cold stream fraction after exchanger i, j, k in period p, ◦ C temperature of hot stream i at hot end of stage k in period p, exchanger inlet, ◦ C temperature of hot stream i at hot end of stage k in period p, exchanger outlet, ◦ C temperature of hot stream fraction after exchanger i, j, k in period p, ◦ C temperature of cold stream j at hot end of stage k in period p, exchanger outlet, ◦C temperature of cold stream j at hot end of stage k in period p, exchanger inlet, ◦C

Binary variables zcu (i) zhu (j )

existence of cold utility match for stream j , dimensionless existence of hot utility match for stream i, dimensionless

W. Verheyen, N. Zhang / Chemical Engineering Science 61 (2006) 7730 – 7753

z(i, j, k)

existence of match between i and j in stage k, dimensionless

References Aaltola, J., 2002. Simultaneous synthesis of flexible heat exchanger network. Applied Thermal Engineering 22, 907–918. Biegler, L.T., Grossmann, I.E., Westenberg, A.W., 1997. Systematic Methods of Chemical Process Design. Prentice-Hall Inc., NJ. Bjork, K.-M., Westerlund, T., 2002. Global optimization of heat exchanger network synthesis problems with and without the isothermal mixing assumption. Computers & Chemical Engineering 26 (11), 1581–1593. Briones, V., Kokossis, A.C., 1999. Hypertargets: a conceptual programming approach for the optimization of industrial heat exchanger networks—I, Grassroots design and network complexity. Chemical Engineering Science 54 (4), 519–539. Cerda, J., Westenberg, A.W., Mason, D., Linhoff, B., 1983. Minimum utility usage in heat exchanger network synthesis. A transportation problem. Chemical Engineering Science 38 (3), 373–387. Chen, J.J.J., 1987. Comments on improvements on a replacement for the logarithmic mean. Letters to the editor. Chemical Engineering Science 42 (10), 2488–2489. Chen, C.L., Hung, P.S., 2004. Simultaneous synthesis of flexible heatexchange networks with uncertain source-stream temperatures and flow rates. Industrial & Engineering Chemistry Research 43, 5916–5928. Ciric, A.R., Floudas, C.A., 1991. Heat exchanger network synthesis without decomposition. Computers & Chemical Engineering 15 (6), 385–396. Floudas, C.A., 1995. Nonlinear and Mixed-Integer Optimization. Oxford University Press, New York, pp. 359–372. Floudas, C.A., Grossmann, I.E., 1986. Synthesis of flexible heat exchanger networks for multiperiod operation. Computers & Chemical Engineering 10 (2), 153–168. Floudas, D.A., Grossmann, I.E., 1987. Automatic generation of multiperiod heat exchanger network configurations. Computers & Chemical Engineering 11 (2), 123–142. Floudas, C.A., Ciric, A.R., Grossmann, I.E., 1986. Automatic synthesis of optimum heat exchanger network configurations. A.I.Ch.E. Journal 32 (2), 276–290. Hohmann, E.C., 1971. Optimum networks for heat exchange. Ph.D. Thesis, University of Southern California. Linnhoff, B., 1979. Thermodynamic analysis in the design of process networks. Ph.D. Thesis, The University of Leeds, England. Linnhoff, B., Flower, J.R., 1978a. Synthesis of heat exchanger networks I: systematic generation of energy optimal networks. A.I.Ch.E. Journal 24 (4), 633–642. Linnhoff, B., Flower, J.R., 1978b. Synthesis of heat exchanger networks II: evolutionary generation of networks with various criteria of optimality. A.I.Ch.E. Journal 24 (4), 633–642. Marselle, D.F., Morari, M., Rudd, D.F., 1982. Design of resilient processing plants—II, design and control of energy management systems. Chemical Engineering Science 37 (2), 259–270. Masso, A.H., Rudd, D.F., 1969. Synthesis of system design II—heuristic structuring. A.I.Ch.E. Journal 15 (1), 10–17.

7753

Papalexandri, K.P., Pistikopoulos, E.N., 1994a. Synthesis and retrofit design of operable heat exchanger networks—I, flexibility and structural controllability aspects. Industrial & Engineering Chemistry Research 33 (7), 1718–1737. Papalexandri, K.P., Pistikopoulos, E.N., 1994b. Synthesis and retrofit design of operable heat exchanger networks—II, dynamics and control structure considerations. Industrial & Engineering Chemistry Research 33 (7), 1738–1755. Papoulias, S.A., Grossmann, I.E., 1983. A structural optimisation approach in process synthesis—II, heat recovery networks. Computers & Chemical Engineering 7 (6), 707–721. Paterson, W.R., 1984. A replacement for the logarithmic mean, shorter communication. Chemical Engineering Science 39 (11), 1635–1636. Saboo, A.K., Morari, M., Woodcock, D.C., 1985. Design of resilient processing plants-VIII, a resilience index for heat exchanger networks. Computers & Chemical Engineering 40 (8), 1553–1565. Smith, R., 2005. Chemical Process Design and Integration, second ed. Wiley, New York, pp. 5–8, 357–438. Swaney, R.E., Grossmann, I.E., 1985a. An index for operational flexibility in chemical process design, part I: formulation and theory. A.I.Ch.E. Journal 31 (4), 621–630. Swaney, R.E., Grossmann, I.E., 1985b. An index for operational flexibility in chemical process design, part II: computational algorithms. A.I.Ch.E. Journal 31 (4), 631–641. Tantimuratha, L., Asteris, G., Antonopoulos, D.K., Kokossis, A.C., 2001. A conceptual programming approach for the design of flexible heat exchanger networks. Computers & Chemical Engineering 25 (4–6), 887–892. Viswanathan, J., Grossmann, I.E., 1990. A combined penalty function and outer-approximation algorithm for MINLP optimization. Computers & Chemical Engineering 14 (7), 769–782. Yee, T.F., Grossmann, I.E., 1990. Simultaneous optimisation models for heat integration—II, heat exchanger network synthesis. Computer & Chemical Engineering 14 (10), 1165–1184. Yee, T.F., Grossmann, I.E., Kravanja, Z., 1990. Simultaneous optimisation models for heat integration—I, area and energy targeting and modelling of multi-stream exchangers. Computers & Chemical Engineering 14 (10), 1151–1164. Zhu, X.X., 1995. Automated synthesis of HENs using block decomposition and heuristic rules. Computers & Chemical Engineering 19 (Suppl.), S155–S160. Zhu, X.X., O’Neil, B.K., Roach, J.R., Wood, R.M., 1995a. A new method for heat exchanger synthesis using area targeting procedures. Computers & Chemical Engineering 19 (2), 197–222. Zhu, X.X., O’Neil, B.K., Roach, J.R., Wood, R.M., 1995b. A method for automated heat exchanger synthesis using block decomposition and nonlinear optimisation. Transactions of the Institution of Chemical Engineers 73 (A), 919–930. Zhu, X.X., O’Neil, B.K., Roach, J.R., Wood, R.M., 1995c. Area targeting methods for the direct synthesis of heat exchanger networks with unequal film coefficients. Computers & Chemical Engineering 19 (2), 223–239. Zhu, X.X., 1997. Automated design method for heat exchanger network using block decomposition and heuristic rules. Computers & Chemical Engineering 21 (10), 1095–1104.