Structure and Area Optimization of Flexible Heat Exchanger Networks

Jun 16, 2014 - synthesized at the nominal operating point and renewed by the topological union with the structure of the critical point and the improv...
1 downloads 0 Views 5MB Size
Download PDF
Article pubs.acs.org/IECR

Structure and Area Optimization of Flexible Heat Exchanger Networks Jilong Li, Jian Du,* Zongchang Zhao, and Pingjing Yao State Key Laboratory of Fine Chemicals, Chemical Engineering Department, Dalian University of Technology, Dalian 116024, Liaoning, China S Supporting Information *

ABSTRACT: A novel method is outlined for flexible heat exchanger network synthesis including nonconvex problems. The presented method is sequentially implemented by two main steps: structure synthesis and area optimization; nevertheless, the optimization of heat exchanger areas can still react on the structure to gain the global optimal solution. The structure is initially synthesized at the nominal operating point and renewed by the topological union with the structure of the critical point and the improved heat transfer loops disconnection strategy. For area optimization, an iterative approach with strong robustness is proposed based on the influences of heat exchanger areas on flexibility index and total annual cost, respectively. The direction matrix method is employed to provide the operational flexibility of the network and the critical operating point. Two examples with nonconvex feasible regions have been studied, and the results well demonstrate the effectiveness of the proposed approach. resulting network is not strictly guaranteed since only a finite number of operating conditions are considered during the design especially for nonconvex problem. Verheyen and Zhang16 retrofitted the Aaltola’s model with a combined MINLP-NLP model which uses maximum area per period in the area cost calculation of the MINLP objective function and removes slack variables and weighed parameters from the existing NLP model. Konukman et al.17 also presented a noniterative, superstructure-based MILP formulation for simultaneous flexibility targeting and minimum-utility. Nevertheless this formulation still presumes that the feasible region in the space of uncertain input parameters is convex. Taking advantage of the previous methods, Chen and Hung18 proposed a three-step iterative approach to design flexible multiperiod HENs based on the mathematical formulation of Aaltola15 and the active set strategy13 for testing the flexibility of HEN. In their work, the maximum area was obtained with a discontinuous function according to some specific operating conditions, which may probably miss the optimal solution. The iterations with superstructure also limit the problem scale. On the other hand, Ma et al.19 generated the flexibility-qualified structure by combining the HEN topologies of all the selected operating periods and then improved the network to overcome the problem of oversynthesis. Nevertheless, all the multiperiod-based methods have a common defect that the formulation, especially the objective function, is limited by the selected operating points. This causes most multiperiod based methods to solve only the problems with a convex feasible region. Although some methods include the flexibility tests to guarantee the feasibility for nonconvex problems, such as that of Floudas and Grossmann11 or Chen

1. INTRODUCTION As a heat exchanger network (HEN) consumes a great fraction of energy in a chemical process system, the network’s synthesis and optimization are vitally important. The conventional HEN synthesis is performed under nominal conditions with fixed operating parameters. However, the network configuration with high integration always makes the process more difficult to control and operate1 with variable conditions (for example, feed temperatures, flow rates, cost data, heat transfer coefficients, etc.) which is inevitable in practical operating environments. Therefore, maintaining operable and optimum synthesis methods of the flexible HENs under uncertain parameters is very important. Since Marselle et al.2 first developed the concept of resilient HENs that can tolerate uncertainties in temperatures and flow rates, a significant amount of synthesis methods of flexible HENs3−9 has been developed. Reviews of the research of flexible HENs can be found in the literature.1,10 Here, a timeline of innovation and major discoveries of flexible HEN synthesis are also sketched. Floudas and Grossmann11 proposed a sequential synthesis method for a flexible HEN, which predicts stream matches with the multiperiod mixed-integer linear programming (MILP) trans-shipment12 and deduces the network configuration with the active set strategy13 to guarantee the desired flexibility. The decomposition of the problem into different stages significantly reduces the size of the problem, but does not take rigorously into account the trade-off between area, number of units, and energy. Papalexandri and Pistikopoulos14 introduced a twostage stochastic mathematical programming for the synthesis and retrofit of flexible and structurally controllable HENs while simultaneously measuring design feasibility. Aaltola15 developed a multiperiod mixed-integer nonlinear programming (MINLP) model based on the stage-wise superstructure representation. This method simultaneously handles the minimum total annual cost (TAC) and the flexibility without relying on sequential decomposition. However, the operational feasibility of the © 2014 American Chemical Society

Received: Revised: Accepted: Published: 11779

March 26, 2014 May 29, 2014 June 16, 2014 June 16, 2014 dx.doi.org/10.1021/ie501278c | Ind. Eng. Chem. Res. 2014, 53, 11779−11793

Industrial & Engineering Chemistry Research

Article

and Hung,18 the feasible synthesis cannot be accomplished accurately and globally in the whole feasible region. Pistikopoulos and Ierapetritou20 approximated the objective with flexible constraints through the Gaussian quadrature formula. However, this formulation yields a very complex model that is limited to small problems. Pintarič and Kravanja21 discussed the objective functions of process synthesis under uncertainty and developed the reduced dimensional stochastic strategy for the approximation of the expected objective function. Though it is still calculated over a reduced set of the extreme points, the approach theoretically treats the flexible HEN synthesis problems via dispensing the multiperiod presumption. Pintarič and Kravanja22 further expounded the mathematical model for flexible design problems and provided a synthesis method of flexible processes based on the identification of critical points. Their approach identifies critical values of uncertain parameters a priori by the separate maximization of each design variable, together with simultaneous optimization of the economic objective function. Three alternative methods are proposed respectively: the Karush− Kuhn−Tucher (KKT) formulation, the iterative two-level method, and the approximate one-level method. However, the former two are both impractical for large and complex problems, while the last one somewhat oversimplify the optimization of the design variables. Moreover due to the matrix-including partial derivative, the method for determination of critical points is hard to be solved automatically with a program, which further limits its application range. In this study, a novel stepwise method for flexible HEN synthesis is proposed. The mathematical formulation of the flexible HEN is retrofitted on the basis of the general model of Pintarič and Kravanja22 without multiperiod limitation. The maximum heat exchanger areas are defined as independent variables for global optimization instead of the conventional maximized operation. Considering the calculation complexity, the proposed method decomposed the flexible HEN synthesis problem into two main steps: structure synthesis and area optimization. The structure synthesis is accomplished through an iterative method among critical points with a flexibility test similar to the approach of Chen and Hung.18 However, the new structure is generated with the topological union19 and the improved disconnection strategy of heat transfer loops for flexible problems. The critical points are provided by the direction matrix strategy in our previous work23 instead of the empirical strategy. The area optimization is implemented after the flexible-qualified structure is obtained. A quantified iterative formulation with strong robustness is developed on the basis of the influences of the heat exchanger areas on the flexibility index and the TAC. As the presented method is established directly on the general formulation of flexible HEN synthesis without vertices assumption, constraints requirement or linearization, it can effectively synthesize flexible HENs with nonconvex feasible regions and all kinds of uncertain parameters, even the one defined with a formulation. Examples from literature are studied with this new method, and the results are quite satisfactory.

min E TAC d,z,x θ

s.t. h(d, z , x , θ) = 0 g (d , z , x , θ ) ≤ 0 d, z, x ∈ R, θ L ≤ θ ≤ θU

(1)

where E denotes the expectation of the TAC in the domain of operating points; h and g are the vectors of equations (such as energy and mass balance) and inequalities (typically physical operating limits or product specifications). All the variables in the formulation are expressed in four categories: d is the vector of design variables that define the network configuration and the sizes of equipment. It contains both continuous and binary variables. θ is the vector of uncertain parameters such as temperatures and flow rates; z, the vector of control variables, represents the degrees of freedom during the operation, which can be adjusted for different realization of θ; x is the vector of state variables, a subset of the remaining variables with the same dimension as h. For a process, it can be expressed as functions of other variables (d, z, θ) using the equality constraints. Equation 1 is an infinite problem that cannot be solved directly unless the expression of the expectation is definite. Approaches for the expressions have been discussed in the literature and they can generally be classified into the following three categories. (a) The simplest one is taking the nominal point as the expectation, ETAC = TAC(θ N)

(2)

θ

where N means the nominal point. (b) Another widely used method approximates the expectation with a number of uncertain parameter realizations (eq 3). It is also regarded as the multiperiod problem. E TAC = θ

∑

p(n) TAC(θ(n)) (3)

∀ n ∈ VT

where n denotes the operating point n. VT = {0, 1, ..., nv} is the set of the operating conditions considered, where 0 means the nominal condition. p(n) represents the probabilities corresponding to the operating points. (c) To describe the flexible problem exactly, the probabilistic distribution functions of the uncertain parameters are employed for the expected value. ND

E TAC = θ

∫θ∈Ω p(θ) TAC(θ) dθ ≈ ∑ viTAC(θi) i=1

(4)

where Ω is value space of the uncertain parameter vector θ. However, in practice the integral still needs a discrete approximation through some numerical methods such as the Gaussian integration formula,20 for the rigorous analytic solution is hardly to be obtained. vi is the corresponding coefficients determined by the numerical method, and ND is the number of the points needed. On the other hand, as the TAC of the HEN consists of the utility cost (operating cost) and the equipment cost (capital cost), the objective function in eq 1 can be described as

2. PROBLEM FORMULATION Taking the total annual cost as the objective, a mathematical model for synthesis problems under uncertainty can be generally expressed in the following formulation: 11780

dx.doi.org/10.1021/ie501278c | Ind. Eng. Chem. Res. 2014, 53, 11779−11793

Industrial & Engineering Chemistry Research

Article

3. SYNTHESIS OF FLEXIBLE HEN Unlike the flexibility analysis problems (Grossmann and Floudas, 13 Swaney and Grossmann, 24 Konukman and Akman, 25 Petracci et al.,26 Li and Chang,27 Adi and Chang28), the goal of flexible synthesis is to provide the optimal design over the uncertain space rather than identify the maximum violation of uncertain constraints for fixed designs. Though the mathematical model eq 9 is established for the optimization target, the equations h and the inequalities g are necessary for the solution, and they are mainly decided by the structure. Two categories of approaches can solve this problem. The first one is building a superstructure that directly includes the flexible target. It contains not only all the possible structures but also all the detailed flexible constraints, and deals with them simultaneously. So it is also called the simultaneously optimization method. This approach is rigorous and theoretically can find the optimal solution; however, it is too complex to realize. Even for the HENs with fixed conditions, the superstructure method is often so complex that some simplifications have to be made. The other approach treats the structure and the area design stepwise. It starts with a candidate structure derived by finite number of operating conditions, and then the flexible constraints are established and further optimization of structures and areas is made sequentially until the target is accomplished. Since the candidate structure simplifies the flexible constraints, the search space can be drastically reduced. In this study, the second approach is employed. The problem is solved by first finding a feasible structure for the uncertain constraints and then making the area optimization based on the obtained structure. 3.1. Structure Synthesis. The target of structure synthesis is to find a relatively optimal structure with enough flexibility. As discussed above, it begins with finite operating points and keeps iterating until the flexibility requirement is satisfied. The iterative steps are listed as follows and also illustrated in Figure 1.

min E TAC = E(UC(d(n), z(n), x(n), θ(n))) d,z,x θ

θ

+ E EC(d(n), z(n), x(n), θ(n)) θ

n ∈ VT

(5)

where UC indicates the utility cost and EC means the equipment cost, with the following discretized expressions, UC(d(n), z(n), x(n), θ(n)) 1 ( ∑ (CCU = NV + 1 ∀ n ∈ VT + CHU

∑

∑

Q CU , i(d(n), z(n), x(n), θ(n))

i ∈ HP

(n)

Q HU , j(d , z(n), x(n), θ(n)))) (6)

j ∈ CP

EC(d(n), z(n), x(n), θ(n))

∑

=

C EX( max Ak (d(n), z(n), x(n), θ(n))) ∀ n ∈ VT

k ∈ EP

(7)

where HP, CP, and EP are the set of hot process streams, cold process streams, and heat exchangers (including utility heaters and coolers), respectively. CEX denotes the heat exchanger cost function. As a min−max problem with uncertain constraints, it is very difficult to be solved directly especially when the number of the operating points is large. In this study, the group of maximum heat exchanger areas are selected from the design variables d, and renamed as the vector Am. Accordingly, the vector A represent the actual areas calculated by the data of the particular operating points. Thus, the maximum constraints eq 7 can be transformed into the following formulation, E EC(d(n), z(n), x(n), θ(n)) = θ

Amk ≥

∑

C EX(Amk )

k ∈ EP

k ∈ EP

Ak (d(n), z(n), x(n), θ(n)) ∀ n ∈ VT

(8)

By extending the operating point range in the inequality constraints of eq 8 to the whole feasible space and considering the expression of the expectation, eq 1 is changed into the following form, which is the model used in this study. min E TAC = E UC(d(n), z(n), x(n), θ(n)) + d,z,x θ

θ

∑

CEX(Amk )

k ∈ EP

n ∈ VT s.t. h(d, z , x , θ) = 0 g (d , z , x , θ) ≤ 0 Am ≥ A (d, z , x , θ) ∀ θ∈Ω

d, z, x ∈ R, θ L ≤ θ ≤ θU (9)

Figure 1. Proposed approach for the structure synthesis of a flexible HEN.

Compared with the traditional multiperiod formula as eq 7, the maximum heat exchanger areas needed (Am) are independent of the operating points in this model. They are used as variables to be optimized. Equation 9 can accurately guarantee the flexibility of the HEN since the feasible region is no longer affected by the selection of operating points. It indicates this new model can treat both convex and nonconvex HEN problems effectively.

(1) Synthesize the HEN at nominal operation point. (2) Evaluate the flexibility index of the HEN. If it is less than 1, then go to step 3, or else the iteration ends. (3) Point out the critical point (if there are more than one, then just select any one) and synthesize the HEN at this operating point. 11781

dx.doi.org/10.1021/ie501278c | Ind. Eng. Chem. Res. 2014, 53, 11779−11793

Industrial & Engineering Chemistry Research

Article

Figure 2. Generation of flexible HEN structure with topological union: (a) topology of operating point 1; (b) topology of operating point 2; (c) topology of the combined HEN.

variable to be optimized. Since the match of hot and cold streams is completed vertically in each enthalpy interval separately, the pattern of partitioning enthalpy intervals affects the network structure a lot. Thus, the merge array is proposed to further optimize the structure of the HEN. The detailed presentation and demonstration of the improved pseudo-T−H diagram method for synthesizing HENs can be found in the literature.29 3.1.2. Topological Union for the New Structure of the HEN. Since the uncertain parameters vary at different operating points, the HEN structure of each operating point and the heat exchanger area for the same match are probably different. To get a feasible HEN for the two selected operating points at each iteration, the topological union approach presented by Ma et al.19 is employed. A brief review of this approach is shown below through a simple example, while the detailed description can be found in the literature.19 As illustrated in Figure 2, structures a and b are two HEN structures with different matches for two operating points. EXi,j,l,n represents the match between hot stream i and cold stream j at interval l in the topology of operating point n. The topological union is the set of all distinct elements in the collection of topological space. In the flexible HEN problems, it means the new structure is constructed by uniting all the matches existing in any of the original ones. EX1,2,1, for instance, only existing in structure a, is involved in the combined structure c. Moreover, the maximum heat exchanger area for the same match is selected to be the value for the combined structure to ensure the new network is feasible for all the original operating points. This group of heat exchanger areas are used as reference values for the latter area optimization as discussed above, if the structure obtained with the topological union is qualified. 3.1.3. Heat Transfer Loops Disconnection. The topological union of the HEN structures would probably bring the structure redundancy reflected in heat transfer loops, so that it is necessary to disconnect them. However, the disconnection strategy is different from that with certain conditions. The conventional disconnection strategy is to remove the heat exchanger with minimum heat load in the loop and distribute the heat load to other heat exchangers sequentially. First there is a simple example composed by 2 hot streams and 2 cold streams with certain conditions. The initial network is shown in Figure 3a. It can be found that there is a heat exchanger loop, denoted as (EX1, EX2, EX4, EX3). The reason it begins with EX1 is that EX1 is the one to be removed for its minimum heat load and it is convenient to distribute the heat load in this sequence. When EX1 is deleted, the heat loads of

(4) Generate the new structure with the network of critical point and that of last iteration by topological union and disconnect the heat transfer loops. (5) Evaluate the flexibility index of the new result. If it is less than 1, then go to step 3, or else the iteration is terminated. The flexibility index F is a quantitative performance indicator for a flexible process. It reveals the maximum scaled deviation of the uncertain parameters from their nominal values with a guarantee of feasible operation. If F ≥ 1, it means the network can satisfy the design requirement, or else it means not. Details about the flexibility index can be found in the literature.24 Note, in the structure synthesis, the flexibility specifically indicates the structure flexibility excluding the heat exchanger area constraints, because the main target of this step is to find a feasible structure for the whole region. The heat exchanger areas obtained at each synthesis only represent the optimal solution for the single operating point. They just provide initial values for the further optimization. Thus, it is assumed that all the heat exchangers can have large enough areas in the flexibility analysis of this step. The detailed description of the structure flexibility and the overall flexibility is shown in section 4. 3.1.1. Synthesis of the HEN with Certain Conditions. A pseudotemperature enthalpy diagram (pseudo-T−H diagram)based method29 is applied to construct the structure of the HEN, since it can solve the complex constraints in heat transfer processes efficiently. As it is named, the pseudo-T−H diagram30,31 is attained with the pseudotemperature of streams consisting of the real temperature and heat transfer temperature difference contribution value. For hot stream i TiP = Ti − ΔTci

(10)

For cold stream j T jP = Tj + ΔTcj

(11)

where superscript P means the pseudotemperature, and ΔTci, ΔTcj are the heat transfer temperature difference contribution values of the hot and cold streams, respectively. In the pseudoT−H diagram, the influences of the temperature levels, the heat transfer coefficients, and the energy losses are all expressed with the heat transfer temperature difference contribution values which are used as decision variables. Besides, the merge array, describing the mergence situation of the enthalpy intervals in the pseudo-T−H diagram, is also employed as the decision 11782

dx.doi.org/10.1021/ie501278c | Ind. Eng. Chem. Res. 2014, 53, 11779−11793

Industrial & Engineering Chemistry Research

Article

Figure 3. Conventional disconnection strategy for HEN: (a) original HEN with the heat transfer loop; (b) disconnection by removing EX1.

EX2 and EX3 need to be increased and that of EX4 needs to be decreased as shown in Figure 3b. Nevertheless, the treatment does not always work in a flexible problem. This is illustrated with a similar example. Figure 4a gives the initial network. In this example, the heat load of H1 (QH1) is the uncertain parameter and a utility heater is added on C1 since its heat load is increased. The numbers under the heat exchangers represent the nominal heat loads of them. With an adjustment of the heat load of the heater (QHU), the process can maintain operation even if QH1 changes in its range (±30). If the network is treated exactly as the former one, the structure in Figure 4b is obtained. In this network, if QH1 >195, the process is infeasible no matter what the value of QHU is. That indicates the flexibility of the network is declined by the disconnection. To keep the original flexibility, a utility cooler should be added as shown in Figure 4c. Therefore, the disconnection strategy in flexible problems should be retrofitted by one more rule: If the stream where the heat exchanger is removed in the disconnection does not contain any utility heat exchanger, then a utility heat exchanger is added to it. Although this new strategy will sometimes generate an overneeded structure, it is easier to handle the redundant utility heat exchangers, for each stream has only one utility heat exchanger but maybe several stream-to-stream heat exchangers. Anyway, this problem will be solved in the next section. If a heat exchanger is actually not needed, its area will be very close to 0 after the area optimization, so that it can be regarded as a nonexistent one. The significance of the new strategy is to make sure the flexibility of the network will absolutely not decline. Without this strategy, it may not find a flexible enough structure in finite iterations. 3.2. Area Optimization. Unlike the multiperiod based methods of finite operating points, the area optimization in the

Figure 4. Proposed disconnection strategy for flexible HEN. (a) Original HEN with the heat transfer loop; (b) disconnection with the conventional strategy; (c) disconnection with the proposed strategy.

formulation eq 9 will require too much computation cost with conventional algorithms due to the global flexibility constraints. Thus, an iterative approach is proposed here to accomplish the further optimization of heat exchanger areas based on the qualified structure obtained. In flexible HEN problems, despite the uncertainty of some parameters, the expected ranges of the uncertain parameters are already known. Considering the definition of the flexibility index, the following property can be obtained (after the structure synthesis, the same hereinafter): If F > 1, then some of the heat exchanger areas must be overneeded; If F < 1, then some of the heat exchanger areas must be not large enough. On the other side, in the HEN a larger heat exchanger area always results in more cost with a decided configuration. Thus, the optimization target here is to find a group of the areas that can just satisfy the flexibility requirement (F = 1) with the least possible total annual cost. In the iterative way, the areas of heat exchangers are increased/decreased according to the flexibility 11783

dx.doi.org/10.1021/ie501278c | Ind. Eng. Chem. Res. 2014, 53, 11779−11793

Industrial & Engineering Chemistry Research

Article

index each time, while the amplitude of the increase/decrease is decided by the influence of the areas on the flexibility index and the total annual cost. The detail iteration rules are as follows: If F > 1, then all the heat exchanger areas (including heaters and coolers) are decreased. The larger F is, the more they should be decreased. If the area-decrease of one heat exchanger influences the flexibility index more than the TAC, then this heat exchanger should have a smaller area-decrease; or else the area of this heat exchanger should be decreased more. If F < 1, then all the heat exchanger areas (including heaters and coolers) are increased. The smaller F is, the more they should be increased. If the area-decrease of one heat exchanger influences the flexibility index more than the TAC, then this heat exchanger should have a larger area-increase; or else the area of this heat exchanger should be increased less. Despite the same variation trend (increase or decrease) of all the heat exchanger areas at each iteration, they can either increase or decrease compared to the initial values after several iterations due to the different variation amplitudes. So the iteration with these rules can play a role of optimization. It will be proved in the case study. According to the above rules, the qualitative formulation of the area iteration is proposed as follows, Amk ,new = Amk ·(1 + Δ)

k ∈ EP

⎧ IAC −(DF)· ·CA F ≥ 1 ⎪ ⎪ IAF Δ=⎨ ⎪ (DF)· IAF ·CA F < 1 ⎪ ⎩ IAC

(12)

where DF means the absolute value of the deviation of F from 1. IAC is the influence of areas on cost. IAF is the influence of areas on flexibility. CA represents an adjusting coefficient to make sure Δ in a reasonable range. For example, if Δ < −1, then a negative new area is obtained, which is incorrect. In addition, the proportional item (IAC/IAF or IAF/IAC) is just a tendency expression. To obtain the quantitative formulation, each item of Δ is analyzed sequentially in the following: 3.2.1. DF. The first item of Δ offers the basic fluctuation of the iteration. It can be easy to think of using (1 − F) as the first item, where F is the flexibility index in the current iteration. It can express not only the degree of the deviation of F from 1, but also its sign. But when F is close to 1, the DF is usually too small to generate enough fluctuation for the new generation (Amk,new). If we employ CA to amplify the fluctuation, the Δ may be out of control (such as Δ could be less than −1) when F is not so close to 1. Thus, this linear relationship is not suitable here. Considering the above discussion, the inverse tangent function is a reasonable choice. tan−1(η·(1 − F ))

Figure 5. Graph of the function eq 13: (a) η = 2.5; (b) η = 5.

satisfactory and adjustable basic fluctuation during the whole iterations. On the basis of the analysis of practical instances, η should be between 2 and 5. An oversmall η could lead to a premature convergence for the iteration. On the other hand, an overlarge η will cause divergence due to the too intense fluctuations. In practice, η usually begins with 3, and is adjusted according to the fluctuation degree of the flexibility index. 3.2.2. IAC/IAF or IAF/IAC. The second item of Δ is employed to provide different change-rates for each design variable (heat exchanger area). This is key reason that the iteration can realize the optimization. To describe the influence of area on cost and flexibility, the following two variables are defined,

(13)

where η is a coefficient to regulate the amplification degree in the region that F is close to 1. The greater η means the more amplified effect. The graph of this function is illustrated in Figure 5, where the straight line denotes the function f(F) = 1 − F. Figure 5a is the one of η = 2.5, while Figure 5b is the one of η = 5. As shown in Figure 5, an appropriate fluctuation is obtained around 1, and the curve is comparatively flat when F is farther from 1. These indicate that function eq 13 can provide a

ACk =

∂C ∂Amk

k ∈ EP

AFk =

∂F ∂Amk

k ∈ EP

(14)

(15)

In the practice, they are calculated with a numerical differentiation method so that it can be completed by computer programs automatically. Obviously, the flexibility index is 11784

dx.doi.org/10.1021/ie501278c | Ind. Eng. Chem. Res. 2014, 53, 11779−11793

Industrial & Engineering Chemistry Research

Article

usually between 0 and 2 and the cost is generally a very large number, so that the two corresponding partial derivatives also vary widely. If they are used directly, the result will probably be meaningless. Therefore, they are normalized separately before the further treatment. AC* and AF* express the variables after the normalization where the superscripted asterisk (∗) denotes the normalization. Moreover, for convenience, a new variable Rcf is introduced to represent the value of the second item, which has different expressions under the conditions F ≥ 1 and F < 1. First, we analyze the simple proportional formulation (only take IAC/IAF as example, IAF/IAC is similar): R cf, k =

AC*k AF*k

k ∈ EP (16)

Its graph is illustrated in Figure 6. Because of the normalization, ACk* and AFk* are all between 0 and 1. The graph indicates

Figure 7. Graph of the function eq 17 with α = 5.

Figure 6. Graph of the function eq 16.

Figure 8. Graph of the function eq 18 with α = 5 and ω = 10.

exchanger reduces the flexibility of the network a lot (such as AF* > 0.6), then the area of this heat exchanger should scarcely be decreased no matter how much the cost would be cut down. Only when the area affects the flexibility slightly (such as AF* < 0.6), is the change-rate of the area mainly decided by the influence of area on cost. This treatment can guarantee the convergence of the flexibility index in the iteration. On the basis of the calculation of cases, ω usually can be set as 10. The larger ω values enhance the influence of AF* and also result in a shorter iteration but probably a larger cost. On the other side, when F < 1, the formulation of Rcf can be expressed as eq 19 by a similar derivation as above.

that Rcf barely varies in most part of the domain but it varies sharply as AF* tends to 0. This formulation cannot provide a desired change-rate. Moreover when AFk* equals 0, error happens in eq 16. Therefore, some measures are needed to adjust Rcf, the change-rate. With use of the power function, eq 16 is changed into the following one, *

R cf, k =

α ACk

k ∈ EP * (17) α AFk α is a coefficient. The larger α leads to the steeper surface. The graph (α = 5) is shown in Figure 7. It can be seen that the new surface is much flatter than that obtained by the former formulation. Furthermore, with cases analyzed, it is found that the flexibility influence (IAF) usually impacts more than the cost influence (IAC) in the iteration. Thus, a weight factor ω is introduced and eq 17 is changed into eq 18.

*

R cf, k =

k ∈ EP * (19) α ACk The graph of function eq 19 is shown in Figure 9. The coefficients are also evaluated as α = 5, ω = 10. Conversely, in this situation the Rcf keeps a small value when AF* is close to 0. It means the heat exchangers with a tiny effect on the flexibility (such as AF* < 0.4) hardly have the area-increase, while the ones having signal influence on the flexibility (such as AF* > 0.4) appropriately have a higher area increase, which is also influenced by AC*. In addition, the α and ω in eq 18 and eq 19 can have different values, distinguished by the subscript 1 and 2 as in eq

*

R cf, k =

α ACk

*

(ωα)AFk

(ωα)AFk

k ∈ EP (18)

The graph (α = 5, ω = 10) is shown in Figure 8. Compared to Figure 7, the flexibility has a stronger effect: when AF* is close to 1, the Rcf is always small no matter what value AC* is. Its practical significance is that if the area-decrease of a heat 11785

dx.doi.org/10.1021/ie501278c | Ind. Eng. Chem. Res. 2014, 53, 11779−11793

Industrial & Engineering Chemistry Research

Article

Figure 9. Graph of the function eq 19 with α = 5 and ω = 10.

21. Besides, Rcf is also normalized as Rcf* for the further calculations. 3.2.3. AC. The third item of Δ is employed to restrict the range of Δ and it is defined as μ. To ensure the new area in the iteration is positive, Δ must be greater than −1. On the other hand, despite lack of clear upper limit of Δ, an overlarge Δ will cause the iteration to be instable and even divergent. On the basis of the above discussion, Δ can be expressed as * k·μ Δ = tan−1(η·(1 − F ))·R cf,

k ∈ EP

(20)

where F is generally between 0 and 2, and Rcfk* is between 0 and 1 due to the normalization. If the coefficients η and μ are confirmed, the graph of eq 20 can be gained. Since η has been valued as 2.5 and 5 in the former analysis (section 3.2.1), μ is set as 0.75 and 0.4, respectively, and the figures (Figure 10a and Figure 10b) are plotted. It can be seen that Δ has already been limited between −1 and 1. However, in practice μ can be appropriately amplified because Rcfk* is barely close to 1 after the normalization and F is hardly too far from 1. The value of μ can normally be set as 2, and it should also be adjusted larger or smaller according to the too intense or too slight fluctuations in the practical problems. Finally the formulation for iteration can be obtained, as eq 21.

Figure 10. Graph of the function eq 20. (a) η = 2.5 and μ = 0.75; (b) η = 5 and μ = 0.4.

the iteration terminates. However, in the case study, to illustrate the iteration trend a few more iterations are computed.

4. FLEXIBILITY ANALYSIS As the above discussion, the flexibilities in HEN are divided into two kinds, the structure flexibility and the overall flexibility. The former only shows the flexibility of the topological structure without considering the heat exchanger area, as expressed in the following formulation,

* k·μ) Amk ,new = Amk ·(1 + tan−1(η·(1 − F ))·R cf,

Rcf , k

⎧ AC*k ⎪ α1 F≥1 AF* ⎪ ⎪ (ω1α1) k =⎨ ⎪ (ω α )AF*k ⎪ 2 2* F