CHEMICAL ENGINEERING TRANSACTIONS VOL. 76, 2019 A publication of The Italian Association of Chemical Engineering Online at www.aidic.it/cet Guest Editors: Petar S. Varbanov, Timothy G. Walmsley, Jiří J. Klemeš, Panos Seferlis Copyright © 2019, AIDIC Servizi S.r.l. ISBN 978-88-95608-73-0; ISSN 2283-9216 Mathematical Programming for Heat Exchanger Network Design with Temperature-Dependent Specific Heat Capacity Siwat Valeekiatkul, Kitipat Siemanond* The Petroleum and Petrochemical College, Chulalongkorn University, Bangkok, Thailand Kitipat.s@chula.ac.th Heat Exchanger Network (HEN) optimization by mathematical programming has been developed for half of century. Stage-wise superstructure by Yee and Grossmann (1990) is one of the most popular HEN design models dividing process temperature into stages under assumption of constant specific heat capacity (Cp). To make model more realistic, the model is developed under temperature-dependent Cp which is a cubic equation of temperature. To simplify this circumstance, non-linearity of Cp is linearized by applying partitions of linear equations of Cp at different temperature intervals. HEN from the modified stage-wise superstructure model with variable Cp by GAMS software version 24.2 is validated with crude preheat train case from commercial simulation program (Pro/II software version 9.1). The Cp from Pro/II software as a function of temperature is divided into three linear equations of Cp at three temperature intervals. A case study of crude preheat train without HEN consists of five hot product streams and one cold crude oil stream, consuming hot and cold utilities of 175,353.53 kW and 107,945.89 kW, respectively. The new model synthesizes HEN for the crude preheat train, consuming less hot and cold utilities of 82,100.93 kW and 14,654.38 kW. The synthesized HEN is validated by Pro/II software, giving duties of hot and cold utilities of 81,835.30 kW and 14,417.20 kW. Percent error of hot and cold utilities of HEN from GAMS model are only 0.31 % and 1.65 %, compared to ones from Pro/II software. Compared to Pro/II results, percent error of process HEN area and utility HEN area from new model are 1.22 % and 2.91 %, respectively. New model with variable Cp assumption gives a good agreement with validated HEN by Pro/II. Compared to examples from other publications, the model shows that it generates HEN design with less TAC and less computational time. 1. Introduction Previously, Heat Exchanger Network (HEN) design is based on heuristic technique and has been developed since late nineteenth century. Linnhoff and Flower (1978) proposed Problem Table Algorithm which is important for heat analysis and become one of major keys in Pinch Analysis. However, Pinch Design method is time- consuming. Mathematical programming for HEN is applied to overcome this limitation and generates good solution as well. There were many HEN design/retrofit models invented like one from Escobar and Trierweiler (2013) who studied mathematical programming models. One of the most interesting model is a stage-wise superstructure model (or SYNHEAT) invented by Yee and Grossmann (1990) which has been studied with fouling effects by many authors such as Yanyongsak and Siemanond (2018). This SYNHEAT model has temperature as variable along each temperature stage. The model contains only one non-linear equation from objective function of TAC, subject to robust linear constraint. The limitations of SYNHEAT model are from assumptions of isothermal mixing and constant heat capacity flow rate. Zhu and Asante (1999) and Ayotte- Sauvé et al. (2017) developed their model for HEN retrofit with linearization of enthalpy by keeping Cp constant in each temperature interval using high number of interval for more accuracy. Li et al. (2012) studied SYNHEAT model using non-linear Cp as cubic equation. They synthesized non-linear cubic equation into one linear equation combined with genetic algorithm and showed that it was possible to run SYNHEAT model with variable Cp. Kim and Bagajewicz (2017) used their new developed model with cubic equation of Cp for HEN synthesis. However, their models used large computational time because they use substage over another stage. This work presents novel technique for temperature-dependent Cp, called partitioning technique, is applied with stage- wise superstructure model to synthesize HEN. DOI: 10.3303/CET1976059 Paper Received: 16/03/2019; Revised: 02/07/2019; Accepted: 02/07/2019 Please cite this article as: Valeekiatkul S., Siemanond K., 2019, Mathematical Programming for Heat Exchanger Network Design with Temperature-Dependent Specific Heat Capacity, Chemical Engineering Transactions, 76, 349-354 DOI:10.3303/CET1976059 349 2. Problem statement For this study, crude preheat train with minimum annual cost is designed using the modified stage-wise superstructure model with linearization form of non-linear Cp of cold crude oil and five hot products of naphtha, kerosene, diesel, gasoil and residue; as shown in Figure 1. In this case, inlet/outlet temperature and flow rate of process streams are specified by product specification. Cp is collected from simulated crude distillation unit in commercial simulation program (Pro/II software). Hot/cold utilities are steam and cooling water. Exchanger minimum approach temperature (EMAT) is set at 10 °C. For the case assumption, stream splitting is allowed with maximum of three splits at each stage interval for cold stream but not allowed for hot streams. Finally, economic analysis for HEN design is performed using cost parameters from published paper to target the optimum solution. (a) (b) Figure 1: Cp of crude oil: (a) temperature-dependent Cp (b) linearization form of temperature-dependent Cp 3. Methodology To design HEN, stage-wise superstructure from Yee and Grossmann (1990) is modified under temperature- dependent Cp using mixed integer non-linear programming (MINLP) from general algebraic modeling system (GAMS). 3.1 Equations When Cp is varied with temperature, overall energy balance equation and stage energy balance of stage-wise superstructure must be redefined as average Cp and flow rate. Average Cp can be calculated carefully from stage temperature thi,k and thi,k+1 based on linear equation in each partitions. Weight average calculation is selected to average Cp of the stage for more accuracy. Piece-wise linearization on Cp by Zhu and Asante (1999) and Ayotte-Sauvé et al. (2017) used constant average Cp in each segments as shown in Figure 2b. However, new technique uses weight average calculation as shown in Figure 2a. It shows that new Cp linearization technique give higher accuracy for Cp calculation. (a) (b) Figure 2: Example calculation of Cp linearization technique: (a) new Cp linearization technique (b) Cp linearization technique from Zhu and Asante (1999) and Ayotte-Sauvé et al. (2017). 350 Process stream equations Eq(1) and Eq(2) calculate average Cp of stage by weighted average calculation from each partition. Specific heat capacity is normally calculated by average mean temperature and they are set to calculate by linear relationship from each partition. Variable Actin,i,k is applied to specify used partition. It is positive when temperature stage k and k+1 are located in used partition, else negative. 𝐶𝑃𝐻𝑜𝑡𝑖,𝑘 = ∑ [𝐶𝑃𝐻𝑛,𝑖,𝑘 × max(0,𝐴𝑐𝑡𝑖𝐻𝑛,𝑖,𝑘)] 𝑁𝑂𝑃 𝑛=1 /∑ [max(0,𝐴𝑐𝑡𝑖𝐻𝑛,𝑖,𝑘)] 𝑁𝑂𝑃 𝑛=1 (1) 𝐶𝑃𝐶𝑜𝑙𝑑𝑗,𝑘 = ∑ [𝐶𝑃𝐶𝑛,𝑗,𝑘 × max⁡(0,𝐴𝑐𝑡𝑖𝐶𝑛,𝑗,𝑘)] 𝑁𝑂𝑃 𝑛=1 /∑ [max(0,𝐴𝑐𝑡𝑖𝐶𝑛,𝑗,𝑘)] 𝑁𝑂𝑃 𝑛=1 (2) Utility streams equations Eq(3) and Eq(4) are separated from Eq(1) and Eq(2) because they need specific average Cp for last stage target temperature of cold utility and first stage target temperature of hot utility. 𝐶𝑃𝐻𝑜𝑡𝐿𝑖,𝑁𝑂𝐾+1 = ∑ [𝐶𝑃𝐻𝐿𝑛,𝑖,𝑁𝑂𝐾+1 × max⁡(0,𝐴𝑐𝑡𝑖𝐻𝐿𝑛,𝑖,𝑁𝑂𝐾+1] 𝑁𝑂𝑃 𝑛=1 /∑ [max(0,𝐴𝑐𝑡𝑖𝐻𝐿𝑛,𝑖,𝑁𝑂𝐾+1)] 𝑁𝑂𝑃 𝑛=1 (3) 𝐶𝑃𝐶𝑜𝑙𝑑𝐹𝑗,1 = ∑ [𝐶𝑃𝐶𝐹𝑛,𝑗,1 × max⁡(0,𝐴𝑐𝑡𝑖𝐶𝐹𝑛,𝑗,1)] 𝑁𝑂𝑃 𝑛=1 / ∑ [max(0,𝐴𝑐𝑡𝑖𝐶𝐹𝑛,𝑗,1)] 𝑁𝑂𝑃 𝑛=1 (4) 3.2 Concept of partitioning technique to approximate Cp Figure 3 shows example of Cp calculation for hot process stream i in stage k by partitioning technique. Suppose that a hot process stream has supply temperature (thi,k) and target temperature (thi,k+1) varied in the temperature range from partition 1 to partition 2, as shown in Figure 3. First, ActiHn,i,k of each process temperature (thi,k to thi,k+1) occupied partition becomes positive value and temperature difference of the partition is calculated. They show that enthalpy; ActiH1,i,k (FixtempH1 - thi,k+1), and enthalpy; ActiH2,i,k (thi,k - FixtempH1), are positive except enthalpy; ActiH3,i,k (FixtempH2 - thi,k+1), is negative value. Maximum operation force negative value to zero value (enthalpy of partition 3). Next, CPHn,i,k are determined by average mean temperature calculation between thi,k and thi,k+1. Finally, CPHoti,k is calculated from CPH1,i,k and CPH2,i,k by weighted average calculation concept as shown in Figure 3. The advantages of new technique are linear equation fitting very well with data and using small number of partitions (only 3 partition is enough for this case study). Moreover, average Cp calculation for each stage has high accuracy and it is free to calculate without any fixed Cp data from user. Figure 3: Example of partitioning technique for hot process stream 3.3 Case study First case study is crude distillation unit from Pro/II library selected to design crude preheat train. Figure 1a shows non-linearity of Cp of crude oil preheated from 50 to 376.8 °C. Temperature range of each process stream is linearized by partitioning to three linear sections (Cpn = AnTmean + Bn) shown in Figure 1b and Table 1. All three partitions of linear Cp are set to fit the curvature of Cp profile to reduce complexity of the model. This example aims to validate HEN design results from GAMS by Pro/II software using heat flow rate of each exchanger by fixed heat duty and HEN topology from GAMS results. Then, utility cost, capital cost and total annual cost (TAC) are compared between Pro/II and GAMS results. Grid diagrams of HEN design by GAMS and HEN validation by Pro/II are shown in Figure 4. Table 2 shows HEN designed by GAMS giving slightly small error in utility usage compared to validated one by Pro/II because of linearization with three partitioning technique in the GAMS model. 351 Table 1: Specific heat capacity linearization (Cpn = An×Tmean + Bn). Streams Partition Number (n) An Bn Fixed Temp 1 (°C) Fixed Temp 2 (°C) Cpaverage (kJ/kg °C) H1 1 0 2.01157288 - - 2.01157288 H2 1 0.00394085 1.74992049 145.23 76.27 2.19753245 2 0.00402349 1.74851827 3 0.00389833 1.75594557 H3 1 0.00348348 1.78714691 195.21 95.08 2.28077031 2 0.00396272 1.70138615 3 0.00371589 1.70777188 H4 1 0.00326872 1.81109354 243.24 114.86 2.42342119 2 0.00377460 1.69563864 3 0.00402867 1.66442583 H5 1 0.00301501 1.83347265 262.47 122.78 2.41611128 2 0.00359289 1.67211858 3 0.00408324 1.60983833 C1 1 0.00339671 2.13223103 168.18 112.12 2.76662173 2 0.00908254 1.23185053 3 0.00411179 1.69174821 Figure 4: Validated HEN by Pro/II for cast study 1 Table 2: Duty and area data comparison between GAMS and Pro/II of partition technique Heat Exchanger Duty (kW) GAMS Duty (kW) Pro/II Percent Error (%) Area (m2) GAMS Area (m2) Pro/II Percent Error (%) E1 10,086.30 10,086.30 - 271.91 273.10 0.44 E2 43,176.10 43,176.10 - 1,175.95 1,181.75 0.49 E3 5,561.60 5,561.60 - 139.37 140.01 0.46 E4 3,660.20 3,660.20 - 248.75 254.51 2.26 E5 6,325.80 6,325.80 - 287.14 288.38 0.43 E6 4,584.50 4,584.50 - 255.83 259.38 1.37 E7 20,124.00 20,124.00 - 1,023.36 1,047.05 2.26 CU1 1,378.39 1,378.40 0.00 131.50 131.49 0.01 CU2 2,089.76 2,074.50 0.74 130.63 129.45 0.91 CU3 2,149.99 2,131.50 0.87 130.36 129.16 0.93 CU4 4,868.91 4,831.00 0.78 130.45 129.46 0.76 CU5 4,167.33 4,001.80 4.14 168.01 162.68 3.28 HU1 82,100.93 81,846.20 0.31 904.99 902.03 0.33 352 Table 3: Total annual cost data comparison between GAMS and Pro/II results Cases GAMS TAC ($/y) Pro/II TAC ($/y) TAC Error (%) Capital Cost Error (%) Utility Cost Error (%) 3 Partitions 6,554,118 6,536,373 0.27 0.08 0.33 Table 4: Utility usage and area data comparison between GAMS and Pro/II results Cases Cold Utility (kW) Hot Utility (kW) Process Area (m2) Utility Area (m2) GAMS Pro/II Error GAMS Pro/II Error GAMS Pro/II Error GAMS Pro/II Error 3 Partitions 14,654 14,417 1.65 % 82,100 81,846 0.31 % 3,402 3,444 1.22 % 1,630 1,584 2.91 % On the other hand, area of each exchanger and Cp calculation from HEN by new GAMS model and validated one by Pro/II program are not much different (Table 2) because of model assumption of temperature-dependent Cp. Thus, TAC from GAMS and Pro/II program is 6,554,118 $/y and 6,536,373 $/y. TAC error is only 0.27 % illustrated in Table 3. Table 4 compares results of cold utility, hot utility, process area and utility area between GAMS and Pro/II program. It can be observed that the error is really small. Therefore, new model is always getting more accurate results in every result like temperature of each partition, area of exchanger and Cp calculation of each partition. All of this accuracy variables are needed in real HEN design that could be found only in temperature-dependent specific heat capacity model. Generally, Cp is an empirical correlation in form of cubic equation as a function of temperature calculated from fitting experimental data but it makes model more non-linear and difficult to be solved by MINLP solver. Another case study from Kim and Bagajewicz (2017) includes 11 hot streams, 2 cold streams and it is used to compared between our partitioning technique and ordinary non-linear cubic equation technique. EMAT is set at 10 °C. No splitting is allowed for hot streams and four splitting is allowed for cold streams. Figure 5: New HEN design result from Kim and Bagajewicz (2017) case study. 13 heat exchangers, 1 hot utility and 4 cold utilities exchanger are installed in this case study which show HEN in Figure 5. TAC of new solution is 3,229,804 $/y from capital cost 878,538 $/y and utility cost 2,351,265 $/y. Compared to previous solution, 2 more installed exchangers increase capital cost from 874,980 to 878,538 $/y 353 but the results show that utility cost decrease from 2,576,605 to 2,351,265 $/y. It reduces TAC around 6.42 % (3,451,585 to 3,229,804 $/y) as shown in Table 5. Table 5: Total annual cost data comparison between GAMS and Pro/II results Cases TAC ($/y) Capital Cost ($/y) Utility Cost ($/y) Previous Solution (Kim and Bagajewicz, 2017)) 3,451,585 874,980 2,576,605 New Solution 3,229,804 878,538 2,351,265 4. Conclusions Cp should be set as temperature-dependent variable. To reduce complexity of cubic equation and improving accuracy of Cp calculation by previous piece-wise technique for specific heat capacity, general solution from other publisher, linearization by dividing the specific heat capacity into partition is used here. It shows that new model reduces error in TAC (representative of overall error) to 0.27 % and temperature outlet of every heat exchanger is almost the same effect on higher accuracy of area calculation. Moreover, another example from other publications (Kim and Bagajewicz, 2017)) are used to calculate compare with this new model and it gives slightly lower TAC (improving 6.42 %) and computational time. HEN retrofit can be applied by this technique in the future work. Nomenclature Indices Parameters i Index for hot process stream A Slope coefficient of Cp j Index for cold process stream B Y-intercept of Cp k Index for stages (1, …, NOK) Fixtemp Temperature at edge of partition, °C n Number of partition (1, …, NOP) Variables ActiH Activated variable of hot stream (Positive when exchanger locate in the partition), °C ActiC Activated variable of cold stream (Positive when exchanger locate in the partition), °C ActiHL Activated variable of cold utility (Positive when exchanger locate in the partition), °C ActiCF Activated variable of hot utility (Positive when exchanger locate in the partition), °C CPH Average Cp calculation in the partition of hot stream based on average temperature, kJ/(kg°C) CPC Average Cp calculation in the partition of cold stream based on average temperature, kJ/(kg°C) CPHL Average Cp calculation in the partition of cold utility based on average temperature, kJ/(kg°C) CPCF Average Cp calculation in the partition of hot utility based on average temperature, kJ/(kg°C) CPHot Cp at stage of hot stream, kJ/(kg°C) CPHotL Cp at stage of cold utility, kJ/(kg°C) CPCold Cp at stage of cold stream, kJ/(kg°C) CPColdF Cp at stage of hot utility, kJ/(kg°C) Acknowledgments I am very grateful for Chulalongkorn university’s Rachadapisaek Sompote Fund (2017) for this research and I sincerely thankful to Mr. Natchanon Angsutorn for his guidance. References Ayotte-Sauvé E., Ashrafi O., Bédard S., Rohani N., 2017, Optimal retrofit of heat exchanger networks: A stepwise approach, Computers & Chemical Engineering, 106, 243-268. Escobar M. and Trierweiler J. O., 2013, Optimal heat exchanger network synthesis: A case study comparison, Applied Thermal Engineering, 51(1), 801-826. Kim S. Y. and Bagajewicz M., 2017, Global Optimization of Heat Exchanger Networks. Part 2: Stages/Substages Superstructure with Variable Cp, Industrial & Engineering Chemistry Research, 56(20), 5958-5969. Li G., Luo Y., Xia Y., Hua B., 2012, Improvement on the Simultaneous Optimization Approach for Heat Exchanger Network Synthesis, Industrial & Engineering Chemistry Research, 51(18), 6455-6460. Linnhoff B. and Flower J. R., 1978, Synthesis of heat exchanger networks: I. Systematic generation of energy optimal networks, AIChE Journal, 24(4), 633-642. Yanyongsak J. and Siemanond K., 2018, Heat Exchanger Network Retrofit Under Fouling Effects with Cleaning Schedule, Chemical Engineering Transactions, 70, 1549-1554. Yee T. F. and Grossmann I. E., 1990, Simultaneous optimization models for heat integration—II. Heat exchanger network synthesis, Computers & Chemical Engineering, 14(10), 1165-1184. 354