A computational tool for guiding retrofit projects of industrial heat recovery systems subject to variation in operating conditions

Heat exchanger networks (HEN) in industrial heat recovery systems often consist of large and complex subsystems. Usually, such HENs are subject to variation in operating conditions, such as varying inlet conditions or changing heat capacity flow rates. Additionally, complexities such as stream splits and recycle loops are commonly present in industrial HENs. Therefore, extensive modelling and/or analytical calculations may be necessary when analyzing different retrofit proposals. Furthermore, retrofit opportunities in industrial heat recovery systems are often constrained by operability considerations, i.e. retrofit actions are supposed to have as little impact as possible on the production process to maintain the quality of the core product. In this work, a computational analysis tool is proposed for effective screening of HEN retrofit design proposals at an early stage in the design process. The proposed tool enables fast evaluation of the network’s response, i.e. temperatures and heat loads, when operating conditions change and/or operational settings are manipulated, and it is applicable for a wide range of HEN structures. The practical use of the analysis tool is demonstrated in a case study on the HENs of a large modern Kraft pulp mill.


Introduction
In accordance with the 2018 amendment of the EU's Directive on Energy Efficiency, the European Commission demands their member states to increase the efficiency of final energy use by 32.5% by 2030, compared to 2007 [1]. In Sweden, the national energy agency has been commissioned by the government to ensure a reduction in energy intensity (primary energy consumption per unit of GDP) by 50% by 2030 compared to 2005 [2]. In 2017, the total Swedish final energy use was reported to be 378 TWh of which 143 TWh in industry [2]. The process industry was responsible for a significant share of the final industrial energy use in Sweden, and especially the pulp and paper industry accounted for more than 50% [2]. Moreover, a number of studies have shown that there is a substantial potential for energy savings in the pulp and paper industry [3][4][5][6][7][8][9][10][11][12][13][14][15]. Increased heat recovery is one important measure for increasing the energy efficiency in the pulp and paper industry. It has been estimated that about 25% reduction in steam demand can be achieved by heat integration measures in a typical chemical market pulp mill [11].
Industrial heat recovery systems often consist of large and complex heat exchanger networks (HEN) to transfer heat between different process streams. In pulp and paper mills, secondary heating systems are often implemented as a complement to direct process-to-process heat exchange. In these secondary heating systems, a heat transfer medium (commonly water) is used to recover/distribute heat from/to primary processes, in order to avoid operability issues resulting from the direct interconnection of different process streams. One example of such operability issues is that cooling and heating demands may vary over https://doi.org/10.1016/j.applthermaleng.2020.115648 Received 6 March 2020; Received in revised form 23 April 2020; Accepted 22 June 2020 time. The heat transfer medium can be stored in tanks to enable decoupling of the cooling and heating demands. The resulting hot and warm water systems are used for secondary heat recovery, but also for distributing process water for use in, e.g. dilution and washing steps. The interconnections between the secondary heating system and the primary processes are commonly characterized by a high degree of complexity, e.g. through stream splits, recycle and closed circulation loops.
Industrial plants are subjected to different kinds of time-dependent variations [5]. These variations can occur on very different timescales and are often characterized as short-, medium-or long-term variations [5]. Variations in the short term usually relate to disturbances in the process. These disturbances can have various reasons and one example is the daily variation of the ambient air temperature. In comparison to short-term variations, variations in the medium term relate to seasonal changes. In the case of long-term variations, an example is changes in legislation which can affect strategic business decisions and result in operational but also design changes of the mill. An example for how a mill's secondary heating system can be affected by different kinds of time-dependent variations is shown in Fig. 1. Fig. 1 shows measurements of both the ambient air temperature and the temperature of intake river water at a pulp mill in south-east Sweden for the year 2017. Both temperatures are inlet conditions for the mill's secondary heating system, i.e. intake water and air taken from the ambient are heated by means of the secondary heating system. In both temperature curves, the seasonal influence can be observed, i.e. medium-term variations in accordance with the description above. Additionally, the high fluctuation of the ambient air temperature implies that the mill's secondary heating system is also exposed to short-term variations.
A typical strategic business decision which usually affects the secondary heating system of a pulp mill in the long term is an increase of the mill's production capacity. Retrofitting of the secondary heating system may be a necessary consequence. Furthermore, flows and temperatures in the secondary heating system may be affected when the mill's production capacity is increased. Legislative incentives for increased energy recovery may be another long-term development affecting the secondary heating system of a mill as revamping of the existing system may be required.
Industrial heat recovery systems are normally designed to cope with expected short-and medium-term variations to ensure operability. It is essential that retrofit measures have as little impact as possible on the operability of the heat recovery system and thus on the core production process to maintain the quality of the manufactured products.
For investigating heat integration opportunities in the pulp and paper industry, graphical approaches (e.g. approaches based on the graphical pinch methodology) are commonly used. The graphical pinch methodology for process integration originates from the 1970s with the works of Linnhoff and Flower [16,17] and has been developed further in numerous publications. A comprehensive overview on the development of the pinch methodology for process integration can be found in the Handbook of Process Integration [18] and in the recent review by Klemeš et al. [19].
Several examples of heat integration studies based on graphical approaches for different kind of pulp mills which are based on (annual) average and steady state values can be found in the literature (see e.g. [3,4,[6][7][8]11]). The results of these studies suggest a potential of up to 25% of reduction in steam use. However, such studies often assume year-round steady-state operating conditions. Persson and Berntsson [5] showed that considering annual average values to calculate the steam savings potential of a retrofit heat integration project in a pulp mill can lead to an overestimation of 15% compared to the steam savings potential calculated with monthly average values (accounting for seasonal variations).
Dealing with variation in operating conditions when designing HENs is a well-studied problem and an overview of available methodologies for the synthesis of flexible HENs is provided by Kang and Liu [20] in a recent review paper. Flexible HENs remain, by definition, feasible for a-priori defined variations of operating conditions. In this context, feasibility is achieved if a-priori defined target values, e.g. stream target temperatures, can be reached for the entire span of variations. For the synthesis of flexible HENs, mathematical programming is a powerful tool due to the possibility to efficiently incorporate multiple operating periods during the synthesis process. Floudas and Grossmann [21] first introduced this multiperiod approach to the HEN synthesis problem which was developed further in several publications in the late 1980s [22,23]. Further work on the multiperiod approach has been done [24][25][26] and it has been applied to industrial case studies [27].
Compared to the numerous publications concerning synthesis of flexible HENs for greenfield problems, very little is found relating to retrofitting HENs which are subject to variation in operating conditions. Papalexandri et al. [28] developed a multiperiod MI(N)LP model which is based on multiperiod hyperstructures in which all retrofit alternatives are accounted for. To obtain the retrofitted HEN, which is operable for a-priori defined variations and yielding minimum total annualized cost, an iterative scheme was developed between the multiperiod MI(N)LP and a flexibility analysis subproblem. Another retrofitting methodology was reported by Kang and Liu [22], who introduced a two-step method to achieve retrofit design proposals which can operate cost-efficiently over multiple periods. In a first step, the multiperiod HEN synthesis model is solved as a greenfield design problem. In the second step, existing exchangers are relocated in order to meet required area demands identified in step one. However, size and complexity of industrial heat recovery systems very often limit the application of methodologies based on mathematical programming to industrial case studies.
In the case of graphical approaches, a common strategy is to develop different optimal design proposals for a number of selected sets of operating conditions [5]. The different design proposals are then either compared or combined to achieve an operable and energy efficient design for all operating points considered. Another strategy is to develop retrofit proposals by applying a graphical retrofitting methodology (e.g. methodologies based on heat transfer enhancement [29], temperature driving force curves [30] or identifying retrofit bridges [31]) for a specific nominal point and analyze the network's response to variations in a separate analysis step. Recently, Lal et al. [32] applied this approach to obtain insights and to identify the best performing retrofit design proposal by means of Monte Carlo simulation. Additionally, Langner et al. [33] proposed to combine flexibility and energy analysis to evaluate graphically derived retrofit proposals with respect to flexibility and energy efficiency. However, to evaluate the different designs over various operating conditions, extensive simulation (compare Lal et al. [32]) and/or iterative calculations (compare Langner et al. [33]) is required to determine the response of the network, i.e. the network temperatures and heat loads of exchangers, when the different designs are exposed to variations and operational settings are varied.

Aim
The aim of this paper is to propose a computational tool based on automated mathematical modelling of HENs that allows for calculation of temperatures and heat loads in a HEN that is exposed to variations in inlet conditions, e.g. supply temperatures. By means of this, the proposed tool can also be utilized to assess the effect of adjustments of operational parameters of a HEN such as split ratios. The tool is intended as an alternative to detailed simulation models and/or iterative (trial-and-error) calculations which can be computationally burdensome, especially for large-scale problems. Furthermore, the tool includes a systematic method for defining complex network topologies. As a result, the complete set of equations necessary to describe the heat and mass balances of a given HEN can be generated automatically. This enables the analysis of HENs of any size, with stream splitting and mixing, and with recirculating streams and closed circulation loops. The proposed analysis tool can be utilized in different use cases when working with HENs, i.e. the analysis tool is of interest if simulation of the HEN response is demanded. One example for such a use case is the screening process of different HEN design alternatives at an early stage in the design process (see Section 2). In this paper, usage of the proposed analysis tool is demonstrated in the screening process of a retrofit design project on the HENs of the secondary heating system of a Kraft pulp mill which is presented in a case study in Section 3 of this paper.

The computational analysis tool
The proposed analysis tool determines temperatures and heat exchanger (HEX) duties including the duties of utility exchangers (i.e. duties of heaters and coolers) of a HEN given the inlet conditions (e.g. supply temperatures), network topology, HEX areas and the operational setting of network parameters, such as split ratios or HEX bypass ratios. As mentioned in Section 1.1, it can be used to screen HEN retrofit measures as well as HEN greenfield design proposals. Such screening processes usually aim at providing information to the designer regarding the performance and/or costs of different design alternatives. Consequently, different design alternatives need to be collected prior to a screening process. The design alternatives can be derived using different methodologies, which is interesting for industrial HENs subject to varying operating conditions. As outlined in Section 1, industrial HENs are commonly characterized by complex network structures which leads to difficulties when applying design methodologies accounting for variations in operating conditions especially for retrofitting problems. On the other hand, graphical retrofitting methodologies which are based on average values are advantageous for structurally complex problems due to their beneficial user interaction. Consequently, instead of accounting for the variations in operating data during the definition of the design alternatives, the computational tool can be utilized (in a subsequent step) to analyze the network's response of each (graphically derived) design alternative when operating conditions vary. For retrofitting an industrial HEN that is subject to varying operating conditions, the following procedure is suggested which is further visualized in Fig. 2 energy efficiency) throughout the entire operating period by means of sensitivity analyses (see Section 2.3); and 5. Identification of the design alternative(s) which fulfills a defined objective (e.g. most cost-efficient and/or most energy-efficient).
In Fig. 2, the final step of the procedure is not illustrated. It should be noted that the proposed analysis tool does not generate design alternatives but rather enables the designer to identify bottlenecks in existing design proposals. Furthermore, the tool enables the designer to conveniently assess the effect of manipulating operational parameters such as split ratios and HEX bypass ratios. The main objective of such manipulation is to ensure operability (e.g. meet a-priori defined operational targets) and to optimize savings (e.g. steam savings, cost savings, etc.) for all (tested) operational data sets. By adopting a systematic methodology for defining any given network design topology, followed by automated network modelling, highly efficient screening of different design alternatives can be achieved at an early stage in the design process. Consequently, the proposed analysis tool allows for fair comparison of more design proposals under various operating conditions than if, e.g. time-consuming simulative approaches are chosen, which increases the time-efficiency of screening processes.

Theoretical background for automated HEN modelling
The purpose of the proposed analysis tool is to calculate all network temperatures and HEX duties (process-to-process and utility exchangers) for a given HEN. The number of unknown network temperatures can be specified a-priori [34]. For a network with N streams, n E process-to-process HEXs, n S stream splits, n M stream mixing points and n U utility exchangers, the number of unknown temperatures N T is: 2· .
The number of unknown temperatures N T is reduced by N if the supply temperatures for all process streams are known.
If the heat capacity flow rate CP i is constant and known for each stream as well as the UA-value (overall heat transfer coefficient multiplied by the heat surface area), type and flow arrangement of each process-to-process exchanger, 2 ⋅ n E linear equations in T (vector of the unknown network temperatures) can be obtained by the P-NTU method. The P-NTU method is used to model and solve rating equations for different types of HEXs [35]. A complete description of the set of equations can be found in Appendix A. Summarizing, the P-NTU method was used to model process-to-process HEXs.
In comparison to process-to-process HEXs, utility exchangers were modelled differently. Essentially, it was assumed that utility exchangers are not limited by heat transfer area. This assumption is in good accordance with secondary heating systems of pulp mills where e.g. steam can be injected to heat process (water) streams. In such a case no measurable heat transfer area is present. Therefore, the modelling of utility exchangers was done based on energy balances on the process side of the utility exchangers. If the heat load Q u or outlet temperature T out,u,i of the respective process stream for each utility exchanger is given, n U additional equations linear in T are obtained via energy balances on the process side of the utility exchangers. The energy balance for a given utility exchanger is: In Equation (2), CP i is the CP-value of process stream i and T in,u is the temperature of the process stream at the inlet of the utility exchanger. The heat load Q u in Equation (2) will be positive if the respective process stream is cooled. If the process stream is heated, the sign of the temperatures needs to be switched in order to obtain a positive value for Q u .
For the two unknown temperatures after a stream split (T out,1 and T out,2 ), trivial equations can be derived since isothermal splitting can be assumed. Thus, 2 ⋅ n S additional linear equations in T are obtained. The trivial equations can be expressed as follows: Mixing of two streams gives one additional energy balance equation at the mixing point which is linear in T if the CP-value after mixing CP mix,m and the CP-values of the mixed streams (CP in,1 and CP in,2 ) are known. If the energy balance is solved for the outlet temperature after the mixing point T out,m , it is expressed in the following form: In Equation (4), T mix,1 and T mix,2 are the temperatures of the two streams which are mixed. Overall, the number of equations is equal to the number of unknown temperatures N T as defined in Equation (1).
Recirculation of streams or closed circulation loops result in an identity change of the respective stream. In this context, the identity of a stream is defined by whether the stream releases heat (hot stream) or receives heat (cold stream) in the HEXs on the stream. Thus, recirculation of a stream or a closed circulation loop results in a stream that changes from cold to hot or from hot to cold, respectively. If a stream's identity changes, the temperature T out,c and the CP-value CP out,c after the identity change are theoretically unknown. However, trivial pairs of equations describing the relation between the temperatures (T out,c and T in,c ) as well as the CP-values (CP out,c and CP in,c ) over the identity change c ∈ IC can be derived.
The relation between the temperatures and the CP-values over a stream's identity change is visualized in Fig. 3. A switch model known from flowsheet modelling is used. A well-known example is the switch model commonly used in flowsheet modelling software, e.g. Aspen HYSYS, to model recycle flows. Summarizing, recirculation of streams or closed circulation loops add the same number of unknown temperatures as of linear equations in temperature to the problem.
Under certain assumptions (known and constant CP-values and UAvalues), it is therefore possible to express the unknown network temperatures by means of linear equations which can be solved as a linear equation system.
As the overall heat transfer coefficients (U-values) of the HEXs usually vary with temperatures and/or CP-values, assuming a constant UA-value when variations in temperatures and flow rates are introduced introduces a simplification of the actual problem. To avoid this simplification, thermodynamic modelling and iterative processes may be necessary. However, it is assumed that in early design stage screening, a certain degree of simplification is a valid trade-off between accuracy and computational expense. Therefore, temperature-dependencies of U-values are neglected. On the other side, there exist simple approaches (e.g. correction factors) to account for the dependency of the U-value if given CP-values differ from design values. In the proposed tool, a correction factor which was introduced and studied by Persson and Berntsson in [5] has been implemented. This correction factor uses a simplified correlation between the heat transfer coefficient and and the flow rate of a stream. Further details can be found in [5].

Network modelling
To model a HEN, the network structure or topology needs to be described. This can be challenging for large and complex networks due to the presence of a high number of units, multiple stream splits, recirculation of streams and closed circulation loops. A general methodology has been developed to model the topology in a straightforward way. In a first step, each stream of the network is numbered. To handle stream splits, recirculation of streams and closed circulation loops, streams are numbered with respect to their CP-values and identities. This means that a stream in the modelling sense needs to have a constant CP-value and fixed identity. A stream split where one stream is split into two streams is, thus, counted as three individual streams. If the two streams resulting from the split are mixed back together a fourth stream is counted. Fig. 4 shows the classical grid representation of an example network. The streams are numbered according to the proposed numbering procedure. Furthermore, all unknown network temperatures with respect to Equation (1) are marked.
For recirculation of streams and closed circulation loops, switches are used which define the locations where the identity of streams changes. Consequently, with each switch a new stream is introduced. The general counting procedure is shown in Fig. 3, while a specific example can be seen in Fig. 4.
Besides stream numbering, a generalized way to refer to the location of a specific process-to-process HEX (P2P) or utility exchanger (utility) is needed. Basically, this includes the (process) streams which are exchanging heat by means of the exchanger as well as the location of the specific exchanger on those streams. For this reason, a system to specify the position of the exchangers is developed which is based on the stream numbering. For each exchanger, the hot and cold stream as well as the position of the exchanger on these streams must be known. For utility exchangers, target temperatures are specified in order to calculate the duty of each utility exchanger (see Section 2.1). For the network depicted in Fig. 4, the following location matrix for the HEXs can be derived (see Table 1).
In Table 1 and Fig. 4, the position of the HEXs is determined by looking from left to right which in this case is the grid flow direction of the hot streams. In general, different ways are possible and can be implemented in the tool. A similar positioning system is used to describe the location of stream splits and mixing of streams. It is also based on the introduced stream numbering system. Table 2 presents the location matrix for the stream splits and mixes for the example network in Fig. 4.
To describe the location of switches, the ingoing and the outgoing streams must be specified, and a similar location matrix can be derived. Table 3 shows the location matrix for switches of the example network in Fig. 4.
Based on the exchanger location matrix, the split/mix location matrix and the switch location matrix, the vector of the unknown network temperatures, T, can be sorted in a data structure called temperature matrix. In this data structure, the unknown network temperatures are allocated to the process streams. The allocation process is automated in the proposed analysis tool. The basis for this automatization is Equation (1). In this context, one element of T (unknown network temperatures, see Fig. 4) is allocated to a stream if the stream is connected to a HEX (P2P or utility). Furthermore, one element of T is allocated to each stream resulting from a split, mix or switch. This way, all unknown network temperatures can be allocated to the different streams. Therefore, one dimension in the temperature matrix represents the streams present in the HEN.
The temperature matrix for the example network is shown in Table 4. As mentioned previously, the rows of the temperature matrix represent the different process streams (row number equivalent to stream number in Fig. 4), while the number of columns represent the number of unknown network temperatures of the corresponding process stream. As the number of unknown network temperatures may be different for the individual process streams, a data structure which accounts for this must be used (e.g. "Cell Array" in MATLAB or "List of Lists" in PYTHON).
The temperature matrix is the basis for the automated HEN modelling strategy as it allows for allocating the elements of T (unknown network temperatures, see Fig. 4) to specific HEXs (P2P and utility), stream splits, mixing points and switches. The allocation of the unknown network temperatures is vital to be able to derive the set of energy and mass balances presented in Section 2.1 automatically with the correct elements of T. This is illustrated by means of HEX 1 of the example network shown in Fig. 4: Based on the exchanger location matrix (Table 1), the two process streams which are connected by HEX 1 are hot stream 5 and cold stream 11. Additionally, the position of HEX 1 on these two streams is specified in the exchanger location matrix. For hot stream 5, HEX 1 is the first exchanger (counting from left to right) on this stream. Thus, the hot stream inlet temperature of HEX 1 is the inlet temperature of stream 5, T in,5 , while the hot stream outlet temperature of HEX 1 is T7 (elements 1 and 2 for stream 5 in the temperature matrix, Table 4). For cold stream 11, HEX 1 is the second exchanger on this stream. Thus, the cold stream outlet temperature of HEX 1 is T22 while the cold stream inlet temperature of HEX 1 is T21 (elements 2 and 3 for stream 11 in the temperature matrix, Table 4). Consequently, the unique elements of T which are present in Equations (A.1) and (A.2) for HEX 1 can be identified.
As the example for HEX 1 shows, by means of the exchanger location matrix and the temperature matrix, specific elements of T are allocated to each exchanger. Similar allocations can be achieved by means of the split/mix or the switch location matrix to allocate unknown network temperatures to the corresponding split, mix or switch. Consequently, the set of equations to model an arbitrary HEN as a linear equation system (see Section 2.1) can be derived automatically.

Sensitivity analysis via sensitivity tables
The proposed analysis tool can be utilized to generate sensitivity tables as a well-arranged tool to analyze the influence of varying operating data on the network's response. These tables are generated by repeatedly calculating the network's response (i.e. temperatures and duties of HEXs) for defined variations in operating data but also design parameters (e.g. area of a new HEX) or operational settings (e.g. split ratio). The concept of HEN analysis based on sensitivity tables was introduced by Kotjabasakis and Linnhoff in [36]. In comparison to the manual derivation of the sensitivity tables introduced in [36], sensitivity tables can be generated automatically by means of the proposed analysis tool. An example of how sensitivity analysis based on  T2   T3  T4   T7   T6   T8  T9   T5   T14  T15   T16   T17  T18   T19   T20   T21   T22   T23   6   T12  8  T10  T11   10  sensitivity tables can support the designer during a retrofit screening process is given in Section 3.5.

Input data and limitations of the computational tool
All necessary input data required to analyze a given HEN by means of the proposed analysis tool is summarized in Table 5. Providing the data listed in Table 5, a mathematical model for a HEN can be derived automatically which is linear in the network temperatures and can be solved as a linear equation system (see Section 2.1).
Based on this input data, the proposed analysis tool enables the design engineer to calculate all network temperatures as well as duties of process-to-process and utility HEXs. As the above listed input data • Phase-changing streams are not considered; • The tool is limited to calculating the network's response to variations in network supply temperatures, heat capacity flow rates and UA-values as well as operational parameters such as HEX bypass ratios, split ratios and flow rates of utility streams (loads of utility exchangers); and • The tool can be used for conducting sensitivity analyses on all input parameters (see Section 3.5); optimization algorithms to identify the optimal setting of operational parameters to a certain objective, e.g. minimize steam demand, are not incorporated.

Case study
This section illustrates usage of the proposed automated HEN modelling strategy and the computational analysis tool for screening different retrofit design proposals for utilizing more excess heat in the secondary heating system of one of the largest Kraft pulp mills in Sweden. The respective mill has an annual pulp production capacity of 750,000 air dried tons. As raw material softwood and hardwood are used while around 75% of the annual pulp production is based on softwood. Besides pulp, other bio-based products such as turpentine and tall oil are produced and sold. Additionally, several back-pressure turbines and a condensing turbine are utilized to generate 800 GWh/a of electricity of which around one forth is exported to the grid. The secondary heating system of the pulp mill is based on water as a heat transfer medium and three main temperature levels are present (cold water: CW; warm water: WW; hot water: HW).
Prior to this study, a Pinch analysis was carried out based on average values for spring operating conditions and the heating savings targets for the pulp mill of interest were quantified (27.7 MW) as well as the magnitude of major Pinch rule violations in the existing network. The heating savings targets correspond roughly to 10% of the mill's current total heating demand. Spring conditions were chosen to avoid over-or underestimation of the potential savings which may occur when using more extreme summer or winter conditions. The two largest Pinch rule violations identified involved steam heating below the Pinch for primary boiler air heating (6 MW) and boiler feedwater heating (4.1 MW), representing 36.5% of the total identified Pinch rule violations. The two corresponding process streams (primary air to boiler and boiler feedwater) are located in the same independent subsystem within the secondary heating system of the pulp mill, which is shown in Fig. 5. In this subsystem, a third process stream (internal heating network) is also present.

Table 1
Exchanger location matrix of the example HEN for defining the position of all process-to-process (P2P) and utility exchangers and the target temperatures of streams with utility exchangers.    Table 4 Temperature matrix to allocate unknown network temperatures of example network to specific exchangers (process-to-process and utility), stream splits, mixing points and switches.  2  T1  T2  3  T3  T4  4  T5  T6  5  T in,5  T7  T8  T9  6  T10  T11  T12  7  T13  T in,7  8  T15  T14  9  T18  T17  T16  10  T19  11  T23  T22  T21  T20  T in, The three process streams are connected through the secondary heating system which is fed by CW from the corresponding tank and by a small amount of blowdown steam which is released from a boiler. The CW is heated by hot mist (small droplets of water suspended in air) from the recovery boiler system in HEX 2 yielding HW. The HW is then distributed to heat the primary boiler air through a closed circulation loop (HEX 1 and 3) and over a series of stream splits to the boiler feedwater heating (HEX 6 and 7) and the internal heating network (HEX 4 and 5). Additionally, a small amount of HW is extracted from the secondary heating system to the HW tank (stream split 2) from where it is later distributed to other processes, e.g. to the washing and dilution steps in the bleaching section of the pulp mill. The heat transfer medium in the closed circulation loop is also secondary heating water.
The closed circulation loop is a safety measurement to avoid a high mass flow of secondary heating water to the boiler in case of a leakage in HEX 1. After being heat exchanged with the three process streams in the respective HEXs, the HW has been cooled to the temperature level of WW and is recirculated to HEX 2.
In the studied subsystem, five locations were identified where the system is connected to other (sub)systems of the mill. For the analysis, the conditions of the streams entering the studied system at these locations were defined as inlet conditions. The locations are marked with circled numbers in Fig. 5.
To reduce the Pinch rule violations in the subsystem of interest (see Fig. 5), a retrofitting project was considered. Different retrofit design alternatives were suggested (based on graphical Pinch considerations) Table 5 Necessary input data for applying the computational analysis tool to a HEN.

For each heat exchanger
• UA-value (for process-to-process heat exchanger) • Type of heat exchanger (process-to-process or utility) • Hot and cold streams connected by the exchanger • Position number of heat exchanger on hot and cold stream • Individual heat transfer coefficients and CP values for the design case to calculate corrected UA-values by means of the chosen correction factor (for process-to-process heat  and the computational analysis tool was applied to screen the different design alternatives. Generally, the procedure presented in Section 2 was followed. However, as mentioned in Sections 2.1 and 2.4, the computational tool and the automated HEN modelling strategy rely on assumptions in order to achieve a linear equation system. These assumptions imply that the design specifications of the HEXs in the system (UA-values and individual film heat transfer coefficients) are known. Thus, design specifications of the HEXs in the subsystem of interest (see Fig. 5) were collected and the computational tool was used to calculate the temperatures (and flow rates resulting of stream splits) which were then compared to available measurement data. This validation step was done prior to the steps of the procedure described in Section 2 and is outlined in Section 3.1.

Validation of the HEX design data
The HEX design specifications which were derived in collaboration with mill experts are presented in Table 6. The accuracy of this data was validated using measurement data for spring conditions. Additionally, HEX 1 in Fig. 5 aggregates a number of different exchangers in the real plant layout. Therefore, the UA-value of HEX 1 was fitted according measurement data during the validation process. The automatically generated HEN model was solved using the computational tool and the calculated values were compared to the measurement data. Fig. 6 shows the measurement points in the studied subsystem. Table 7 presents the measured values and the calculated values.
In general, there is a good agreement between the measured and calculated data. This is indication that the available input data (including HEX design specifications) and the made assumptions were valid.

Retrofit measures to increase energy efficiency
In the next step, different retrofit proposals were defined based on graphical Pinch analysis considerations. The main objective was to reduce the Pinch rule violations identified in the Pinch analysis carried out prior to this study. It should be noted that investment costs were only considered qualitatively, i.e. it was assumed that new HEXs are more costly than repiping of existing HEXs.
One option is to implement two new HEXs (New1 and New2) which recover excess heat from the digester section of the pulp mill. Pressurized secondary heating water is heat exchanged in the digester section of the mill to generate very hot secondary heating water (VHW) which is then heat exchanged with the process streams in the studied subsystem. Two locations in the studied subsystem were identified for HEXs New1 and New2. With respect to the flow arrangement of the pressurized VHW, this resulted in two different design alternatives (designs 1 and 2), as shown in Fig. 7. By means of the proposed retrofit design alternatives, a sixth location is introduced at which operational conditions can be interpreted as inlet conditions for the studied subsystem. This location is marked accordingly.
Additionally, it was considered to reduce the investment cost by investing in only one new HEX (New1) but use the pressurized VHW in the existing HEX 6 and HEX 7 instead of the HW generated by HEX 2. With respect to the flow arrangement of HEX New1 and HEX 6/7, two additional design alternatives (designs 3 and 4) were identified which are shown in Fig. 8. In these two design alternatives, the HW, which was originally split in stream split 3 (see Fig. 5) to feed the hot sides of HEX 6 and HEX 7, bypasses the feedwater heating (HEX 6 and HEX 7) completely and flows directly to stream split 4. This is not shown in Fig. 8.

Identification of operational data sets
In collaboration with mill experts, consistent time periods for the operating year 2017/2018 were identified which represent typical conditions (e.g. ambient temperature) for the four annual seasons and during which softwood pulp was produced at full plant capacity. After identification of the operational periods, the seasonal variations of the inlet conditions were investigated (marked by circled numbers in Fig. 5, Fig. 7 and Fig. 8). In Table 8, average values of the inlet conditions for the identified operational periods are shown. Additionally, the average values of the sum of the steam heater duties (steam heaters are shown in Fig. 5) for each operational period is shown in Table 8. It should be noted that the total steam demand is presented as a heat flow with the unit MW since it represents the average value of the steam heater duties for each operational period.
It should be noted that the operating conditions considered for spring were similar to the operating conditions used for the Pinch analysis. In practice, spring and autumn conditions were sufficiently similar to consider the same data set for both seasons.

Generation of mathematical model
The necessary input data for applying the computational analysis tool is listed in Table 5. Flow rates and supply temperatures of streams (inlet conditions) are listed in Table 8 and the network's topology is also known (see Fig. 5 as well as Fig. 7 and Fig. 8 for retrofit design alternatives). Heat capacity and density of the respective flow rates were assumed to be constant. Furthermore, the design specifications of the existing HEXs are listed in Table 6. In accordance with the procedure proposed in Section 2, operational constraints were defined which constrain the operability of the system, i.e. feasible operation for each set of operating data. Furthermore, adjustable operational parameters were defined which can be controlled in order to achieve feasible operation. This was done in collaboration with mill process experts and further information regarding this can be found in Section 3.4.1.
In accordance with Table 5, design specifications for the new HEXs (New1 and New2) were also needed to be able to fully compare the different design alternatives. The analysis tool was used to perform a sensitivity analysis on the UA-values to determine their influence on the steam demand of the studied subsystem as well as on the operational constraints defined in Section 3.4.1 (see Section 3.4.2). Again, investment costs were not considered thoroughly in this analysis. Instead, the major focus of the analysis was on maximizing the steam savings. This was done to avoid limiting the analysis results to possibly not representative cost values. On the other hand, as there obviously exists a trade-off between steam savings and investment cost, it was decided to limit the total increase of the new HEXs by a preliminary defined value of 500 kW/K, i.e. the maximum sum of the UA-values of the HEXs New1 and New2 is fixed to 500 kW/K (see Section 3.4.2). This way the focus of the analysis remained on energy efficiency, while unrealistic high investment costs were avoided. The derived design specifications were used to define a base case design which was then subject to further analysis, i.e. it was exposed to all sets of operating data and the influence of the operational parameters on the steam demand was analyzed (see Section 3.5).

Definition of operational constraints and identification of operational parameters
Seven operational constraints were defined in collaboration with mill process experts in order to guarantee the operability of the studied subsystem and to some extent also of other processes such as the washing and dilution steps which use HW produced in the subsystem. The operational constraints are listed below, and the corresponding numerical values are presented in Table 9: 1. The heat transferred in HEX 2 is limited by the available heat and  the HEX size -This can be expressed by three constraints: a) Maximum duty of HEX 2; b) Specified outlet temperature of the HW from HEX 2; c) Maximum flow of secondary heating water through HEX 2; 2. The return temperature of the pressurized VHW has a lower limit to satisfy the energy balance around the HW tank; 3. Maximum flow rate on hot side of HEX 3; 4. Maximum flow rates on the hot sides of HEX 6 (a) and HEX 7 (b); and 5. The flow rate of HW to the HW tank is specified (stream split 2).
In order to ensure that the identified operational constraints are not violated when the subsystem is exposed to variations, the following adjustable operational parameters were identified: • Split ratios (SR) 1-5 (6 for designs 3 and 4); • Flow rate of heat transfer medium (secondary heating water) in closed circulation loop (hereafter referred to as flow_loop); and • Bypass ratio of HEXs New1 or New2 (only one of the HEXs needs to be bypassed to control the return temperature of the VHW to the HW tank).

Base case design
The analysis tool was used to perform a sensitivity analysis on the UA-values of the two new HEXs New1 and New2 to identify the influence of these UA-values on the steam demand of the studied subsystem as well as on the operational constraints defined in Section 3.4.1. This sensitivity analysis was done for one set of operating data (spring/autumn operating conditions). In order to conduct the sensitivity analysis, the identified adjustable operational parameters were fixed. Therefore, additional constraints were defined for the splits ratios (SR) corresponding to stream splits 1, 2 and 4: • SR 1 was set so that flow of HW through HEX 3 is maximized (see constraint 3); • SR 2 was set so that the correct amount of HW is sent to the HW tank (see constraint 5); and • SR 4 was set to avoid any steam demand in the internal heating network (if possible: SR 4 < = 1).
Additionally, the flow of secondary heating water through HEX 2 was maximized respecting constraint 1c. For SR 3, SR 5 and flow_loop, the settings were taken from the historical data for spring/autumn operating conditions (in accordance with constraints 4 a and b). For design alternatives 3 and 4, the split ratio SR 6 was set in accordance with constraint 4 to maximize the flow through HEX 7 due to the larger UA-value of HEX 7 compared to HEX 6 (compare Table 6).
For design alternatives 1 and 2, the sensitivity analysis was conducted individually for HEX New1 and HEX New2, i.e. while the UAvalue of HEX New1 was varied, the UA-value of HEX New2 was kept at 0 and vice versa. Therefore, the results of the sensitivity analyses are similar for design alternatives 1 and 2. For design alternatives 1 and 2, Fig. 7. First retrofit proposal and the resulting two design alternatives (designs 1 and 2) of the studied subsystem of the pulp mill's secondary heating system. the results of the sensitivity analysis are presented in Fig. 9. Since most of the constraint values are fixed, e.g. flow through HEX 2 (see above), the only values which need to be monitored during the sensitivity analysis are the duty of HEX 2 (constraint 1a) to reach the specified outlet temperature of HEX 2 (constraint 1b) and the return temperature of the VHW (constraint 2). For all investigated UA-values, the minimum return temperature of the VHW to the HW tank (95.8 °C for spring/autumn conditions, see Table 9) was not reached (101.8 °C for maximum increase of HEX New1 and 97.8 °C for maximum increase of HEX New2). Additionally, the maximum duty of HEX 2 (15 MW for spring/autumn conditions, see Table 9) was not reached. For both design alternatives, Fig. 9 illustrates that steam savings increase with increasing UA-values of HEX New1 or HEX New2. As mentioned in Section 3.4, it was decided to consider a maximum total increase of 500 kW/K for the HEXs New1 and New2.
In a next step, sensitivity analyses were conducted to identify the optimal distribution of the total increase of 500 kW/K between HEX New1 and HEX New2 for design alternatives 1 and 2. The results of these sensitivity analyses are shown in Fig. 10 for design alternative 1 and in Fig. 11 for design alternative 2. Fig. 10 and Fig. 11 show the decrease in steam demand as well as the return temperature of the VHW to the HW tank for different combinations of UA-values for the HEXs New1 and HEX New2. The minimum return temperature of the VHW for spring/autumn operating conditions is marked as a grey line.
For design alternative 1, the highest steam savings were achieved for a maximum UA-value of HEX New2 (UA-value of HEX New1 is 0). For design alternative 2, the highest steam savings were achieved for distributed UA among the HEXs New1 and New2 (UA-value HEX New1: 100 kW/K, UA-value HEX New2: 400 kW/K). The minimum return temperature of 95.8 °C is not reached in any of the investigated UAvalue options.
For design alternatives 3 and 4, the situation is different since only one new HEX (New1) is implemented, while the two existing HEXs 6 and 7 are repiped. In this context, a separate sensitivity analysis was conducted for each of these two design alternatives. The results are presented in Fig. 12 and Fig. 13. Again, besides the decrease in steam demand, the corresponding VHW return temperature and its feasible minimum for spring/autumn operating conditions (grey line) are shown.
For design alternatives 3 and 4, the maximum feasible UA-value of HEX New1 (i.e. without violating operational constraints) is lower than 500 kW/K. In case of design alternative 3, the maximum UA-value of HEX New1 is smaller than 175 kW/K (see Fig. 12). In case of design alternative 4, the maximum UA-value of HEX New1 is smaller than 200 kW/K (see Fig. 13). Choosing larger UA-values would result in a too low return temperature of the VHW to the HW tank. Summarizing, the UA-values of HEX New1 and HEX New2 yielding the largest reduction in steam demand for each design alternative (tested with spring/autumn operating conditions, fixed adjustable parameters and limited maximum UA-value increase) were chosen as design specifications for the base case designs. Table 10 presents the design specifications of the HEXs New1 and New2 and the corresponding reduction in steam demand compared to the current steam demand for spring/autumn operating conditions (see Table 8) for the four base case designs.

Comparison of the steam savings potential of the proposed retrofit design changes for the entire operating period
In the next step, the studied subsystem was exposed to each set of operating data identified (compare Table 8) and the system's response was evaluated. Additionally, changes in the adjustable parameters (compare Section 3.4.1) were considered. The main objective was to reduce the steam demand compared to the current operating conditions (compare Table 8) while satisfying all operational constraints. Again, sensitivity analyses were conducted to quantify the influence of the adjustable parameters on the operational constraints as well as the steam demand. The following adjustable parameters were subject to the analyses: • SR 6 for design alternatives 3 and 4, • SR 3 for design alternatives 1 and 2, • SR 5, • Flow_loop (maximum flow rate is 56.7 l/s, i.e. approx. 150% of current operating point). Section 3.5.1 presents detailed example sensitivity analyses on the three adjustable operational parameters of design alternative 1 for spring/autumn conditions. Additionally, the results of the sensitivity analyses on the adjustable parameters of the design alternatives 2-4 for Table 8 Location and description of the seasonal variations in the studied subsystem of the pulp mill's secondary heating system (including the two retrofit proposals) and seasonal steam demand of studied subsystem.  Table 9 Numerical values for the identified operational constraints. spring/autumn operating conditions are summarized in Table 11. In Section 3.5.2, the potential to reduce the steam demand of the studied subsystem by means of the design alternatives 1-4 is compared considering all identified operating periods.

Sensitivity analyses on the adjustable parameters for spring/autumn operating conditions
To analyze the influence of the adjustable parameters on the steam demand, sensitivity analyses based on sensitive tables were done. Table  B1, Table B2 and Table B3 in Appendix B present the sensitivity tables for the operational parameters SR 3, SR 5 and flow_loop for design alternative 1. The changes in steam demand shown in Table B1, Table B2 and Table B3 relate to the steam demand of the base case design of design alternative 1 for spring/autumn conditions (see Table 10). It is worth mentioning that the steam demand of the internal heating system is 0 kW in all base case designs for spring/autumn conditions. It could be identified that changing the split ratios SR 3 or SR 5 leads in an increased total steam demand (compare Table B1 and Table B2). This implies that further optimization of these two split ratios is either not possible (at all) or is not possible considering only the respective split ratio. However, Table B3 shows that the parameter flow_loop could be optimized with respect to the steam demand. Increasing the flow rate by 30% decreases the duty of the steam heater in the primary air heating system by around 100 kW. It can be noted that the largest decrease in duty is not achieved at the highest possible flow rate (see Table B3). Concluding, changes in the parameters SR 3 and 5 did not lead to any decrease in steam demand, while the steam demand could successfully be decreased by increasing the parameter flow_loop by 30%.
Similar sensitivity analyses (based on sensitivity tables) were conducted for design alternative 2 for spring/autumn conditions. For design alternatives 3 and 4, the situation for spring/autumn conditions     was different, since the potential of the VHW is almost entirely utilized (return temperature of VHW close to minimum value of 95.8 °C for the respective chosen design specifications of HEX New1). Any increase in the parameter flow_loop would therefore result in an infeasible low return temperature of the VHW. A trial and error analysis would be necessary to test whether the steam demand could be decreased when increasing the parameter flow_loop, while at the same time reducing the UA-value of HEX New1 to prevent the return temperature from dropping beneath the threshold value. However, this was not investigated. Table 11 presents the achieved reduction in steam demand and a summary of the adjusted operational settings for the four design alternatives considering spring/autumn operating conditions.

Comparison of steam savings of design alternatives 1-4 for the different operating cases
In a final step, the design alternatives were exposed to the other identified sets of operating conditions (summer and winter operating conditions, compare Table 8) and the adjustable parameters and their influence on the steam demand were analyzed (similar to the procedure presented in Section 3.5.1). Sensitivity analyses based on sensitivity tables of the adjustable parameters (SR3/SR6, SR5 and flow_loop) were carried out when there was potential for reducing the steam demand without violating the operational constraints. Fig. 14 shows the reduction in the total steam demand of the studied subsystem for the different operating cases by means of design alternatives 1-4.
From Fig. 14, it can be concluded that for spring/autumn conditions design alternative 2 achieves the largest reduction in steam demand. However, for summer and winter conditions, design alternative 4 achieves the largest reduction in steam demand. With design alternative 3, the smallest reduction in steam demand was achieved for each of the investigated operating period.
Comparing the achieved reduction in steam demand with the identified Pinch rule violations in the subsystem of interest (10.1 MW for spring/autumn operating conditions) reveals that these Pinch rule violations could be decreased by 47.6-54.2 % for spring/autumn operating conditions depending on the chosen design alternative.
It should be noted that the reduction in steam demand gives indication for the performance of the different retrofit design alternatives for the considered operating points. However, a fair comparison of the energy performance of the different retrofit design alternatives must include the holistic performance. Therefore, the seasonal and annual steam savings based on the reduction in steam demand were calculated.
For calculating the seasonal and annual steam savings potential, it was assumed that the annual production time of the pulp mill in question is 8000 h. However, as mentioned in Section 3.3 only the production of softwood pulp was considered while the mill in question produces softwood as well as hardwood pulp in production campaigns (only one production line is in place). Therefore, it was further assumed that the annual production time of softwood pulp is 6000 h (based on the annual production capacity of softwood pulp mentioned in the first paragraph of Section 3), and that the annual seasons are isochronous (equally long time intervals). Based on these assumptions and the results shown in Fig. 14, the seasonal and annual steam savings were calculated for each design alternative and are presented in Table 12.
Based on the steam savings shown in Table 12 and an analysis of qualitative indicators, design alternative 4 may be identified to be favorable since it demands the second smallest investment with respect to new HEX area (repiping and one new HEX unit with a UA-value of 175 kW/K) and achieves the highest steam savings for two of the four annual seasons (compare Table 12).

Discussion of the obtained results
It should be noted that detailed cost analyses are necessary to identify the cost-optimal design alternative considering investment and operating cost. However, as mentioned in the first paragraph of Section 3, the aim of this case study section was to demonstrate usage of the developed computational tool by means of an industrial application. One possibility to consider investment costs more thoroughly is to consider cost functions in the sensitivity analyses presented in Sections 3.4 and 3.5. Consequently, the definition of the base case designs (presented in Section 3.4) and the adjustment of the operational parameters (presented in Section 3.5) could be repeated in an iterative process until a sufficiently good solution (i.e. cost-efficient) is obtained. It should be noted that the proposed computational analysis tool is still likely to increase the time-efficiency of such an iterative procedure compared to more detailed simulation tools since the degree of detail is decreased to a reasonable level for screening processes. For the demonstration of the proposed tool, it was decided to avoid trial and error procedures and focus instead on reducing the identified Pinch rule violations for an a-prior limited investment. Nevertheless, based on the achieved results and qualitative reasoning regarding the difference in investment cost, one of the four design alternatives could be identified to be most favorable as presented in Section 3.5.2.

Conclusion
Industrial heat recovery systems are often based on complex HENs. Regardless of the network complexity, HENs are usually based on similar sets of equations which implies that large parts of the modelling process can be automated. In this paper, an automated network modelling strategy has been proposed which can be implemented in any high-level programming language. The modelling strategy can handle complexities commonly present in industrial HENs such as stream splits, closed loops or recirculation. It is based on a table-based representation of the HEN which can directly be applied to a process flowsheet and a transformation to the commonly used but limited griddiagram representation is not necessary.
The proposed HEN modelling strategy was implemented in a computational tool to guide retrofit projects of HENs subject to varying operating data. Depending on the complexity of the process, the design of a HEN exposed to variations can be a cumbersome problem to solve if state-of-the-art methods are applied (graphical or optimization-based approaches). For approaches based on graphical considerations, full simulation and/or analytical calculations are necessary to evaluate the network's response when operational data changes. This applies also for the evaluation of adjustable operational parameters and (minor) design changes. The analysis tool proposed in this paper overcomes these difficulties for early design stage screening processes by automating the calculations necessary for screening processes. Furthermore, a strategy has been proposed for applying the analysis tool in a screening process based on Pinch-based retrofit design methodologies to handle HEN retrofit problems subject to variation in operating conditions.
The proposed analysis tool was applied to a case study on the HENs of a modern Kraft pulp mill with the aim of increasing the energy efficiency. Since the mill in question operates all year round, sets of operating data were identified which represent typical operating conditions during the different annual seasons. For one set of operational data, four different retrofit design alternatives were developed using Pinch-based design considerations to increase the overall energy efficiency of the mill. The computational tool was then applied to effectively identify promising design specifications of the new HEXs for each design alternative. In a next step, the computational tool was used to evaluate the potential for reducing the steam demand of the mill by means of the different design alternatives when operating conditions change. To allow for a fair comparison between the different design alternatives, the computational tool was used to identify the influence of changes in adjustable operational parameters such as split ratios on the overall potential for decreasing the steam demand by means of sensitivity analyses based on sensitivity tables.

Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.
The P-NTU method is used to model and solve rating equations for different types of heat exchangers. An overview of the model and its fundamental equations is presented below and a more detailed description about the methodology and the equations can be found in [35].
The P-NTU method uses a thermal effectiveness which is individually defined for the hot (index h for hot) and cold (index c) sides of the exchanger: The denominator is the same for both definitions and represents the maximum possible temperature change for any of the fluids in the exchanger. Using the definitions above the heat transfer between the hot and cold fluids can be expressed as: Here, CP ( ) h and CP ( ) c are the heat capacity flow rates for the hot and cold fluid, respectively. When solving for specific heat exchangers the thermal effectiveness can be expressed in the form: R is defined as the ratio between the heat capacity flow rates according to:  [35].

Table B1
Sensitivity analysis on the split ratio of stream split 3 of the studied subsystem for design alternative 1; listed values correspond to steam savings or operational constraints. SR