*Article* **Optimization of Heat Recovery Networks for Energy Savings in Industrial Processes**

**Jui-Yuan Lee \* and Po-Yu Chen**

Department of Chemical Engineering and Biotechnology, National Taipei University of Technology, 1, Sec. 3, Zhongxiao E. Rd., Taipei 10608, Taiwan

**\*** Correspondence: juiyuan@ntut.edu.tw; Tel.: +886-2-2771-2171 (ext. 2524)

**Abstract:** Among the pillars of decarbonization of the global energy system, energy efficiency plays a key role in reducing energy consumption across end-use (industry, transport and buildings) sectors. In industrial processes, energy efficiency can be improved by exploiting heat recovery via heat exchange between process streams. This paper develops a stage-wise superstructure-based mathematical programming model for the optimization of heat exchanger networks. The model incorporates rigorous formulation to handle process streams with phase change (condensation or evaporation), and is applied to a case study of an ethylene glycol production plant in Taiwan for minimizing utility consumption. The results show a compromise between steam savings and process feasibility, as well as how the model is modified to reflect practical considerations. In the preliminary analysis, with a substantial potential steam saving of 15,476 kW (28%), the solution involves forbidden matches that pose a hazard to the process and cannot be implemented. In the further analysis without process streams that cause forbidden matches, although the space limitation in the plant renders the best solution infeasible, the compromise solution can achieve a considerable steam saving of up to 8448 kW (91%) and is being evaluated by the plant managers and operators.

**Keywords:** energy conservation; heat integration; heat exchanger network synthesis; retrofit; mathematical programming; superstructure

### **1. Introduction**

To tackle the climate crisis and avoid the worst effects of climate change, many countries in the world have pledged to reduce global CO<sup>2</sup> emissions to net zero by 2050, while limiting the long-term increase in average global temperatures to 1.5 ◦C. Achieving the rapid reduction in CO<sup>2</sup> emissions in the net-zero emissions scenario (NZE) requires a range of technologies. Energy efficiency plays a key role in reducing energy consumption across the industry, buildings and transport sectors, and is among the key pillars of decarbonization of the global energy system, delivering around 50% of emissions reductions to 2030 in the NZE together with wind and solar [1].

In the industry sector, energy efficiency can be improved by exploiting process integration options such as waste heat recovery via heat exchange. The idea is to reduce the need for external heating and cooling using utilities (e.g., steam and cooling water) by transferring heat from process streams that need cooling to those that need heating in heat exchangers. However, synthesizing heat exchanger networks (HENs) using heuristics could lead to suboptimal solutions as a result of simply matching process streams near each other and losing sight of the big picture. Improving and optimizing the heat exchange between process streams involves systematic design procedures. Many techniques for heat integration and HEN synthesis have been developed since the 1970s [2], following two main approaches: pinch analysis [3] and mathematical programming [4]. Although pinch analysis allows a visualization of the process and provides insights into the energy use, mathematical programming would be preferred when addressing process constraints

**Citation:** Lee, J.-Y.; Chen, P.-Y. Optimization of Heat Recovery Networks for Energy Savings in Industrial Processes. *Processes* **2023**, *11*, 321. https://doi.org/10.3390/ pr11020321

Academic Editors: Pei Liu, Ming Liu and Xiao Wu

Received: 28 November 2022 Revised: 10 January 2023 Accepted: 14 January 2023 Published: 18 January 2023

**Copyright:** © 2023 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

(e.g., forbidden matches) and cost trade-offs. Arguably, all the problems that can be solved by pinch analysis can equally well be solved by mathematical programming, but not vice versa.

Recent works on HEN synthesis using mathematical programming have improved the quality of solutions [5–8] and considered retrofits [9,10], multiple plants [11–13] and multiperiod operations [14,15]. Huang and Karimi [5] proposed a multistage, matchcentric superstructure and a stage-less, exchanger-centric superstructure for HEN synthesis; the corresponding mixed integer nonlinear programming (MINLP) models consider a wider range of network configurations and allow solutions with lower total annualized costs (TACs) to be obtained. Beck and Hofmann [6] proposed the use of convex linear approximations to reformulate the original MINLP problem for HEN synthesis into a mixed integer linear programming (MILP) problem. Ziyatdinov et al. [7] introduced an iterative sequential method for optimal HEN synthesis. Liu et al. [8] presented an extended stage-wise superstructure with intermediate placement of utilities; the resulting complex mathematical model was solved using a genetic algorithm-based solution approach. For retrofit design, Cuˇcek et al. [ ˇ 9] presented approaches developed for the retrofit of existing HENs within individual processes and industrial sites to achieve energy savings, cost savings and emissions reductions. Pavão et al. [10] presented a superstructure-based model and a meta-heuristic solution approach for the retrofit of HENs. For multi-plant and multiperiod heat integration, Cheng et al. [11] developed a sequential design procedure to generate the inter-plant heat integration scheme in industrial parks. Nair et al. [12] proposed an MINLP model for configuring an eco-industrial park-wide multi-enterprise HEN with a practical and rational strategy. Kachacha et al. [13] used a MILP model for the synthesis of heat exchange and transport networks between processes considering the variability of streams. Pavão et al. [14] used a meta-heuristic two-level method to handle multiperiod HEN optimization with a post-optimization strategy to improve the results. Hofmann et al. [15] carried out the optimization of HENs and storage integration considering operational constraints of existing steam generation units.

Further extensions of HEN models include heat-integrated water integration [16,17], work integration [18–22] and process synthesis [23]. Ibri´c et al. [16] proposed a superstructure for the synthesis of heat-integrated water-using and wastewater treatment networks, and a two-step solution strategy to minimize TAC with upper bounds to freshwater and utility consumption. Ibri´c et al. [17] later addressed the synthesis of non-isothermal water networks thermally integrated with process streams. Huang and Karimi [18] presented an efficient superstructure-based MINLP formulation to synthesize work-heat exchanger networks. Li et al. [19] proposed a building block-based superstructure and formulated an MINLP optimization model for the synthesis of work and heat exchanger networks for minimum TAC. Viske et al. [20] proposed the use of a nonsmooth extension to handle unclassified process streams in the synthesis of work and heat exchange networks. Yang et al. [21] formulated a thermo-economic multi-objective MINLP model to synthesize sub- and above-ambient HENs of constant-pressure and pressure-changing streams with multi-stream compression. Hamedi et al. [22] presented a simulation-based method for work-integrated HEN synthesis without correlations or simplifying assumptions. Li et al. [23] proposed a surrogate-based optimization framework for the simultaneous synthesis of a chemical process and its HEN.

Most previous studies of HEN synthesis assume a global minimum temperature difference for non-isothermal process streams with constant heat capacities and do not consider isothermal process streams with phase change. In addition, most of the recently developed methods and models were applied only to examples and case studies in the literature, with limited practical applications. In this paper, a superstructure-based mathematical programming model for the optimization of HENs is developed and applied to a new industrial case study for minimizing utility consumption. The model handles process streams with phase change rigorously, as in the formulations of Ponce-Ortega et al. [24] and Watanapanich et al. [25], and can be modified to accommodate various practical considera-

tions by adjusting model parameters and adding appropriate constraints. The rest of the paper is organized as follows. A brief problem statement is given in Section 2. The model formulation is presented in Section 3, followed by the industrial case study in Section 4 to demonstrate its application. Conclusions and prospects for future work are given in Section 5. of the paper is organized as follows. A brief problem statement is given in Section 2. The model formulation is presented in Section 3, followed by the industrial case study in Section 4 to demonstrate its application. Conclusions and prospects for future work are given in Section 5. **2. Problem Statement**

and Watanapanich et al. [25], and can be modified to accommodate various practical considerations by adjusting model parameters and adding appropriate constraints. The rest

*Processes* **2023**, *11*, 321 3 of 17

#### **2. Problem Statement** Given a process with a set of hot streams ∈ **I** to be cooled and a set of cold streams

Given a process with a set of hot streams *i* ∈ **I** to be cooled and a set of cold streams *j* ∈ **J** to be heated, some of the process streams involve only sensible heat (temperature change), some involve only latent heat (phase change), and some involve both. The data (e.g., flowrates, inlet and outlet temperatures, heat capacities, latent heats and phase change fractions) for these process streams are known. The objective is to synthesize an optimal HEN of the process for minimum utility consumption, while satisfying the heating and cooling demands as well as practical constraints on HEN design. ∈ **J** to be heated, some of the process streams involve only sensible heat (temperature change), some involve only latent heat (phase change), and some involve both. The data (e.g., flowrates, inlet and outlet temperatures, heat capacities, latent heats and phase change fractions) for these process streams are known. The objective is to synthesize an optimal HEN of the process for minimum utility consumption, while satisfying the heating and cooling demands as well as practical constraints on HEN design.

### **3. Model Formulation 3. Model Formulation**

Figure 1 shows a condensed stage-wise HEN superstructure, which is divided into a number *K* of stages by *K* + 1 temperature locations. A hot process stream may be split and matched with cold process streams for heat exchange in each stage. If its cooling demand is not satisfied by heat exchange through the stages, a cold utility will be used at the outlet of the superstructure to provide the unsatisfied cooling. Likewise, a cold process stream may be split and matched with hot process streams in each stage. If its heating demand is not satisfied through the stages, a hot utility will be used at the outlet of the superstructure to provide the unsatisfied heating. Although only a match between one hot and one cold process stream is shown in Figure 1, the condensed superstructure can be expanded for the number of process streams. For split streams, isothermal mixing [4] is assumed to simplify the model formulation. It is also assumed that the heat losses from the system are small (with exchangers properly lagged) and can be neglected. Based on the superstructure in Figure 1, the HEN model is formulated and presented in the following subsections. Notation used is given in the Nomenclature in Appendix A. Figure 1 shows a condensed stage-wise HEN superstructure, which is divided into a number *K* of stages by *K* + 1 temperature locations. A hot process stream may be split and matched with cold process streams for heat exchange in each stage. If its cooling demand is not satisfied by heat exchange through the stages, a cold utility will be used at the outlet of the superstructure to provide the unsatisfied cooling. Likewise, a cold process stream may be split and matched with hot process streams in each stage. If its heating demand is not satisfied through the stages, a hot utility will be used at the outlet of the superstructure to provide the unsatisfied heating. Although only a match between one hot and one cold process stream is shown in Figure 1, the condensed superstructure can be expanded for the number of process streams. For split streams, isothermal mixing [4] is assumed to simplify the model formulation. It is also assumed that the heat losses from the system are small (with exchangers properly lagged) and can be neglected. Based on the superstructure in Figure 1, the HEN model is formulated and presented in the following subsections. Notation used is given in the Nomenclature in Appendix A.

**Figure 1.** Stage-wise superstructure for a HEN. **Figure 1.** Stage-wise superstructure for a HEN.

#### *3.1. Overall Energy Balance for Process Streams 3.1. Overall Energy Balance for Process Streams*

Equations (1) and (2) describe the overall energy balances for hot and cold process streams, respectively. Hot process stream *i* can be cooled by transferring heat to cold process streams *j* in stages *k* (*qijk*) and the cold utility ( cu) at the outlet of the superstructure. Equations (1) and (2) describe the overall energy balances for hot and cold process streams, respectively. Hot process stream *i* can be cooled by transferring heat to cold process streams *j* in stages *k* (*qijk*) and the cold utility (*q* cu *i* ) at the outlet of the superstructure. Cold process stream *j* can be heated by transferring heat from hot process streams *i* in stages *k* (*qijk*) and the hot utility (*q* hu *j* ) at the outlet of the superstructure.

$$F\_i\left(\mathbf{C}p\_i\left(T\_i^{\text{in}} - T\_i^{\text{out}}\right) + \lambda\_i a\_i\right) = \sum\_{j \in \mathbf{J}} \sum\_{k \in \mathbf{ST}} q\_{ijk} + q\_i^{\text{cu}} \quad \forall i \in \mathbf{I} \tag{1}$$

$$F\_{\dot{\boldsymbol{\beta}}}\left(\mathbf{C}\boldsymbol{p}\_{\dot{\boldsymbol{\beta}}}\left(\boldsymbol{T}\_{\dot{\boldsymbol{\beta}}}^{\rm out} - \boldsymbol{T}\_{\dot{\boldsymbol{\beta}}}^{\rm in}\right) + \lambda\_{\dot{\boldsymbol{\beta}}}\boldsymbol{a}\_{\dot{\boldsymbol{\beta}}}\right) = \sum\_{\boldsymbol{i} \in \mathbf{I}} \sum\_{k \in \mathbf{ST}} q\_{\dot{\boldsymbol{i}}\boldsymbol{j}\boldsymbol{k}} + q\_{\dot{\boldsymbol{j}}}^{\rm hu} \quad \forall \boldsymbol{j} \in \mathbf{J} \tag{2}$$

where *F<sup>i</sup>* and *F<sup>j</sup>* are the flowrates, *Cp<sup>i</sup>* and *Cp<sup>j</sup>* the heat capacities, *λ<sup>i</sup>* and *λ<sup>j</sup>* the latent heats, and *α<sup>i</sup>* and *α<sup>j</sup>* the phase change fractions of hot process stream *i* and cold process stream *j*, respectively. For streams involving only sensible heat, the latent heat and phase change fraction are set to zero; for streams involving only latent heat, the heat capacity is set to zero. Thus, Equations (1) and (2) apply to three types of process streams, including those involving both sensible and latent heat.

### *3.2. Energy Balance for Process Streams in Each Stage*

Equations (3) and (4) describe the energy balances in stage *k* for hot (*i* ∈ **I** <sup>1</sup> <sup>∪</sup> **<sup>I</sup>** 3 ) and cold process streams involving sensible heat (*j* ∈ **J** <sup>1</sup> <sup>∪</sup> **<sup>J</sup>** 3 ), respectively. Process streams involving only latent heat (*i* ∈ **I** 2 and *j* ∈ **J** 2 ) are isothermal and can be modelled simply by Equations (1) and (2).

$$F\_i \mathbb{C} \, p\_i(t\_{ik} - t\_{i,k+1}) + q\_{ik}^\Lambda = \sum\_{j \in \mathcal{J}} q\_{ijk} \quad \forall i \in \mathbf{I}^1 \cup \mathbf{I}^3, k \in \mathbf{ST} \tag{3}$$

$$F\_{\hat{\jmath}} \mathbb{C} p\_{\hat{\jmath}} \left( t\_{jk} - t\_{j,k+1} \right) + q\_{jk}^{\Lambda} = \sum\_{i \in \mathbf{I}} q\_{ijk} \quad \forall j \in \mathbf{J}^1 \cup \mathbf{J}^3, k \in \mathbf{ST} \tag{4}$$

where *tik* is the temperature of hot process stream *i* at temperature location *k*, *tjk* the temperature of cold process stream *j* at temperature location *k*, *q* Λ *ik* the condensation heat load of hot process stream *i* in stage *k*, and *q* Λ *jk* the evaporation heat load of cold process stream *j* in stage *k*. For streams involving only sensible heat (*i* ∈ **I** 1 and *j* ∈ **J** 1 ), the latent heat loads (*q* Λ *ik* and *q* Λ *jk*) are set to zero.

### *3.3. Temperature Assignment and Feasibility Constraints*

As shown in Figure 1, hot process streams enter the superstructure from the lefthand side, while cold process streams enter from the right-hand side. Therefore, the inlet temperature of hot process stream *i* (*T* in *i* ) is equal to its temperature at the first temperature location (*ti*,1), and the inlet temperature of cold process stream *j* (*T* in *j* ) is equal to its temperature at the last temperature location (*tj*,*K*+1), as given in Equations (5) and (6). Equations (7)–(10) specify the trends in temperature. For non-isothermal process streams, a monotonic decrease in temperature throughout the superstructure is expected; for isothermal process streams, the temperature remains constant at the inlet temperature. In addition, the outlet temperature of hot process stream *i* (*T* out *i* ) is used as a lower limit for its temperatures in the superstructure, while the outlet temperature of cold process stream *j* (*T* out *j* ) is used as an upper limit, as given in Equations (11) and (12).

$$t\_{i,1} = T\_i^{\text{in}} \; \forall i \in \mathbf{I} \tag{5}$$

$$t\_{j,K+1} = T\_j^{\text{in}} \quad \forall j \in \mathbf{J} \tag{6}$$

$$t\_{ik} \ge t\_{i,k+1} \quad \forall i \in \mathbf{I}^1 \cup \mathbf{I}^3, k \in \mathbf{ST} \tag{7}$$

$$t\_{ik} = T\_i^{\text{in}} \quad \forall i \in \mathbf{I}^2, k \in \mathbf{ST} \tag{8}$$

$$t\_{jk} \ge t\_{j,k+1} \quad \forall j \in \mathbf{J}^1 \cup \mathbf{J}^3, k \in \mathbf{ST} \tag{9}$$

$$t\_{jk} = T\_j^{\text{in}} \quad \forall j \in \mathbf{J}^2, k \in \mathbf{ST} \tag{10}$$

*tik* ≥ *T* out *i* ∀*i* ∈ **I**, *k* ∈ **ST** (11)

$$t\_{jk} \le T\_j^{\text{out}} \quad \forall j \in \mathbf{J}, k \in \mathbf{ST} \tag{12}$$

### *3.4. Utility Requirements and Latent Heat Balance*

Equations (13) and (14) are used to calculate the hot and cold utility requirements for non-isothermal process streams.

$$q\_i^{\mathbf{cu}} = F\_i \mathbb{C} p\_i \left( t\_{i, \mathbf{K} + 1} - T\_i^{\text{out}} \right) + q\_i^{\Lambda, \text{cu}} \; \forall i \in \mathbf{I}^1 \cup \mathbf{I}^3 \tag{13}$$

$$q\_j^{\text{hu}} = F\_j \mathbb{C} p\_j \left( T\_j^{\text{out}} - t\_{j,1} \right) + q\_j^{\Lambda, \text{hu}} \quad \forall j \in \mathbf{J}^1 \cup \mathbf{J}^3 \tag{14}$$

where *q* Λ,cu *i* is the condensation heat load of hot process stream *i* in the utility cooler and *q* Λ,hu *j* is the evaporation heat load of cold process stream *j* in the utility heater. Again, for streams involving only sensible heat, the latent heat loads are set to zero.

Equations (15) and (16) describe the latent heat balances for process streams involving both sensible and latent heat (*i* ∈ **I** 3 and *j* ∈ **J** 3 ).

$$F\_{\bar{i}}\lambda\_{\bar{i}}\alpha\_{\bar{i}} = \sum\_{k \in \mathbf{ST}} q\_{\bar{i}k}^{\Lambda} + q\_{\bar{i}}^{\Lambda, \text{cu}} \quad \forall i \in \mathbf{I}^3 \tag{15}$$

$$F\_{\dot{\jmath}}\lambda\_{\dot{\jmath}}\alpha\_{\dot{\jmath}} = \sum\_{k \in \mathbf{ST}} q\_{\dot{\jmath}k}^{\Lambda} + q\_{\dot{\jmath}}^{\Lambda, \text{hu}} \quad \forall \dot{\jmath} \in \mathbf{J}^3 \tag{16}$$

### *3.5. Feasibility of Latent Heat Exchange*

Equations (17)–(20) ensure feasible latent heat exchange for type-3 process streams. Equations (17) and (18) are used for hot process streams that are condensed at the inlet (supply) temperature and cooled to the outlet (target) temperature. If hot process stream *i* enters stage *k* at a temperature lower than the condensation temperature (*tik* < *T* in *i* ), binary variable *yik* has to be zero, and *q* Λ *ik* is forced to zero, meaning no condensation of hot process stream *i* in stage *k*. Equations (19) and (20) are used for cold process streams that are heated from the inlet temperature to the outlet temperature and evaporated. If cold process stream *j* leaves stage *k* at a temperature lower than the evaporation temperature (*tjk* < *T* out *j* ), binary variable *yjk* has to be zero, and *q* Λ *jk* is forced to zero, meaning no evaporation of cold process stream *j* in stage *k*.

$$
\varepsilon - \Gamma y\_{ik} \le T\_i^{\text{in}} - t\_{ik} \le \Gamma (1 - y\_{ik}) \quad \forall i \in \Gamma^3, k \in \mathbf{ST} \tag{17}
$$

$$q\_{ik}^{\Lambda} \le Q\_i^{\mathbf{U}} y\_{ik} \quad \forall i \in \mathbf{I}^3, k \in \mathbf{ST} \tag{18}$$

$$
\varepsilon-\Gamma y\_{jk}\leq T\_j^{\text{out}}-t\_{jk}\leq \Gamma\left(1-y\_{jk}\right) \quad \forall j\in \mathbf{J}^3, k\in \mathbf{ST} \tag{19}
$$

$$q\_{jk}^{\Lambda} \le Q\_j^{\mathbf{U}} y\_{jk} \quad \forall j \in \mathbf{J}^3, k \in \mathbf{ST} \tag{20}$$

where *ε* is a sufficiently small value, Γ a sufficiently large value, *yik* indicates if *tik* ≥ *T* in *i* , and *yjk* indicates if *tjk* ≥ *T* out *j* .

### *3.6. Logical Constraints and Temperature Differences*

Equations (21)–(23) set upper limits on the heat loads of exchangers (*Q*<sup>U</sup> *ij*), coolers (*Q*<sup>U</sup> *i* ) and heaters (*Q*<sup>U</sup> *j* ), and force the heat load to zero if the match does not exist.

$$q\_{\rm ijk} \le Q\_{ij}^{\rm U} z\_{\rm ijk} \quad \forall i \in \mathbf{I}, j \in \mathbf{J}, k \in \mathbf{ST} \tag{21}$$

$$\mathbf{q}\_{i}^{\mathbf{cu}} \le \mathbf{Q}\_{i}^{\mathbf{U}} \mathbf{z}\_{i}^{\mathbf{cu}} \quad \forall i \in \mathbf{I} \tag{22}$$

$$\mathbf{q}\_{j}^{\text{hu}} \le \mathbf{Q}\_{j}^{\text{U}} z\_{j}^{\text{hu}} \quad \forall i \in \mathbf{I} \tag{23}$$

where *zijk*, *z* cu *i* and *z* hu *j* are binary variables indicating the existence of matches between hot process stream *i* and cold process stream *j* in stage *k*, between hot process stream *i* and a cold utility, and between a hot utility and cold process stream *j*, respectively.

Equations (24)–(27) ensure feasible driving forces for existing matches. Different minimum temperature differences can be assigned to different matches.

$$
\Delta T\_{ij}^{\min} - \Gamma \left( 1 - z\_{ijk} \right) \le t\_{ik} - t\_{jk} \quad \forall i \in \mathbf{I}, j \in \mathbf{J}, k \in \mathbf{ST} \tag{24}
$$

$$
\Delta T\_{ij}^{\min} - \Gamma \left( 1 - z\_{ijk} \right) \le t\_{i,k+1} - t\_{j,k+1} \quad \forall i \in \mathbf{I}, j \in \mathbf{J}, k \in \mathbf{ST} \tag{25}
$$

$$
\Delta T\_{i, \text{cu}}^{\text{min}} - \Gamma (1 - z\_i^{\text{cu}}) \le t\_{i, \text{K}+1} - T\_{\text{cu}}^{\text{out}} \quad \forall i \in \mathbf{I} \tag{26}
$$

$$
\Delta T\_{\rm hu,j}^{\rm min} - \Gamma \left( 1 - z\_j^{\rm hu} \right) \le T\_{\rm hu}^{\rm out} - t\_{j,1} \quad \forall j \in \mathbf{J} \tag{27}
$$

where ∆*T* min *ij* , ∆*T* min *<sup>i</sup>*,cu and ∆*T* min hu,*j* are the minimum temperature differences for matches between hot process stream *i* and cold process stream *j*, between hot process stream *i* and the cold utility, and between the hot utility and cold process stream *j*, respectively; *T* out cu and *T* out hu are the outlet temperatures of the cold and hot utilities, respectively. Equations (26) and (27) may not be needed assuming the hot utility to be sufficiently hot and the cold utility sufficiently cold. To prevent temperature inversion in heat exchangers involving condensation or evaporation of process streams, additional constraints can be used [25]. Alternatively, the minimum temperature difference can be manipulated (i.e., increased) for specific matches.

### *3.7. Objective Function*

For maximum heat recovery in the process, the objective function is to minimize the hot utility requirement (*f* HU), as given in Equation (28).

$$\min f\_{\text{HU}} = \sum\_{j \in \mathbf{J}} q\_j^{\text{hu}} \tag{28}$$

A further step may be taken to obtain simpler HEN designs and keep down the capital cost by minimizing the number of matches or units (*f* NU), as given in Equation (29), subject to the minimum hot utility requirement (*f* ∗ HU) in Equation (30).

$$\min f\_{\text{NU}} = \sum\_{i \in \mathcal{I}} \sum\_{j \in \mathcal{J}} \sum\_{k \in \mathbf{ST}} z\_{ijk} + \sum\_{i \in \mathcal{I}} z\_i^{\text{cu}} + \sum\_{j \in \mathcal{J}} z\_j^{\text{hu}} \tag{29}$$

$$\sum\_{j \in \mathcal{J}} q\_j^{\text{hu}} \le f\_{\text{HU}}^\* \tag{30}$$

Equations (1)–(30) constitute a MILP model, for which global optimality can be guaranteed. The model is solved to determine the optimal HEN. The main results are stream temperatures and heat duties of the matches between process and utility streams. These details can then be used in the economic assessment involving the selection of heat transfer equipment and the calculation of heat transfer area and exchanger cost.

In the next section, an industrial case study is presented to illustrate the application of the proposed MILP model, which is solved in GAMS using solver CPLEX. Solutions were found with reasonable processing time. It is worth noting that numerical uncertainty comes from numerical errors and approximations in the implementation of the mathematical model. In this work, the accuracy of the proposed model was confirmed by comparing the base-case results with the actual conditions. Numerical errors are small enough to be neglected.

### **4. Industrial Case Study**

This case study considers an ethylene glycol (EG) production plant located in southern Taiwan. The EG process involves a reactor, two absorbers, two strippers, an eightstage evaporator, a drying tower and a distillation column followed by three separation columns [26]. Tables 1 and 2 show the stream data extracted from the material and energy balance. The EG process flowsheet is complicated and is not shown for confidentiality reasons. Two levels of saturated steam at 213 ◦C and 180 ◦C are available as hot utilities. Cooling water as the cold utility is available at 29 ◦C with a return temperature of 39 ◦C.



MEG, monoethylene glycol; DEG, diethylene glycol.

### **Table 2.** Cold process stream data.


EO, ethylene oxide.

It should be noted that the HEN model presented in Section 3 does not consider multiple hot or cold utilities. For two levels of steam, this model is solved to determine the minimum total steam demand first. The loads on the two steam levels are then determined according to the outlet temperatures of cold streams that require steam.

### *4.1. Preliminary Analysis*

Currently, the EG process has a HEN with nine matches between H1–H9 and C1–C6, as shown in Table 3, where the heat loads were determined by solving the HEN model considering only the existing matches. More details (e.g., intermediate temperatures of the process streams) are given in Table S1. These heat loads and stream temperatures agree with the operating conditions in the plant, thus validating the proposed model. The condensers

(H10–H14), evaporator (C7) and reboilers (C8-C12) were not included for heat integration. With the current HEN, the process requires 40,567 kW of 213 ◦C steam, 14,920 kW of 180 ◦C steam and 56,325 kW of cooling water. This is taken as the base case.

**Table 3.** Existing HEN of the EG process.


Heat loads of the matches are given in kW.

To retrofit the existing HEN for improved heat recovery and reduced steam requirements, all possible matches are considered. Table 4 shows the minimum permissible temperature differences for matches between process and utility streams.


**Table 4.** Minimum temperature differences (◦C) for different matches.

For process–process matches, a minimum temperature difference of 10 ◦C is assumed. However, a smaller value is used if the existing match has a smaller temperature difference. For such matches (H2–C2, H4–C3, H5–C3 and H9–C6), the minimum temperature difference is set to the smaller of the temperature differences at the cold and hot ends of the exchanger. For process–utility matches, the minimum temperature differences are set according to the current operating conditions or based on practical considerations [26]. To avoid a complex HEN retrofit, the following assumptions are made:


These constraints are given in Equations (31)–(35).

$$\sum\_{j \in \mathbf{J}} z\_{ijk} \le 1 \quad \forall i \in \mathbf{I}, k \in \mathbf{ST} \tag{31}$$

$$\sum\_{i \in \mathbf{I}} z\_{ijk} \le 1 \quad \forall j \in \mathbf{J}, k \in \mathbf{ST} \tag{32}$$

$$\sum\_{k \in \mathbf{ST}} z\_{ijk} \le 1 \quad \forall i \in \mathbf{I}, j \in \mathbf{J} \tag{33}$$

$$\sum\_{j \in \mathbf{J}} \sum\_{k \in \mathbf{ST}} z\_{ijk} \le 1 \quad \forall i \in \{\mathbf{H}10\text{-H14}\} \tag{34}$$

$$\sum\_{i \in \mathbf{I}} \sum\_{k \in \mathbf{ST}} z\_{ijk} \le 1 \quad \forall j \in \{\mathbf{C7-C12}\} \tag{35}$$

Solving the HEN model considering all possible matches in seven stages yields the results in Table 5. More details of this solution are given in Table S2. The minimum utility requirements are determined to be 13,710 kW of 213 ◦C steam, 26,302 kW of 180 ◦C steam and 40,849 kW of cooling water. This can be achieved in the minimum number of units of 33 (20 exchangers, five heaters and eight coolers). The main feature of this optimized HEN is that the outlet recycle gas (H1) is used to heat the evaporator (C7) and reboiler streams C8 and C9. In addition, condenser streams H11-H14 are used to heat the inlet recycle gas (C1) and recycled water (C2). Compared with the base case, the total steam demand decreases by 28% (despite a 76% increase in the load of 180 ◦C steam), while the cooling water demand decreases by 27%.

**Table 5.** Optimized HEN of the EG process.


Heat loads of the matches are given in kW.

Figure 2 shows the variation of minimum utility requirements with the number of stages in the superstructure. Increasing the number of stages does not result in significant reductions in steam and cooling water demands. The optimized HENs have the common feature that H1 is used to heat C7–C9. Therefore, the design with the fewest (33) heat exchange units (i.e., the one in Table 5) was selected for further evaluation. To achieve the theoretical reduction in total steam demand of 15,476 kW, the HEN retrofit entails 14 new matches and two expanded matches. Making all these changes to the current HEN might be prohibitively expensive. Moreover, it was found that several of the new matches (H1–C7, H1–C8, H1–C9, H6–C1, H13–C1 and H14–C1) could cause EO penetration, which

will adversely affect the MEG product. Such matches are deemed risky and should be forbidden, thus rendering the solution infeasible. It was then suggested that the condensers, evaporator and reboilers should be excluded from the heat integration, despite the large potential steam saving. den, thus rendering the solution infeasible. It was then suggested that the condensers, evaporator and reboilers should be excluded from the heat integration, despite the large potential steam saving.

Figure 2 shows the variation of minimum utility requirements with the number of stages in the superstructure. Increasing the number of stages does not result in significant reductions in steam and cooling water demands. The optimized HENs have the common feature that H1 is used to heat C7–C9. Therefore, the design with the fewest (33) heat exchange units (i.e., the one in Table 5) was selected for further evaluation. To achieve the theoretical reduction in total steam demand of 15,476 kW, the HEN retrofit entails 14 new matches and two expanded matches. Making all these changes to the current HEN might be prohibitively expensive. Moreover, it was found that several of the new matches (H1– C7, H1–C8, H1–C9, H6–C1, H13–C1 and H14–C1) could cause EO penetration, which will adversely affect the MEG product. Such matches are deemed risky and should be forbid-

**Figure 2.** Minimum utility requirements versus the number of stages. **Figure 2.** Minimum utility requirements versus the number of stages.

### *4.2. Analysis without Condensers, Evaporator and Reboilers*

*Processes* **2023**, *11*, 321 10 of 17

*4.2. Analysis without Condensers, Evaporator and Reboilers* Having excluded H10–H14 and C7–C12 from the analysis, the process stream data are updated as in Table 6. It was found that the heat capacity of the recycle gas varies significantly for large temperature changes and cannot be taken as a constant, as shown in Figure 3. The variation in heat capacity of recycle gas streams (*Cp*†) with temperature Having excluded H10–H14 and C7–C12 from the analysis, the process stream data are updated as in Table 6. It was found that the heat capacity of the recycle gas varies significantly for large temperature changes and cannot be taken as a constant, as shown in Figure 3. The variation in heat capacity of recycle gas streams (*Cp*†) with temperature (*T*) can be correlated as:

$$\mathbf{C}p\_{\dagger} = a\_{0,\dagger} + a\_{1,\dagger}T \quad \dagger \in \{\text{H1,C1}\} \tag{36}$$

 = , + , †∈ {H1,C1} (36) where *a*0,† and *a*1,† are constants determined by correlating the heat capacity data: *a*0,H1 = 0.396767 and *a*1,H1 = 0.000578 for 60 °C ≤ *T* ≤ 222 °C (*R*<sup>2</sup> = 0.999); *a*0,C1 = 0.399877 and *a*1,C1 = where *a*0,† and *a*1,† are constants determined by correlating the heat capacity data: *a*0,H1 = 0.396767 and *a*1,H1 = 0.000578 for 60 ◦C ≤ *T* ≤ 222 ◦C (*R* <sup>2</sup> = 0.999); *a*0,C1 = 0.399877 and *a*1,C1 = 0.000555 for 45 ◦C ≤ *T* ≤ 193 ◦C (*R* <sup>2</sup> = 0.997). Thus,

$$
\Delta H\_{\rm \dagger} = \int\_{T\_1}^{T\_2} \mathbb{C} p\_{\rm \dagger} dT = \left[ a\_{0, \rm \dagger} T + \frac{a\_{1, \rm \dagger} T^2}{2} \right]\_{T\_1}^{T\_2} \dagger \in \{\text{H1, C1}\} \tag{37}
$$

 where ∆*H*† is the enthalpy change from *T*<sup>1</sup> to *T*2.

(*T*) can be correlated as:

where ∆*H*† is the enthalpy change from *T*<sup>1</sup> to *T*2. Equations (1)–(4), (13) and (14) are then modified for H1 and C1 by replacing "*Cp*∆*T*" with "∆*H*," which involves nonlinear terms in Equations (3), (4), (13) and (14). Thus, the HEN model becomes an MINLP, which is solved using solver BARON in GAMS. Some of the process stream data have also been revised in response to the change in heat capac-Equations (1)–(4), (13) and (14) are then modified for H1 and C1 by replacing "*Cp*∆*T*" with "∆*H*," which involves nonlinear terms in Equations (3), (4), (13) and (14). Thus, the HEN model becomes an MINLP, which is solved using solver BARON in GAMS. Some of the process stream data have also been revised in response to the change in heat capacity of the recycle gas.

ity of the recycle gas. Figure 4 shows the revised base-case conditions. The details were determined by solving the HEN model considering only the existing matches. The resulting heat loads and stream temperatures agree with the updated operating conditions in the plant. Excluding the utility requirements of H10–H14 and C7–C12, the process requires 9277 kW of 180 ◦C steam and 46,089 kW of cooling water. For improved heat recovery and reduced steam demand, the existing HEN is optimized considering all feasible matches. Forbidden matches such as H6-C1, H8-C1 and H9-C1 can be included using Equation (38).

$$z\_{i\bar{j}k} \le M\_{\bar{i}\bar{j}} \quad \forall i \in \mathbf{I}, j \in \mathbf{J}, k \in \mathbf{ST} \tag{38}$$

where *Mij* is a binary parameter denoting allowable (*Mij* = 1) and forbidden matches (*Mij* = 0). To avoid a complex HEN retrofit, the following assumptions hold:

• No stream splitting is allowed (Equations (31) and (32)); **Inlet Outlet Heat Latent**

**(°C)**

• Any pair of hot and cold process streams has at most one match (Equation (33)). **Stream Description Flowrate Temperature Temperature Capacity (J/g/°C) Heat Phase Change**

**(J/g)**

**Fraction**


**Table 6.** Updated process stream data. **(°C)** H1 Outlet recycle gas 981 222 60 *f*H1(*T*) - 0

**Table 6.** Updated process stream data.

**(t/h)**

Figure 4 shows the revised base-case conditions. The details were determined by solving the HEN model considering only the existing matches. The resulting heat loads In addition, to make use of the existing matches and keep down the number of network modifications, Equation (39) is used.

$$\begin{array}{ll} \sum z\_{i\parallel\overline{k}} = 1 \quad \forall (i, j) \\ \sum \{ (\text{H}1, \text{C}1), (\text{H}1, \text{C}2), (\text{H}2, \text{C}2), (\text{H}3, \text{C}4), (\text{H}4, \text{C}3), (\text{H}5, \text{C}3), (\text{H}6, \text{C}3), (\text{H}8, \text{C}5), (\text{H}9, \text{C}6) \} \end{array} \tag{39}$$

steam demand, the existing HEN is optimized considering all feasible matches. Forbidden matches such as H6-C1, H8-C1 and H9-C1 can be included using Equation (38). This constraint enables the existing matches to be preferred when determining the minimum number of units, if the same steam saving can be achieved.

 ≤ ∀ ∈ **I**, ∈ **J**, ∈ **ST** (38) where *Mij* is a binary parameter denoting allowable (*Mij* = 1) and forbidden matches (*Mij* = 0). To avoid a complex HEN retrofit, the following assumptions hold: No stream splitting is allowed (Equations (31) and (32)); Any pair of hot and cold process streams has at most one match (Equation (33)). Solving the HEN model considering all feasible matches in six stages yields the results in Figure 5. The minimum utility requirements are determined to be 0 kW of 180 ◦C steam and 36,812 kW of cooling water. This can be achieved in the minimum number of units of 15 (10 exchangers, no heater and five coolers). The main feature of this optimized HEN is that instead of using steam, the outlet recycle gas (H1) is used to heat the EO solution (C3). Compared with the base case in Figure 4, the steam demand decreases by 100%, while the cooling water demand decreases by 20%. The cost benefit associated with this steam saving is calculated to be 7.5 million USD/y, assuming steam costs 0.101 USD/kWh and the plant operates for 8000 h/y [26].

∈ {(H1,C1), (H1,C2), (H2,C2), (H3,C4), (H4,C3), (H5,C3), (H6,C3), (H8,C5), (H9,C6)}

minimum number of units, if the same steam saving can be achieved.

work modifications, Equation (39) is used.

In addition, to make use of the existing matches and keep down the number of net-

This constraint enables the existing matches to be preferred when determining the

then suggested seeking an alternative solution that does not require additional heat trans-

(39)

**Figure 4.** Existing HEN of the EG process (updated). **Figure 4.** Existing HEN of the EG process (updated). fer area for match H1–C1.

 ∈**ST**

= 1 ∀(,)

**Figure 5.** Optimized HEN of the EG process (without condensers, evaporator and reboilers). **Figure 5.** Optimized HEN of the EG process (without condensers, evaporator and reboilers).

To avoid additional heat transfer area required for match H1–C1, the minimum temperature difference between H1 and C1 is set to 29.3 °C (instead of 10 °C) according to the base-case conditions. Additionally, it is desirable to have no more than three hot process streams matched with C1, if its heating demand cannot be satisfied solely by H1. This constraint is given in Equation (40). ,C1, ∈I ∈ST ≤ 3 (40) Solving the HEN model with the revised ∆H1,C1 min and the constraint on the number of matches of C1 gives the results in Figure 6. The minimum utility requirements are determined to be 829 kW of 180 °C steam and 37,641 kW of cooling water. This can be achieved in the minimum number of units of 18 (11 exchangers, one heater and six coolers). In this HEN design, H1 is used to heat C3 instead of using H5, which is used to heat C1 together with H4 and H1. Due to practical constraints, the steam saving decreases from 100 to 91%, while the cooling water demand increases by 2% (from 36,812 to 37,641 kW). It should be noted that match H5–C1 is a gas–gas match, for which the current minimum This solution was further evaluated for its feasibility. To achieve the theoretical reduction in steam demand of 9277 kW, the HEN retrofit entails adding a new match (H1–C3) and modifying two existing ones for decreased temperature differences (H1–C1) and increased heat transfer (H2–C2). The total retrofit cost is estimated at 24.2 million USD by the process engineers. Dividing the total retrofit cost by the annual cost saving gives a payback period of 3.2 y. On the other hand, the decrease in the logarithmic mean temperature difference of match H1–C1 (from 32 to 13.3 ◦C) requires a three-times larger heat transfer area for the same heat duty. Although the design in Figure 5 does not pose a hazard to the EG process, the plant has no space to accommodate the required area increase. Therefore, such a retrofit project cannot be implemented. The process engineers then suggested seeking an alternative solution that does not require additional heat transfer area for match H1–C1.

area. Therefore, an additional case is considered where ∆H5,C1

temperature difference of 10 °C may be too small and may lead to a large heat transfer

°C to avoid a too large area requirement. Solving the HEN model with the revised ∆H5,C1

gives the results in Figure S1. With the same HEN structure as in Figure 6 but changes in some heat loads and temperatures, the steam saving decreases from 91 to 50%, while the

cates the trade-off between energy and capital. So far, both solutions in Figures 6 and S1

min is revised upwards to 20

min

To avoid additional heat transfer area required for match H1–C1, the minimum temperature difference between H1 and C1 is set to 29.3 ◦C (instead of 10 ◦C) according to the base-case conditions. Additionally, it is desirable to have no more than three hot process streams matched with C1, if its heating demand cannot be satisfied solely by H1. This constraint is given in Equation (40).

$$\sum\_{i \in \mathcal{I}} \sum\_{k \in \mathcal{ST}} z\_{i, \mathcal{C}1, k} \le 3 \tag{40}$$

Solving the HEN model with the revised ∆*T* min H1,C1 and the constraint on the number of matches of C1 gives the results in Figure 6. The minimum utility requirements are determined to be 829 kW of 180 ◦C steam and 37,641 kW of cooling water. This can be achieved in the minimum number of units of 18 (11 exchangers, one heater and six coolers). In this HEN design, H1 is used to heat C3 instead of using H5, which is used to heat C1 together with H4 and H1. Due to practical constraints, the steam saving decreases from 100 to 91%, while the cooling water demand increases by 2% (from 36,812 to 37,641 kW). It should be noted that match H5–C1 is a gas–gas match, for which the current minimum temperature difference of 10 ◦C may be too small and may lead to a large heat transfer area. Therefore, an additional case is considered where ∆*T* min H5,C1 is revised upwards to 20 ◦C to avoid a too large area requirement. Solving the HEN model with the revised ∆*T* min H5,C1 gives the results in Figure S1. With the same HEN structure as in Figure 6 but changes in some heat loads and temperatures, the steam saving decreases from 91 to 50%, while the cooling water demand increases from 37,641 to 41,441 kW. The results comparison indicates the trade-off between energy and capital. So far, both solutions in Figure 6 and Figure S1 have been taken as alternatives being evaluated by the process engineers for technical and economic feasibility. Although the complete economic assessments are not available yet, the payback period is likely to be shorter than 3.2 y because modifications for match H1-C1 are no longer needed. The capital cost of modifying match H1–C1 (including the equipment and auxiliaries) is 20.1 million USD, accounting for 83% of the total retrofit cost for the design in Figure 5. *Processes* **2023**, *11*, 321 14 of 17 have been taken as alternatives being evaluated by the process engineers for technical and economic feasibility. Although the complete economic assessments are not available yet, the payback period is likely to be shorter than 3.2 y because modifications for match H1- C1 are no longer needed. The capital cost of modifying match H1–C1 (including the equipment and auxiliaries) is 20.1 million USD, accounting for 83% of the total retrofit cost for the design in Figure 5.

**Figure 6.** Alternative HEN of the EG process. **Figure 6.** Alternative HEN of the EG process.

**5. Conclusions**

### **5. Conclusions**

A superstructure-based mathematical optimization model for the synthesis of HENs has been presented in this paper. This HEN model considers three types of process streams involving temperature and phase changes, and has been applied to an industrial case study of optimizing the heat recovery system of an EG process. The results obtained show a compromise between steam savings and process feasibility, and how the proposed model can be modified to reflect practical considerations such as stream matching constraints and variable heat capacities as a function of temperature. Although in the prelim-A superstructure-based mathematical optimization model for the synthesis of HENs has been presented in this paper. This HEN model considers three types of process streams involving temperature and phase changes, and has been applied to an industrial case study of optimizing the heat recovery system of an EG process. The results obtained show a compromise between steam savings and process feasibility, and how the proposed model

matches and temperatures of the optimized HEN (preliminary analysis).

inary analysis a substantial steam saving of 15,476 kW (28%) can be achieved, which is equivalent to a cost saving of 12.5 million USD/y, the solution involves forbidden matches that pose a hazard to the MEG product and is not going to be implemented. In the analysis

are currently being evaluated in more detail. The corresponding cost saving could double if carbon taxes are considered. Future work will extend the optimization approach to other EG plants on the same site for improved heat recovery. The focus will also be on the optimization of utility systems, particularly those with time-dependent supply and demand,

**Supplementary Materials:** The following supporting information can be downloaded at: https://www.mdpi.com/article/10.3390/pr11020321/s1, Figure S1: alternative HEN of the EG process

min = 20℃); Table S1: stream matches and temperatures of the existing HEN; Table S2: stream

(∆H5,C1

for which multiperiod planning is indicated.

can be modified to reflect practical considerations such as stream matching constraints and variable heat capacities as a function of temperature. Although in the preliminary analysis a substantial steam saving of 15,476 kW (28%) can be achieved, which is equivalent to a cost saving of 12.5 million USD/y, the solution involves forbidden matches that pose a hazard to the MEG product and is not going to be implemented. In the analysis without the condensers, evaporator and reboilers, despite the space limitation in the plant, the achievable steam saving of up to 8448 kW (91%) is still considerable and the solutions are currently being evaluated in more detail. The corresponding cost saving could double if carbon taxes are considered. Future work will extend the optimization approach to other EG plants on the same site for improved heat recovery. The focus will also be on the optimization of utility systems, particularly those with time-dependent supply and demand, for which multiperiod planning is indicated.

**Supplementary Materials:** The following supporting information can be downloaded at: https://www.mdpi.com/article/10.3390/pr11020321/s1, Figure S1: alternative HEN of the EG process (∆*T* min H5,C1 = 20 ◦C); Table S1: stream matches and temperatures of the existing HEN; Table S2: stream matches and temperatures of the optimized HEN (preliminary analysis).

**Author Contributions:** Conceptualization, J.-Y.L.; methodology, J.-Y.L.; software, J.-Y.L. and P.-Y.C.; validation, J.-Y.L. and P.-Y.C.; formal analysis, J.-Y.L.; investigation, J.-Y.L.; resources, J.-Y.L.; data curation, J.-Y.L.; writing—original draft preparation, J.-Y.L.; writing—review and editing, J.-Y.L. and P.-Y.C.; visualization, J.-Y.L. and P.-Y.C.; supervision, J.-Y.L.; project administration, J.-Y.L.; funding acquisition, J.-Y.L. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by the National Science and Technology Council (NSTC) of R.O.C., grant number 111-2221-E-027-003.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** The data presented in this study are available on request from the corresponding author.

**Acknowledgments:** The authors thank Nan Ya Plastics Corporation for financial support from an industry– academia collaboration project (grant number 211A058) with National Taipei University of Technology.

**Conflicts of Interest:** The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript; or in the decision to publish the results.

### **Appendix A. Nomenclature**

Notation used in the HEN model is given below.

### *Appendix A.1. Indices and Sets*


### *Appendix A.2. Parameters*


### *Appendix A.3. Variables*


### **References**


**Disclaimer/Publisher's Note:** The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.
