**Heat Exchangers for Waste Heat Recovery**

Special Issue Editors **Petr Stehlik Zdenek Jegla**

MDPI • Basel • Beijing • Wuhan • Barcelona • Belgrade • Manchester • Tokyo • Cluj • Tianjin

*Special Issue Editors* Petr Stehlik Brno University of Technology Czech Republic

Zdenek Jegla Brno University of Technology Czech Republic

*Editorial Office* MDPI St. Alban-Anlage 66 4052 Basel, Switzerland

This is a reprint of articles from the Special Issue published online in the open access journal *Energies* (ISSN 1996-1073) (available at: https://www.mdpi.com/journal/energies/special issues/ heat exchangers for waste heat recovery).

For citation purposes, cite each article independently as indicated on the article page online and as indicated below:

LastName, A.A.; LastName, B.B.; LastName, C.C. Article Title. *Journal Name* **Year**, *Article Number*, Page Range.

**ISBN 978-3-03936-475-6 (Hbk) ISBN 978-3-03936-476-3 (PDF)**

c 2020 by the authors. Articles in this book are Open Access and distributed under the Creative Commons Attribution (CC BY) license, which allows users to download, copy and build upon published articles, as long as the author and publisher are properly credited, which ensures maximum dissemination and a wider impact of our publications.

The book as a whole is distributed by MDPI under the terms and conditions of the Creative Commons license CC BY-NC-ND.

## **Contents**


## **About the Special Issue Editors**

**Petr Stehlik** (Prof. Dr.) is a Professor and Director of the Institute of Process Engineering of the Faculty of Mechanical Engineering at the Brno University of Technology. He is the Vice President of the Czech Society of Chemical Engineers and member of AIChE and ASME. He is executive editor or guest editor of several international journals; member or chairmen of scientific committees of international conferences; reviewer of PhD theses at reputable foreign universities; coordinator or contractor of international and national research projects; and the author and co-author of numerous publications and several patents. His research, development, and application activities involve heat transfer and its applications (heat exchangers, combustion systems, waste heat utilization); thermal treatment and energy utilization of waste (waste-to-energy); and energy savings and environmental protection.

**Zdenek Jegla** (Assoc. Prof. Dr.) is currently an Associate Professor and Senior Researcher in the field of Process Engineering at the Institute of Process Engineering of the Faculty of Mechanical Engineering (FME) at the Brno University of Technology (BUT) and a Head of the Heat Transfer Research Group at the NETME Research Centre FME BUT. He has many years of research and industrial experience, and he is author or co-author of numerous papers and contributions in international journals and conferences. His research and development as well as application activities are aimed at applied and enhanced heat transfer; waste heat recovery systems and equipment; process fired heaters, boilers, and heat exchangers; simulation, optimization, and CFD applications in the process and power industry; and process and equipment design and integration for energy savings and emissions reduction.

## **Preface to "Heat Exchangers for Waste Heat Recovery"**

The economical and efficient recovery of waste heat produced by industrial processes (such as chemical and petrochemical processing, food, pharmaceutical, energetics) and processes and applications in the municipal sphere (such as waste incinerators, heating plants, laundries, hospitals, server rooms) are priorities and challenges. This Special Issue focuses on heat exchangers as key and essential equipment for the practical realization of these challenges. The purpose of this Special Issue is to outline the latest insights and innovative and/or enhanced solutions from the design, production, operation, and maintenance points of view of heat exchangers in different applications of effective waste heat utilization.

> **Petr Stehlik, Zdenek Jegla** *Special Issue Editors*

## *Article* **Heat Transfer Characteristics of Heat Exchangers for Waste Heat Recovery from a Billet Casting Process**

#### **Ju O Kang and Sung Chul Kim \***

School of Mechanical Engineering, Yeungnam University, 280 Daehak-ro, Gyeongsan-si, Gyeongbuk 712-749, Korea

**\*** Correspondence: sungkim@ynu.ac.kr; Tel.: +82-53-810-2572; Fax: +82-53810-4627

Received: 21 June 2019; Accepted: 12 July 2019; Published: 15 July 2019

**Abstract:** The application of the thermoelectric generator (TEG) system to various industrial facilities has been explored to reduce greenhouse gas emissions and improve the efficiency of such industrial facilities. In this study, numerical analysis was conducted according to the types and geometry of heat exchangers and manufacture process conditions to recover waste heat from a billet casting process using the TEG system. The total heat absorption increased by up to 10.0% depending on the geometry of the heat exchanger. Under natural convection conditions, the total heat absorption increased by up to 45.5%. As the minimum temperature increased, the effective area increased by five times. When a copper heat exchanger of direct conduction type was used, the difference between the maximum and minimum temperatures was significantly reduced compared to when a stainless steel heat exchanger was used. This confirmed that the copper heat exchanger is more favorable for securing a uniform heat exchanger temperature. A prototype TEG system, including a thermosyphon heat exchanger, was installed and a maximum power of 8.0 W and power density of 740 W/m<sup>2</sup> was achieved at a hot side temperature of 130 ◦C. The results suggest the possibility of recovering waste heat from billet casting processes.

**Keywords:** waste heat recovery; cylindrical shape heat source; thermoelectric generator; radiative heat exchanger; numerical analysis; industrial experiment

#### **1. Introduction**

The global annual energy consumption is 474.1 PJ, of which 52% is discharged as waste heat in exhaust gas and effluents. Particularly in industrial areas, 22% of the total energy is being used annually but only half the amount is being used efficiently [1]. The low efficiency of industrial facilities affects greenhouse gas emissions, which are the main cause of global warming; therefore, the recovery of waste heat is expected to contribute to increasing economic efficiency and decreasing greenhouse gas emissions by improving the efficiency of facilities [2]. For this purpose, waste heat recovery technologies, including heat pumps, boilers, refrigeration cycles, and heat exchangers, have been developed. Among them, the application of the heat pipe and heat exchanger with a thermoelectric generator (TEG) and organic Rankine cycle (ORC) have been studied recently. The ORC system is a waste heat recovery technology that uses an organic refrigerant with low evaporation temperature to produce electric energy through the turbine. Peris et al. [3] installed an ORC system in a ceramic furnace and produced 21.8 kW of power with a maximum efficiency of 12.8% at 287 ◦C heat source temperature using exhaust gases. Ramirez et al. [4] installed an ORC system on a steel mill and obtained a 21.7% efficiency at a heat source of 529.6 ◦C. However, the ORC-based waste heat recovery system includes a heat exchanger, pump, and turbine, which limits its application due to its wide installation area. TEG is an eco-friendly energy conversion device that uses the Seebeck effect, in which current is generated from the temperature difference between high-temperature and low-temperature

sections at both ends of two semiconductors. Although its current energy conversion efficiency is only 2–5%, a conversion efficiency of 15% or higher is expected to be attainable through the development of material technologies [5,6]. The application of the TEG system to internal combustion engines for waste heat recovery from high-temperature exhaust gas has already been investigated in the transport field. Wang et al. [7] fabricated a prototype TEG system to recover waste heat from the exhaust gas of hybrid vehicles and predicted potential improvement in fuel efficiency by up to 3.6% through experiments on power generation capacity and energy optimization. Moreover, Eddine et al. [8] performed an experiment using a TEG system for recovering waste heat from ship exhaust gas and predicted that maximum energy conversion efficiencies of 0.9% and 1.3% can be achieved with the use of Bi2Te3 at a maximum operating temperature of 250.0 ◦C and Si80Ge20 at 300.0 ◦C, respectively. In the industrial field, large amounts of energy can be generated using waste heat recovery systems even with low efficiency because there are various waste heat generation patterns and the scale of waste heat generation is large. For this reason, some studies have investigated the application of the TEG system for recovering industrial waste heat recently. Ebling et al. [9] installed a TEG system for the cooling process of steel forging products using a copper absorber plate to recover radiative waste heat, and produced 388 W of electric power at a heat source temperature of 1300 ◦C with an efficiency of 2.6%. Yazawa et al. [10] conducted research on the application of the TEG system to a glass melting furnace to recover waste heat. They analyzed various design parameters under the assumption that the TEG module of a direct conduction type was installed on the wall of the furnace, and predicted that a maximum electric power of 55.6 kW can be produced.

For the recovery of industrial waste heat, the TEG system has been applied to processes with high heat source temperatures (600.0 ◦C or higher), such as steel and glass manufacturing processes. In the case of non-ferrous metal industries, such as aluminum, brass, and copper; however, the temperature of the heat source at which waste heat is generated is in the mid-range, i.e., 300.0–600.0 ◦C, and thus the power generated by TEG is significantly reduced. Therefore, for processes with mid-range heat source temperature, the energy conversion efficiency of the TEG module as well as the heat transfer efficiency from the heat source to the TEG module must be considered in order to produce a large amount of energy. In addition, it is necessary to conduct research on the geometry and types of heat exchangers as important design parameters of the TEG system according to the waste heat generation patterns. In particular, for radiative energy, plate-type heat exchangers have been mostly used in the TEG system for heat sources located on horizontal surfaces, such as conveyor belts [9,11]. In actual processes, however, waste heat is generated from more complex geometry and plate-type heat exchangers cannot recover heat efficiently. Thus, some studies have been conducted to recover waste heat using heat pipes. A heat pipe consists of an evaporation unit, an insulation unit, and a condensation unit. The heat exchanger transfers heat through phase change via evaporation and condensation of the internal working fluid. The applicable temperature range of the heat source is wide and the heat transfer efficiency is high. For the application of heat pipes to industrial facilities, various geometries and performance parameters have been investigated through experiments and numerical analysis [12–14]. Regarding studies on the use of heat pipes in the TEG system, Huang et al. [15] confirmed a maximum output of 6.25 W in a TEG module using a loop heat pipe, and Wu et al. [16] showed that the use of heat pipes increased the TEG power generation efficiency by up to 50% in a comparative experiment using heat pipes and plates as heat spreaders. Heat pipe pyroelectric system with energy harvesting technology is similar to heat pipe TEG system. However, pyroelectric devices are more suitable where the temperature of the heat source changes over time [17,18]. The TEG system seems to be suitable because the heat source of this study is less temperature change with time.

This study attempted to minimize heat loss in a brass billet cooling process by applying a heat exchanger with the geometry of the area surrounding the billet such that heat could be more efficiently absorbed than plate-type heat exchangers. The TEG system was applied for the first time to the billet cooling process with a cylindrical heat source that generates waste heat of 600.0 ◦C or less. For the application of the heat exchanger, the difference in heat absorption and surface temperature

distribution were analyzed according to the geometry of the heat exchanger and process conditions through numerical analysis. In addition, the reliability of the numerical analysis of the target heat source was verified through experiments. A prototype TEG system with a thermosyphon heat pipe was installed in the billet casting cooling process. Experiments on the output power and measurement of the temperatures of the high-temperature and low-temperature sections of the heat pipe and TEG module were performed.

#### **2. Experimental Setup and Numerical Analysis**

#### *2.1. Experimental Setup*

In this study, the brass billet cooling process in Gwangmyeong-si, Korea, was set as the target of the TEG system for waste heat recovery. The side exposed to the ambient after the billet continuous casting process was set as the target position, as shown in Figure 1. The billet diameter was 240 mm and the billet length exposed to the outside was 600 mm. An air injector was installed at the outlet of the casting machine to cool the billet surface, and the billet entered the water injector through an air cooling section. Heat loss from the billet at the waste heat recovery position was caused by natural convection due to the high-temperature surface, forced convection due to the air injector, and radiation around the billet. Regarding the heat exchanger for waste heat recovery from the billet, a non-contact heat exchanger is required because the billet moves after casting. As the flow around the billet is also not constant due to the air injector, the waste heat is recovered through radiative heat exchange.

**Figure 1.** Brass billet continuous casting process.

A temperature measurement experiment was performed to verify the analysis. Figure 2 shows the schematic diagram of the experimental apparatus. An air injector was installed at the outlet of the casting machine to inject air towards the billet and a centrifugal fan capable of providing a volumetric flow rate of up to 63 m3/min was connected. A SUS 316 absorber plate was installed at 20 mm from the billet to measure the radiative heat temperature between the water injector and the air injector. The length, height, and thickness of the plate were 500, 290, and 5 mm, respectively. For temperature measurement, three thermocouples were attached to the center of the specimen at 105 mm height intervals. An additional thermocouple was attached to measure the temperature distribution on the billet surface. As the billet was cast and moved, the thermocouples also moved and the surface temperatures of the billet were measured at a distance from the outlet of the casting machine. Table 1 shows the specifications of the experimental apparatus. In the experiment, the temperature of the specimen was judged to reach a steady state when temperature changes in the three thermocouples attached to the absorber plate were ±3.0 ◦C or less considering the errors of the thermocouples and the data logger. The specimen reached a steady state in approximately 30 min. The average temperatures

at each position were calculated by collecting temperature data for more than 10 min and they were used as result data.

**Figure 2.** Schematic diagram of the experimental apparatus for analysis verification.


**Table 1.** Specifications of the experimental apparatus for analysis verification.

Figure 3a shows the schematic diagram of the radiative waste heat recovery TEG system for the experiment on output power. Radiative waste heat from the billet was recovered using the SUS 316 heat pipe. The working fluid of the heat pipe evaporated on the surface facing the billet and condensed on the top plate with the TEG module. Three K-type thermocouples were attached to measure the temperatures of the heat pipe condensation section as well as the high-temperature and low-temperature sections of the TEG module. The charging controller of the maximum power point tracking (MPPT) type was used to charge the battery with electrical energy generated by the TEG module. The battery charging power from the MPPT controller and the temperatures measured from the heat pipe and TEG module were transferred to the monitoring system through the main controller to enable real-time data collection and verification. Bi2Te3 was used in the installed TEG module. This material affords a maximum conversion efficiency of 3.4% when the maximum operating temperature of the high-temperature section is 300.0 ◦C and the temperature difference between both ends is 250 ◦C. Three modules with width, depth, and height of 60, 60, and 3.7 mm, respectively, were used. Figure 3b is a schematic diagram of the heat pipe composed of a single wick. The radius of the heat absorption plate of the wick, Ra, is 130 mm, the radius of the heat insulation plate, Ri, is 160 mm, the height of the wick portion, H, is 160 mm, and the total axial length is 500 mm. The filling ratio of the internal working fluid was set to 30%. In the heat pipe, DOWTHERM-A, a mixture of C12H10 and C12H10O, was used as the working fluid, with the operating temperature ranging from 15.0 to 400.0 ◦C. Table 2 shows all the specifications of the TEG system.

**Figure 3.** Schematic diagrams of the thermoelectric generator (TEG) system and heat pipe (**a**) TEG system; (**b**) heat pipe.


**Table 2.** Specifications of the TEG system.

#### *2.2. Numerical Analysis Model and Boundary Conditions*

The heat exchanger for waste heat recovery was modeled in three geometrical shapes, as shown in Figure 4. The first model is the n-type shown in Figure 4a. It targeted at heat transfer to the TEG module installed at the top using the heat pipe. Flat plates were applied to the sides such that the fixing jig could be installed at the bottom. The second model is the o-type shown in Figure 4b. It is a model for comparing differences in heat absorption and temperature distribution according to the heat exchanger geometry. The flat plates on the side of the n-type were bent towards the billet such that radiative waste heat could be absorbed effectively. The third model is the d-type shown in Figure 4c. It targeted at heat transfer through direct conduction after the recovery of radiative waste heat. The TEG module can be installed both on the side and at the top, unlike the n-type and o-type. SUS 316 was used for the n-type and o-type considering the mechanical durability of the heat pipe, and they had the same heat exchange area of 0.4 m2. Copper was used for the d-type for effective heat conduction through direct conduction. Its heat exchange area is the same as those of the n-type and o-type. Table 3 shows the specifications of all heat exchanger models.

**Figure 4.** Geometry of the waste heat exchanger absorber plate. (**a**) n-Type; (**b**) o-Type; (**c**) d-Type.

**Table 3.** Specifications of heat exchanger model.


ANSYS Fluent (ver. 19.2), a commercial software program, was used for heat flow analysis in this study. Table 4 shows the analysis method. The analysis was conducted on the domain, including the billet, air injector, and heat exchanger. Air was used as an incompressible ideal gas under steady-state, incompressible flow, and turbulent conditions. Moreover, the convergence of the residual was increased using the coupled algorithm and pseudo transient approach. The realizable k-epsilon model was used to simulate the turbulent flow around the billet and heat exchanger caused by the air injector. The scalable wall function was used to compensate for the grid quality degradation occurring from the narrow space between the billet and heat exchanger. The pressure was set as the body force weighted scheme considering the occurrence of the natural convection effect on the surface of the billet. The radiative heat exchange between the billet and heat exchanger is the main heat exchange method in this analysis and the surface to surface (S2S) model was used as an analytical model. The S2S model calculates radiative heat exchange using the view factor and assumes the heat exchange medium as a gray body and a diffusive surface. The absorption coefficient and the scattering coefficient could be neglected in the model, and the model was suitable for the analysis because there was no permeable medium in the analytical domain of this study [19]. For a somewhat larger scale with geometry similar to that of this study, Mirhosseini et al. [20] conducted 2D numerical analysis on the radiative heat exchange of an arc absorber around a cylindrical cement kiln using an S2S model. Table 5 shows the boundary conditions used for heat flow analysis. Emissivity, a surface property, was set to 0.7 for the brass billet, 0.4 for SUS 316, and 0.5 for copper [21–23]. Figure 5 shows the results of measuring the billet surface temperature. Considering the process, the average temperature of 530 ◦C at the position where the heat exchanger can be installed was set as the surface temperature of the billet. To compare the heat absorption and temperature distribution under forced convection and natural convection conditions, blower flow rates of 0 and 0.18 kg/s were added to each geometrical shape.


**Table 4.** Analysis method.

**Table 5.** Analysis of boundary conditions.


**Figure 5.** Temperature distribution on the billet surface.

#### **3. Results and Discussion**

Approximately 15 million grids were used to ensure grid independence of the analytical model. Convergence was determined when the residual value of energy dropped below 10−<sup>6</sup> and the remaining residual value dropped below 10<sup>−</sup>4. The reliability of the analytical model was confirmed through the experiment for analysis verification and the maximum error of 2.4%. Table 6 shows temperatures and errors at each measuring point obtained from the numerical analysis and the experiment.

**Table 6.** Temperature difference between experimental and numerical analysis.


Figure 6 shows the temperature distribution according to the geometry of the SUS 316 heat exchanger. Under forced convection conditions, the high-temperature section on the heat exchanger surface was formed at the front for all geometrical shapes. Figure 6a shows the surface temperature distribution of the n-type under the forced convection condition. The maximum temperature was 272.5 ◦C, whereas the minimum temperature was 75.5 ◦C. Figure 6c shows the surface temperature distribution of the o-type under the same condition. The maximum temperature was 274.1 ◦C, whereas the minimum temperature was 94.1 ◦C. The difference in the maximum temperature between the n-type and o-type was small (1.6 ◦C), but the difference in the minimum temperature was relatively large (18.6 ◦C) with the o-type exhibiting higher minimum temperature. The asymmetric temperature distribution of the symmetric heat exchanger could be predicted because the position of the injection hole on the air injector is not symmetrical. Figure 6b and d show the heat exchanger surface temperature distributions of the n-type and o-type, respectively, under the natural convection condition. The maximum temperature was 348.4 ◦C for the n-type and 346.4 ◦C for the o-type, showing an increase of up to 75.9 ◦C compared to the maximum temperatures obtained under the forced convection condition. This can be attributed to the increase in the temperature of the heat exchanger with the weakening of the effect of forced convection, reducing heat loss on the surface. Moreover, the minimum temperature was 198.1 ◦C for the n-type and 286.3 ◦C for the o-type, showing an increase of up to 192.2 ◦C. The temperature difference between the n-type and o-type was obvious under the natural convection condition. Figure 6e shows the surface temperature distribution of the d-type under the forced convection condition. The maximum temperature was 221.6 ◦C, whereas the minimum temperature was 148.4 ◦C. Figure 6f shows the surface temperature distribution of the d-type under the natural convection condition. The maximum and minimum temperatures were 280.2 ◦C and 243.0 ◦C, respectively. While the difference between the maximum and minimum temperatures was up to 197 ◦C for the n-type and o-type, that of the d-type was up to 73.2 ◦C, which is significantly low. The difference between the maximum and minimum temperatures of the d-type appeared to decrease because of the high thermal conductivity of copper. For the d-type, a uniform temperature distribution is required for the heat exchanger because heat is transferred to the TEG module through direct conduction. The analysis results showed that the use of the copper heat exchanger is appropriate because it exhibits a relatively uniform temperature distribution.

As shown in Figure 6g,h, a convection layer was formed between the billet and the heat exchanger. The high-temperature section was concentrated at the front because the cooling air that entered the rear end of the heat exchanger close to the air injector was slowly heated as it passed between the billet and the heat exchanger. Figure 7 shows the temperature distribution along the heat exchanger circumference. The maximum temperature was not significantly different depending on the geometry, but the temperature significantly increased in the o-type due to the bent sides.

Although stainless steel was selected considering the mechanical durability of the heat pipe, the temperature difference between the high-temperature and low-temperature sections was large due to its low thermal conductivity. To determine the area where the working fluid can evaporate and heat absorbed when heat pipes are used, the effective area and waste heat absorption were analyzed according to the geometrical and flow conditions for the basic design of the heat exchanger. The effective area was defined as the area where the temperature of the heat absorbing section was 250 ◦C or higher considering that the maximum operating temperature of the high-temperature section of the TEG module was 300 ◦C. The effective area and waste heat absorption, according to the geometry and flow conditions, are shown in Figure 8. Under the forced convection condition, a total of 1.0 kW waste heat was absorbed by the heat exchanger for the n-type. For the o-type under the same condition, the radiative heat energy absorbed from the billet increased due to the bent sides and 1.1 kW waste heat was absorbed. Under the natural convection condition, the energy loss attributable to convection decreased and thus the n-type and o-type absorbed 1.5 kW and 1.6 kW waste heat, respectively. The effective area did not significantly differ depending on the geometry, and it significantly increased under the natural convection condition. Thus, most of the area exhibited 250 ◦C or higher for both the n-type and o-type, indicating effective heat transfer. The numerical analysis results confirmed that the geometry contributed to 10% improvement of waste heat absorption but it

did not significantly affect the maximum temperature and effective area. Under the natural convection condition, however, waste heat absorption was increased by up to 45.5% and the effective area was increased up to five times.

**Figure 6.** Temperature distribution with different heat exchanger shapes and air flow and velocity distribution between the billet and the heat exchanger (**a**) n-Type forced convection; (**b**) n-Type natural convection; (**c**) o-Type forced convection; (**d**) o-Type natural convection; (**e**) d-Type forced convection; (**f**) d-Type natural convection; (**g**) n-Type forced convection; (**h**) o-Type forced convection.

**Figure 7.** Surface temperature distribution along the heat exchanger circumference.

**Figure 8.** Heat absorption with different heat exchanger shapes and air flow.

The TEG system with an n-type prototype thermosyphon was tested under the forced convection condition. Table 7 shows the temperatures and output power of the TEG system in a steady state. In the actual experiment, the temperature of the condensation section of the heat exchanger was 263.3 ◦C, meeting the proper operating temperature of the TEG module. In the numerical analysis, the temperature of the high-temperature section at the top of the n-type was similar to that at the top of the thermosyphon with an error of 3.5%. The temperature of the high-temperature section of the TEG module was 130 ◦C, which was also appropriate as it was in the operating temperature range. The three TEG modules provided an output power of 8.0 W when the temperature of the low-temperature section was 49.8 ◦C and the temperature difference between both ends was 80.2 ◦C. The temperature drop that occurred at the condensation section of the heat exchanger and the high-temperature section of the TEG module; however, was a major obstacle to the system output. This temperature drop can be attributed to the insufficient waste heat recovery and low thermal conductivity of the TEG system. Regarding the evaporation section, as confirmed by the numerical analysis, the evaporation of the working fluid was not sufficient because the thermosyphon evaporation section could not secure sufficient temperature and heat absorption, thereby, degrading the performance of the heat exchanger itself. As confirmed through the numerical analysis, heat absorption was increased by 50.0% under the natural convection condition. Therefore, if the air injector is removed from the TEG system as a means of securing output and an experiment is performed under the natural convection condition, the output of the TEG system

is expected to increase. The heat exchanger output can also be increased by increasing the emissivity of the thermosyphon through surface treatment, such as surface blackening using carbon coating.


**Table 7.** Temperatures and output power of the TEG system.

#### **4. Conclusions**

In this study, the heat absorption and temperature distribution of a TEG system for recovering waste heat from a brass billet casting process were analyzed according to the heat exchanger type, geometry, and process condition using numerical analysis. Heat absorption could be increased by up to 10.0% depending on the geometrical variables, but the maximum temperature did not significantly change. Moreover, when the forced convection condition was changed to the natural convection condition, of which the blower flow rate was 0 kg/s, heat absorption could be increased by up to 45.5%, the maximum temperature was also increased by 75.9 ◦C, and the effective area increased five times, thereby, exhibiting more effective heat transfer. The heat exchanger of a direct conduction type also met the proper operating temperature of the TEG module and exhibited a uniform temperature distribution, thereby, suggesting its feasibility for practical application. Moreover, the experiment on the output power of the TEG system with the n-type thermosyphon under the forced convection condition confirmed the output power of 8.0 W.

The results of this study confirmed the possibility of recovering waste heat from industrial facilities. The application of the TEG system is expected to improve the efficiency of facilities and contribute to the selection of heat exchangers for the TEG system according to the type of waste heat generated. The system outputted only 8 W during the actual process because the heat transfer between the heat exchanger and TEG module has not been optimized yet. In the future, for the optimal design of heat exchangers and improvement of the output power, various types of heat exchangers and heat exchangers with coated surfaces for improving heat transfer performance will be applied to the TEG system and experiments as well as numerical analysis on power generation performance and economic assessment of its applicability will be conducted.

**Author Contributions:** J.O.K. performed the experimental/numerical analysis and drafted the manuscript. S.C.K. organized the overall evaluation and reviewed the manuscript.

**Funding:** This work was conducted as a part of the Energy Technology Development project sponsored by the Ministry of Trade, Industry and Energy (No. 20172010000760).

**Acknowledgments:** The support received from the Livingcare company is greatly appreciated.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **Nomenclatures**


#### **Subscripts**


#### **References**


© 2019 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 (http://creativecommons.org/licenses/by/4.0/).

## *Article* **Modelling and Evaluation of Waste Heat Recovery Systems in the Case of a Heavy-Duty Diesel Engine**

## **Amin Mahmoudzadeh Andwari 1,2,\*, Apostolos Pesyridis 1, Vahid Esfahanian 2, Ali Salavati-Zadeh <sup>3</sup> and Alireza Hajialimohammadi <sup>4</sup>**


Received: 3 March 2019; Accepted: 26 March 2019; Published: 11 April 2019

**Abstract:** In the present study, the effects of Organic Rankine Cycle (ORC) and turbo-compound (T/C) system integration on a heavy-duty diesel engine (HDDE) is investigated. An inline six-cylinder turbocharged 11.5 liter compression ignition (CI) engine employing two waste heat recovery (WHR) strategies is modelled, simulated, and analyzed through a 1-D engine code called GT-Power. The WHR systems are evaluated by their ability to utilize the exhaust excess energy at the downstream of the primary turbocharger turbine, resulting in brake specific fuel consumption (BSFC) reduction. This excess energy is dependent on the mass flow rate and the temperature of engine exhaust gas. However, this energy varies with engine operational conditions, such as speed, load, etc. Therefore, the investigation is carried out at six engine major operating conditions consisting engine idling, minimum BFSC, part load, maximum torque, maximum power, and maximum exhaust flow rate. The results for the ORC and T/C systems indicated a 4.8% and 2.3% total average reduction in BSFC and also maximum thermal efficiencies of 8% and 10%, respectively. Unlike the ORC system, the T/C system was modelled as a secondary turbine arrangement, instead of an independent unit. This in turn deteriorated BSFC by 5.5%, mostly during low speed operation, due to the increased exhaust backpressure. It was further concluded that the T/C system performed superiorly to the ORC counterpart during top end engine speeds, however. The ORC presented a balanced and consistent operation across the engines speed and load range.

**Keywords:** waste heat recovery; Organic Rankine Cycle; turbo-compound; brake specific fuel consumption; engine thermal efficiency

#### **1. Introduction**

In recent years, the growing worldwide population and industrial development has seen an equally increasing demand in energy. The internal combustion engine (ICE) has, by far, grown to be the most popular means of transport since the second half of the 20th century. Unfortunately, a typical ICE will only manage to convert approximately 30–35% of the total provided chemical energy into effective mechanical work, as illustrated in Figure 1 [1–6].

**Figure 1.** Variation total fuel energy content in ICE.

As with any system that produces work in real life, it is very challenging to achieve an adiabatic process during which thermal expansion occurs. Even in present ICEs, approximately 60–70% of the energy discharged by the fuel is wasted predominately in the form of heat [7–13]. During the last two decades, typical engine specific power output ratios are in the region of 1.5. At the same time, emission levels have reached a factor of 10. This is mainly due to the restrictions imposed by EURO, forcing technology to progress in the production of cleaner and more efficient ICE [14–20]. Even a conventional turbocharger only takes advantage of a portion of the exhaust gas energy in the shape of kinetics and pressure, which constitutes a fraction of the energy losses in a naturally aspirated engine. The biggest percentage is heat transfer and exhaust gas enthalpy dissipation, which is accountable for about 50–85% of the outstanding low heating values of the utilized fuel [21–25]. As a result, presently, waste heat recovery (WHR) has been the primary concentration point of research and development departments of engine manufacturers. This is due to the considerable potential energy amount that can be recovered in the form of heat [26–28]. Modern WHR systems amongst others include the following:


All of the above can recover a segment of the exhaust gas energy and can subsequently enhance the engines' thermal efficiency. They all operate on relatively similar thermodynamic principles, but not all of them perform in the same fashion. The beginning of research and development on the feasibility of the Rankine cycle system as a WHR method dates back to the 1970's [29–32]. The fundamental orientation of the studies has been the thermal optimization of heavy duty diesel engines (HDDE) due to their high thermal efficiency potential, ranging, by todays' standards, between 40–45% [33–36]. This makes HDDE the most favorable candidates, as a 10–15% improvement in fuel efficiency is not uncommon for ORC applications [37–40]. The high thermal efficiency figures have also influenced industries to make use of HDDE, with power outputs of up to 600 kW for on and off-highway commercial vehicles. Typical engine displacements vary between 6–12 liter with multiple cylinder numbers and configurations, while the latest era of HDDE, utilizing high-pressure common rail direct injection systems, reach pressures of over 2500 bar. Part of the reason for this high efficiency is the use

of forced induction. Forced induction constitutes a method of increasing engine power output that dates back to the late nineteenth and early twentieth century [41–45].

Supercharging an engine can be done both by either mechanical and/or chemical means. Power output is directly proportional to the mass of fuel and air burned inside the engines' cylinders. To ncrease fuel mass delivery, one must first increase air mass intake to avoid continuous combustion of rich mixture. Additional air supplying components that are gear and belt driven are normally called superchargers or blowers. Turbocharging is when the supercharger is driven by the engines' hot exhaust gases. One of the reasons forced induction became so popular over the past decades, other than the significant performance gains, is because it was a clever way of benefiting from the energy of burned exhaust gas residue that would otherwise be wasted. As a result, the automotive industry perceived it to be an acceptable mix of cost, performance, fuel economy, and reliability [40,46–48]. In turbocharged engines, heat transfer is a very complex entity, which greatly affects turbocharger performance, efficiency, and selection. As exhaust gas flows through the turbine, the turbine housing absorbs a sizable percentage of the total enthalpy by forced convection, due to the temperature difference between the walls and exhaust gases. This heat is then lost to the environment by means of radiation. Heat transfers by forced convection are also evident on the turbine wheel blades, shaft, and, subsequently, on the lubricant because of the exhaust gas expansion [49–54]. On the other hand, the compressor side acts as a heat sink and is subject to heat conduction derived through the bearing and/or turbine housings as well as the engine itself. This heat flux is inbound and may affect the temperature and pressure of the inducted air at low rotating speeds and compression ratios [55,56]. It was discovered that despite past research on WHR using primarily ORC and T/C systems either of mechanical or electrical nature, there is limited information on the comparison of the two systems, which are assessed in this research [1,36,57–59]. This study aims to identify which of the two adopted WHR methods (ORC and T/C) performs in a better manner in terms of improving engine thermal efficiency and BSFC by means of 1-D engine simulation software.

#### *1.1. Waste Heat Recovery*

Almost 60% to 70% of the chemical energy provided to the engine in the form of fuel is wasted, mostly from rejected heat. In a system that carries out work, heat loses are inevitable due to the first law of thermodynamics. The most arbitrary heat loses are that of the exhaust and cooling systems, but losses can also be on behalf of pumping losses, internal friction, drivetrain slippage, and other accessories. In fact, during in-town driving a typical vehicle will utilize, on average, only about 13% of the actual fuel energy to propel forward. To put it into context, a diagram presenting the energy pathway is located in Figure 2. Fortunately, developed technologies allow the conversion of a percentage of waste heat back to usable energy via several harvesting systems. Automobiles and especially heavy-duty commercial on and off highway vehicles have under scrutiny been lately [20,60,61].

**Figure 2.** Typical city driving energy pathways.

In the US alone, tractor trailers, delivery vans, garbage trucks, and more are expected to reduce 25% of their exhaust emissions by 2027, with a potential to avoid up to 1.1 billion metric tons of carbon dioxide emissions. Many automotive industries that design and manufacture WHR systems seem to avoid mentioning the disadvantages of these promising devices. For example, they tend to incur an increase in exhaust backpressure, which has a direct impact on engine fuel consumption. Not many preliminary studies can be found that take exhaust backpressure fluctuation caused by the heat exchanger or turbine into account. As a result, the comprehension data, in consideration of designing and optimizing the components, is very limited. Another negative aspect is that integrated WHR units come alongside a weight penalty. The additional inertia will inevitably cause further fuel consumption and, subsequently, an increase in BSFC. These drawbacks are easy to neglect when heat harvesting percentages and thermal efficiencies of WHR means as well as the power unit are the center of attraction of the investigation [20,62,63].

#### *1.2. Steam*/*Organic Rankine Cycle*

The Rankine cycle (RC) is considered an ideal cycle for vapor power plants. It is comprised by four primary components, a pump, a boiler/evaporator, a power expansion turbine, and a condenser. The pump begins the cycle by pumping the working fluid through the system. The evaporator or boiler applies the recovered waste heat on the fluid, thus raising its temperature and pressure, creating (in some cases) superheated steam. The fluid is then expanded in the turbine utilizing the built up temperature and pressure by generating power through a shaft. The process is continued by the condenser, which condenses the vapor back to liquid form ready to be pumped again for another cycle [6,64,65].

The actual and ideal thermodynamic evaluation of the Rankine cycle operations can be slightly different from each other. In an ideal RC system, the compression and expansion processes in the pump and turbine, respectively, are considered isentropic. In an identical trend, the heat addition and heat rejection processes in the evaporator and condenser, respectively, are regarded as constant. This transliterates these processes as internally reversible. A reversible process is the process that can be reversed without leaving any traces to the environment or surroundings. Since the pump, evaporator, turbine, and condenser are all steady-flow components, the RC can be analyzed as a closed loop, steady-flow process following the steady-flow energy equation. This in turn implies that no heat engine will have a thermal efficiency of 100% because there must be a low temperature sink for the heat to be transferred to. In addition, the actual RC system suffers from losses throughout the systems cyclic function such as friction losses, piping losses, and heat transfer to surroundings. All of these losses cause an irreversible increase in entropy [2,9,10,32].

#### *1.3. Organic Rankine Cycle and Internal Combustion Engines*

The only difference between the conventional RC and the organic Rankine cycle (ORC) is the substance which circulates within the system. A traditional RC is known for using H2O (water) as a working fluid. In the case of ORC, the operating fluid is an organic liquid element accompanied by a greater molecular mass and reduced boiling point compared to that of H2O. These alternative working fluids can demonstrate a number of thermodynamic benefits as they allow the system to operate by using downsized temperatures. Consequently, the operation of the ORC results in greatly reduced thermal efficiency readings owning to the lower temperature transactions. However, this also has a positive impact on the total operational cost, as far less heat energy is required to produce a given power. The promising potential of low and moderate heat operation is the primary reason the ORC has grown to be popular amongst automobile industry research and development departments.

This WHR method, however, is not limited to ICE applications, but is also utilized by a number of other heat rejecting machinery, such as geothermal plans, solar thermal systems, and biomass plants, as well as industrial processes. The unstable transient and remarkably variable operational profiles of automotive vehicles make it more demanding to implement ORC systems and, therefore, the technology is expected to hit the market around 2020. This is mainly due to the non-existent control methods and instruments, which are of great necessity for the safety, performance, reliability, and durability of the power unit and ORC. Correspondingly, ORC patents are primarily developed and promoted for immobilized power production units as well as marine purposes [41,44,48,62]. In an investigation of an ORC system for a heavy-duty truck application and a passenger-car application, a 3.4% estimated reduction in fuel consumption was obtained. This was the result of a turbine efficiency of 58%, after the custom blade was designed on CFD software specifically for the employed ORC system. During the analysis of a diesel-Rankine cycle combination in a different HDDE case [58], it was concluded that during full load conditions (BMEP = 2 Mpa) a BSFC reduction of 2.6–3% was possible. The values of the temperature, pressure and, working fluid flow rate were all estimated by the thermodynamic characteristics of the multifarious utilized substances, namely, water (H2O), methanol (MeOH), toluene (PhMe), pentafluoropropane (R245fa), and tetrafluoroethane (R134a). Moreover, the performance of an ORC system by integrating the WHR system on a 15 L diesel engine was investigated. The simulation assessment was performed by gathering the turbine shaft power from the exhaust enthalpy during steady-state operation. The use of dual finned heat exchangers with identical dimensions and properties was implemented in a parallel sequence after the turbocharger turbine to collect waste energy. The total power output of the engine was increased by an estimated 5%, while the engines' pumping losses were kept at a maximum total of 4 kW. It was also supported that much research based on WHR seem to neglect circumstantial disadvantages. Investigations based on ORC are no exception, with a couple of crucial unmentioned performance characteristics. These include the effect of refrigerant flow rate on ORC performance as well as the effects of pressure drops through the heat exchangers, with the resulting parasitic flow-work losses, two negative aspects which are overcompensated considering that the average theoretical integrated vehicle ORC system yields an increase in thermal efficiency of about 6–15%. As far as ORC fitment is concerned, depending on packaging and weight limitations of the given vehicle, the ORC recovery method can include multiple supplementary components. Typical layout implementations of ORC systems to HDDE can include twin parallel or in-series evaporators, individual for the exhaust and EGR valve. Furthermore, the introduction of a recuperator has found its way into the system due to the possibility of ORC efficiency increments. It is typically positioned between the turbine and the condenser and its functionality is to recuperate some of the heat before it is released to the heat sink by the condenser. Preheating the working substance with the aid of a charge air cooler (CAC) also constitutes an investigation possibility. Other researchers suggest replacing current engine block cooling techniques with ORC working fluids to take advantage of the additional waste heat and improve power regeneration. On the other hand, all of these efforts and aspects tend to increase the systems complexity rather than provide considerable ORC gains. Hence, a straightforward simplistic ORC composition is a more appealing solution for vehicle integration than most [6,40,44,65,66].

#### *1.4. Engine Turbocompounding*

In a conventional turbocharger, the engines' exhaust gas heat and airflow energy is harvested by a turbine, which is connected to a compressor through a common shaft. In the compressor side, air is induced and pressurized in the intake manifold, which increases the total power output of the engine with a small penalty on exhaust backpressure. A Turbo-compounding (T/C) system operates in a similar fashion, with the only difference being that there is no compressor at the end of the turbine shaft. The T/C turbine would be typically placed at the outlet of the primary turbine, therefore being driven by the leftover energy, translating it into a torque. T/C is a potentially prosperous WHR method, which can either be of a mechanical or electrical nature. A T/C system is a relatively less complicated arrangement compared to the ORC and this could potentially result in a lower unit production cost and a lighter component all together. On the other hand, one of the primary disadvantages of T/C implementation is the increment of engine backpressure and pumping losses, even more so than the ORC systems' heat exchanger. Exhaust backpressure is directly proportional to cylinder pumping loses. Hence, during meager engine speeds and loads, the total engine brake power output is prone to suffer. As mentioned before, the power produced from the turbine can be manipulated either electrically or mechanically [10,14,30,67].

#### 1.4.1. Electrical Turbo-Compounding

In an electric T/C system, an alternator/stator converts the turbines' rotational shaft power into electrical power. This electricity then either returns to the main battery to be stored for later usage or is immediately effective to operate various engine/vehicle components such as the starter motor, headlights, etc. Another possibility is the integration of an electrical compressor acting as a supercharger to assist the vehicles acceleration during low engine speeds where turbo-lag is yet to be overcome. Furthermore, the function of an electric motor directly mounted to the engines' crankshaft can also assist engine operation. Apart from throttle response, the techniques described can additionally improve fuel economy. In fact, in an electric T/C system, the estimated indication of reduction in fuel consumption ended up being a maximum of 10%. In addition, the strategy of the motor to crankshaft scheme seemed to enhance drivability and engine flexibility during transient periods. This in turn decreased exhaust gas emissions for an altogether greener engine activity. One of the other main advantages of electrical turbo-compounding is the space saving characteristics. Whether it is implemented as an integrated unit or a separate turbo-generator, it is a neat and tightly packaged component compared to the mechanical T/C system. On the other hand, the downside in the fitment of the electric T/C system is that it would generally require modifications to the existing turbomachinery [2,6,13,49,64]. This means that the already implemented turbocharger would have to be customized to incorporate the addition of a stator and rotor doublet in-between the turbine and the compressor impellers. However, there is the possibility that the electrical system is mounted on a separate turbine downstream of the main power turbine, also called a turbo-generator. This would diminish the need for existing turbo modifications because it is an independent and standalone unit.

#### 1.4.2. Mechanical Turbo-Compounding

Similarly, the mechanical T/C system operates again by the addition of a secondary power turbine mounted sequentially after the principal turbine, to scavenge the surplus energy. The generated power is afterwards transmitted, via a shaft, through a gearbox unit followed by a mechanical coupling to the power units' crankshaft. It is generally a low volume and production cost system thus it is fundamentally applicable for medium and heavy-duty diesel power units. These leading HDDE manufacturers have been investigating the effects of different T/C arrangements with satisfactory results, namely a typical reduction in BSFC on a scale ranging between 3–6%. Investigations proved that the total improvement in incremental fuel consumption strictly due to the turbo-compounding action was a 4.2% to 5.3% estimate depended upon the terrain or mission load factor. Incorporating a mechanical T/C system in favor of an 11 Liter 6 cylinder turbocharged diesel engine resulted in a total of 5% reduction in BSFC during full load operation [42,43,45].

#### **2. Engine Waste Heat Recovery System Modelling**

In this section, the methods and tools used to assess the effects of a T/C system against the effects of an integrated ORC are explained. Important performance parameters, such as BSFC, thermal efficiency, power output, and overall fuel consumption are monitored and examined in conjunction with two WHR methods. The same virtual turbocharged HDDE was utilized for both the adopted techniques in favor of result accuracy. The T/C system was regarded as a secondary turbine arrangement downstream of the primary power turbine. With the aid of 1-D computer software (GT-Power), the engine, ORC, and T/C systems are modelled and optimized via trial and error simulations.

#### *2.1. Engine Modelling and Calibration*

An in-line, 6 cylinder, 11.5 liter, and turbocharged HDDE is assumed as a base research engine. The engine's major specifications are presented in Table 1. Figure 3 illustrates the engine model

developed in GT-Power. The engine model is the same validated engine model as featured in Karvountzis-Kontakiotis et al [58], presented at the SAE World Congress. The model was further modified, in addition to ORC, to be able to simulate turbo-compounding operation.


**Table 1.** Modelled engine specifications.

**Figure 3.** Engine model in GT-Power software.

A typical modern direct injection diesel engine is capable of working over a speed range of 600 rpm to 2600 rpm. This speed range is larger than normal but allows the data presented to be used over any part of that range. Therefore, the engine model was operated through 36 different cases ranging from 850 rpm to 2600 rpm using 50 rpm increments for a better result accuracy. It was decided that the best way to approach the optimization process was to limit the number of variables. This task proved to be challenging because, apart from engine speed, the operational cases varied in terms of engine load, fuel mass flow rate, power target of injector controller, turbine speed, etc.

A solid baseline of results was achieved by running the two WHR rivals at operating points, which the engine spends most of its working time. Figure 4 presents the typical engine speed and

torque time percentage distribution for a crawler loader. The engine tested during this task is the type of engine that could be used for similar off-highway vehicles, as in Figure 4.

**Figure 4.** Typical HDDE time operational profiles.

In order to identify the ideal benchmarks in which the WHR systems will be constructed, optimized, and simulated, it is critical to run the engine simulation under real-life operating circumstances. In general, these optimum operational points are where the power unit produces certain usable benefits, such as maximum exhaust gas flow rates, lowest BSFC etc. Therefore, emphasis has been given to these points, which resemble an operation this engine is most likely to have under a typical working condition, similar to the example illustrated in Figure 4. Therefore, the WHR systems are assessed and compared for 6 dominant engine functioning points, which are determined by the calibration process. These include running at idling, ideal BSFC, part load, maximum torque output, maximum brake power output, and maximum exhaust flow rate. These are denoted by X1, X2, X3, X4, X5, and X6 respectively. As a gauge, the obtained baseline engine specifications are plotted using GT-Power 2-D graphical representations. This includes the BSFC contour map, which also indicates individual BMEP readings as well as power versus torque curves, all plotted in Figures 5 and 6 respectively. For the engine assessment, performance parameters, and comparison accuracy, both WHR methods are implemented on the same engine model.

**Figure 5.** Baseline BSFC contour map.

**Figure 6.** Baseline Power vs. Torque.

#### *2.2. ORC System Modelling*

The ORC system model is created and optimized in GT-Power software according to the predetermined aims and objectives of the study. The ORC model was kept to a minimal level and thus consisted of the following four main components: The pump, the boiler/evaporator, the expansion turbine, and the condenser. The pump and turbine elements were each coupled to a speed governor that sets the speed for each case run. These speed governors allowed turbine and pump speed variations for each operational point and, as explained later, this proved to be of significant value for the determination of the individual points' maximum performance enhancement. However, an increase in pump speed provokes an increase in work input requirement. Thus, to accomplish positive results, the energy recovered by the system will have to unavoidably reimburse the energy cost necessary for the pump operation.

The amount of energy deducted by the pump has a direct impact on the ORC system's efficiency due to the following relation. Another important factor which greatly interferes with the efficiency of the ORC system is the working fluid, so the refrigerant of type R245fa (Pentafluoropropane) is selected due to its advantageous low temperature heat recovery characteristics. The model of ORC system setup is illustrated in Figure 7. The control volume was considered to be adiabatic and, therefore, during all processes there was no heat escaping through the walls and surrounding features. This signifies that the exhaust gas pressure and temperature at the inlet and outlet of the evaporator were equal. Similarly, the coolant pressure and temperature at the inlet and outlet of the condenser are kept at an equal value. The next step is implementing the engines' turbocharger turbine outputs as ORC inputs at the heat exchanger. That includes the exhaust gas mass flow rate and temperature in the exhaust pipe downstream of the turbine for the six major running points of X1, X2, X3, X4, X5, and X6. Since the important aspect was to understand the strengths and weaknesses of each WHR system, the ORC pump and turbine speeds are altered during testing. Thus, the pump and turbine speeds were set according to literature review and benchmarking values were set, starting from 1500 rpm reaching up to 2500 rpm using 250 rpm increments. The design parameters of ORC used in the simulation are presented in Table 2.

**Figure 7.** ORC system model in GT-Power.

**Table 2.** Component design parameters of ORC using R245fa refrigerant.


#### *2.3. Turbocompound System Modelling*

Similar to the engine and ORC system model, the simulation of the T/C system model is performed with the aid of GT-Power software. On an industrial technicality level, the T/C system would have to be modelled on a separate template with the full extent of its integrated components. However, for this investigation and simplicity purposes, a simple secondary turbine is placed posterior to the primary turbine. In the pursuance of supervision and manipulation reasons, the turbine is incorporated with a rotational speed regulator as well as a signal output monitor. This model is a baseline calculation estimate and not a detailed representation of a T/C system. For example, the current T/C system model arrangement does not consider mechanical or electrical losses and thus the simulation will not represent real life expectation conclusively. The level of uncertainty is particularly higher at the developmental variables of the turbine, such as the performance map, which is a necessity for the validation of the modelled T/C system. As a result, literature review provides the essential assumptions and input specifications in order to avoid any potential errors. The model of the integrated T/C system is displayed in Figure 8.

**Figure 8.** Turbo-compound system model.

Identically to the ORC system, benchmarking and literature review were not enough to optimize the T/C system model and thus it has to trail a series of experimental procedures by incorporating variable parameters as a plot of trial and error. In general, the power generated by the turbine will not be linear nor at its peak for all operating conditions. Therefore, the models' turbine is assessed during diverse rotational speeds, ranging from 20,000 rpm to 120,000 rpm using 10,000 rpm intervals, for all six benchmark points, resembling the process followed by the ORC system model. This will allow the identification of the optimum turbine speed for each given case and, hence, achieve maximum power output for the system as a total. The exhaust backpressure is expected to rise owing to the layout of the T/C system, which is placed directly after the turbocharger. A rise in exhaust gas pressure translates to further engine pumping losses during the intake and exhaust strokes. A comparable ORC system will also increase backpressure in the evaporator and this effect has been studied closely and published in a dedicated paper [58] by the Brunel University team. The important assumption here was that the heat exchanger technology, which can be employed, would have a minimal impact on fuel consumption. This is not an unfair assumption in view of the availability of heat exchanger technologies available with minimum impact to the gas exchange process. For the T/C case however, backpressure is heavily dependent on the operative expansion ratio and efficiency of the turbine expander—two parameters which are captured in the present investigation.

#### **3. Results and Discussion**

As mentioned, an in-depth investigation is conducted to assess and compare the advantages and disadvantages of integrated T/C and ORC systems. This section will provide a thorough comprehension for the results comparison of both WHR methods. In the specific engine property, between two configurations, the BSFC value ranges from a maximum value of 298.09 g/kW·h down to a minimum of 205.87 g/kW·h. The variation in power gains caused by the different pump and turbine speeds for the ORC and T/C systems (only the latter is varied for the T/C as there is no pump in the

model) is also reviewed and explained in detail. The superior WHR method will be exclaimed by the BSFC reduction percentages and flexibility in operation.

#### *3.1. Engine Waste Heat*

The primary engine parameter, which is responsible for the performance of the WHR systems, is the available energy after the turbocharger turbine. The harvesting of waste heat is mainly depended on the accessibility of waste energy. The exhaust mass flow rate and exhaust gas temperature of the engine define the available energy for harvesting. By recording the exhaust mass flow rate and exhaust gas temperature values, one can determine the exhaust energy at the desired points of study. In general, more efficient energy gatherings are possible during top end operation. If the power unit is working at high engine speeds, the air mass inducted by the pistons increases, followed by larger amounts of fuel injected in the cylinders. Hence, the exhaust energy is enhanced due to the rise in exhaust gas mass flow rate. Therefore, it can be stated that there is a direct correlation between available exhaust energy, exhaust gas mass flow rate, engine speed, and possibly engine power output. Figure 9 represents the mass flow rate, BMEP, and engine speed, proportionality.

**Figure 9.** Variation of exhaust gas mass flow rate in accordance with engine speed and BMEP.

On the flipside, this correlation is not true for the exhaust gas temperature. In general, a high exhaust gas temperature signifies a deficient engine thermal efficiency. This is because a larger portion of energy is escaping from the combustion chamber, instead of being converted into usable mechanical work. It is observable from the contour plotting on Figure 10 that the exhaust gas temperature displays a higher temperature degree during intermediate BMEP and engine speeds.

**Figure 10.** Variation of exhaust gas temperature in accordance with engines' speed and BMEP.

The lowest operating temperature level is achieved at the 1500 rpm point and around 20 bar of BMEP. This reveals the minimum BSFC value point (X2), which means that at that specific point, the engine is operating around peak thermal efficiency. Therefore, this proclaims that the exhaust temperature will predominately be higher at reduced thermal efficiencies. However, the contrast in exhaust gas temperature between the maximum and minimum thermal efficiency points is not of significant scale (approximately 160 ◦C). As a result, this sets the exhaust residue mass flow rate as the primary responsibility factor of regenerated power volume. This can also be validated by comparing the exhaust gas temperature contour map against that of the thermal efficiency profile, Figures 10 and 11, respectively.

**Figure 11.** Variation of engine thermal efficiency in accordance with engines' speed and BMEP.

Conclusively, it can be suggested that there is a direct relation between WHR performance, engine speed, and reduced thermal efficiency profiles. This means that there is a greater potential to recover the engine's exhaust waste heat during operation at lower engine thermal efficiency.

#### *3.2. Organic Rankine Cycle System*

#### 3.2.1. ORC System Speed Variation

By simulating the ORC system using various pump and turbine speeds it is possible to maximize the positive characteristics for each of the six assessment points. Figure 12 represents the system speed variation to the engine test points.

**Figure 12.** ORC system turbine power output in relation to system speed variation.

It is understandable that the power requirement to drive the pump increases proportionally with the pump speed. However, it is noticeable that the turbine power output fluctuates as the turbine and pump speed varies. On one hand, during the first two points (X1, X2), where the engine speed and thus exhaust gas mass flow rate is mediocre, the system is productive mostly at medium to low speeds. In addition, the variation in pump and turbine speed does not seem to provoke a considerable divergence in power output across the intervals. The reason is that the power output during low engine speeds is relatively low. On the other hand, as the engine speed rises (especially at points X5 and X6) so does the exhaust mass flow rate. Therefore, the amount of waste energy available for recovery increases. This achieves a greater power acquirement per working cycle and, hence, the difference between the power output levels between the system speed intervals grows significant. Table 3 includes the total range of pump and turbine speed performance variations for the ORC system, with the maximum work input and output values highlighted in red.


**Table 3.** ORC system speed variation simulations.

#### 3.2.2. ORC System BSFC Reduction

For comparison purposes, all the calculations are conducted by utilizing the maximum values at each point. Now, BSFC is defined as the amount of fuel used per unit amount of power. By increasing power output, without increasing the fuel mass flow rate, the BSFC is reduced. Therefore, the addition of the engines' and ORC systems' power output would naturally decrease the BSFC value. Figure 13 shows the difference in BSFC.

Overall, the modelling and optimization of the ORC system managed to indicate a total average BSFC reduction of 4.78%, as explained in Table 4. The most substantial percentage value (5.6%) occurred at maximum engine operating speed point (X6). After the introduction of the ORC system, it can be observed that X2 is no longer the lowest BSFC point. A 4.36% BSFC reduction at X3 was enough to shift the ideal thermal efficiency benchmark.

**Figure 13.** Reduction of BSFC in ORC system.



The reduction in BSFC during system operation at point X2 presents a value of 3.83%, which is also the minimum reduction amount for the given engine. This means that the power unit at point X2 is already working at its peak thermal efficiency of approximately 43%, as earlier observed in Figure 11. Therefore, any further increments of this peak value are remarkably challenging to accomplish due to the reduced exhaust mass flow rate and temperature. Oppositely, the reason the BSFC reduction percentage is the greatest at X6 (5.6%) owns mostly to the inflated exhaust gas mass flow rate, which implements the largest impact as proved previously. Figure 14 illustrates the net power output to the BSFC reduction percentages.

**Figure 14.** Variation of BSFC in accordance with ORC net power output.

The thermal efficiency of the ORC system model is plotted in Figure 15, wherein the thermal efficiency is calculated to be between 4% and 8%. It is not surprising that the thermal efficiency of the ORC system (ïtherm) is at a very low mark considering that the pump and turbine were never designed to work in accordance with the specific engine exhaust outlets.

**Figure 15.** Thermal efficiency ORC system.

In fact, the thermal efficiency values would decrease even further. For example, if the system had not been configured as adiabatic, heat losses through the surroundings would, as a supplementary, encourage an even less efficient ORC operation. The low thermal efficiency of the ORC system contributed to the selection of the organic fluid, R245fa or any other for that matter. The decreased temperature input required for operation provokes the decrease in thermal efficiency. In the thermal efficiency graph, shown in Figure 15, there is an apparent inflation after the 2000 rpm mark. Nevertheless, despite the improved efficiency at high engine speeds, during high engine speeds and loads the ORC system is not able to take advantage of the additional and excessive amount of waste heat.

Engine efficiency can be seen to increase by a maximum of 5.69% at 1500 rpm. This is one of the speeds of interest for a heavy-duty engine, with the speed range of 1200–1500 rpm being, generally, the region of interest. The achieved efficiency compared well with engine data taken from the Brunel

University London engine ORC test facility and reported by Alshammari et al [68], at least from the point of view of cycle efficiency (4.3%) against a value of 4.95 to 5.69% in the region of interest in Figure 15.

#### *3.3. Turbo-Compound System*

#### 3.3.1. Turbo-Compound System Speed Variation

Identical to the ORC systems' process, the T/C system is run through the six operating benchmarks (X1, X2, X3, X4, X5, and X6) which were defined during engine calibration. The power output of the T/C system model varies significantly with transitional turbine rotational speeds for a given operational occasion. It was revealed that the T/C system performed diversely for bottom, mid, and top range engine speeds. Figure 16 represents the ideal turbine speed configuration for each assessment point. During low fuel, mass flow rate, and load (point X1), the turbine was more productive by operating at high speeds of over 100,000 rpm. During testing at the highest thermal efficiency (point X2) and part load (point X3), a mediocre turbine speed was ideal, showing peak performance at 60,000 rpm. However, during top end runs (points X5 and X6) it is observable that high turbine speeds achieve the best power outputs. In particular, with a turbine speed setting of 100,000 rpm during max power output and a max engine speed (points X5, X6), the T/C system generates 34.66 kW and 35.68 kW respectively. The full extent of the variable speed turbine results is listed in Table 5. The peak turbine power outputs for each point are highlighted in red.

**Figure 16.** T/C system turbine power output in relation to turbine speed variation.


**Table 5.** Turbine analytical power output to speed variation.

#### 3.3.2. Turbo-Compound System BSFC Reduction

Similar to the ORC system, the comparison purposes require the use of optimum power values, despite the fact that the yields are obtained using different turbine speeds. In addition, the BSFC difference between the baseline and turbo-compound engine is calculated. There is an important difference in the analysis of the T/C and ORC systems. The ORC system is modelled on a separate template whereas the T/C system is placed and assessed directly on the stock engine model as an integrated unit. As explained during the simulation section, due to the incorporation of the secondary turbine, the exhaust backpressure inflated, causing additional engine pumping losses. Figure 17 shows the escalation of exhaust backpressure after the introduction of the T/C system.

**Figure 17.** T/C system backpressure increment.

As expected by the increase in exhaust backpressure, the pumping losses are reflected by a proportional increase in BSFC. In fact, if the generated power output from the T/C system is not taken into account for the calculation of BSFC, the escalation in backpressure alone is enough to deteriorate the BSFC by a total average of almost 5.5% at low engine speeds. However, it is worth mentioning that this deterioration is, mostly, only distinguishable during low engine speeds, specifically as listed in Table 6. The results show an average of 3% increase in BFSC before applying the positive aspects of the T/C system in the standard engine model and, thus, BSFC is decreased, due to the harvested power. It is observable that this increase in BSFC is significant from point X1 to point X3 at an average of 5.5%. On the flipside, the BSFC values for points X4, X5, and X6 remain identical, with an average difference of less than 0.5%. This relation suggests that the higher exhaust backpressure only affects the engine

thermal efficiency at low engine speeds. With an overall consideration, the effects of the T/C system on the engine are represented in Figure 18.


**Table 6.** T/C system analytical BSFC increments.

**Figure 18.** T/C system BSFC reduction.

It is notable that, in excess of 2000 rpm, the effect of exhaust backpressure on BSFC was essentially annihilated. After that moment, the T/C system had the opportunity to commence with the positive characteristics of its nature. The beneficial aspects are also evident when plotting the secondary turbines' power output with the BSFC reduction percentage, plotted on Figure 19. Notice that the BSFC reduction axis has a minimum value of zero. This is done to highlight the major benefits of the T/C system during top end operation, which touch a maximum value of 7%. Table 7 includes the specific reduction values at each operational point.

Identically to the ORC system, the point of maximum engine thermal efficiency, the lowest BSFC point, is shifted to the right-hand side, namely, from the minimum point, X2 (214.773 g/kW·h), to previous point of maximum torque output, X4 (213.884 g/kW·h). Major BSFC reductions are shown after point X4, meaning that if the given engine spends most of its operating time at high speeds the T/C system would have a theoretical reduction in BSFC at an approximate average of 6.5%.


**Table 7.** T/C analytical BSFC reduction.

**Figure 19.** T/C power output vs. BSFC reduction.

#### 3.3.3. Turbo-Compound System Efficiency

The thermal efficiency of the T/C system follows a similar fashion to the ORC system. Therefore, it is calculated by the amount of power gained over the exhaust energy input in the form of surplus heat. However, unlike the work input necessity for the ORC systems' pump, there was no work input required in the T/C system. As far as the exhaust energy is concerned, it was again calculated by the product of exhaust gas specific heat, mass flow rate, and temperature difference between the secondary turbines' inlet and outlet. Figure 20 illustrates the T/C systems' thermal efficiency plot.

**Figure 20.** T/C system thermal efficiency.

A maximum of 10.31% of thermal efficiency may not seem significant, but it is heat that would otherwise be wasted. In addition, unlike the ORC system, of which the efficiency remains relatively stationary after the 2000-rpm mark, the T/C systems shows a continuously progressive increase.

#### *3.4. ORC T*/*C System vs. ORC System Comparison*

By assessing the two WHR methods individually, a solid base of strength and weakness points was set. It was noticeable that at point X6, which is the point for maximum engine speed, both WHR methods were at their peak BSFC reduction percentage. In fact, both produced maximum power output and reached maximum thermal efficiencies at point X6. This is due to the benchmarks' higher exhaust mass flow rate. On the other hand, the increment of exhaust mass flow rate subsequently reduced the exhaust gas temperature. As engine speed was increasing for a given engine load, fuel mass flow rate was also raised. At the same time, the compressor was forcing additional air into the cylinders. As a result, the combustion process improved due to the higher oxygen content within the combustion chamber, featured by the increased air mass flow rate and air to fuel ratio. The improved combustion converts the chemical energy supplied by the diesel to effective mechanical work more effectively, instead of wasting it as energy in the form of heat. However, the reduction in BSFC was naturally greater at higher BSFC values because they are susceptible to permit more room for improvement. For the same reasons, it was observed that for both the T/C and ORC systems, the lowest BSFC point is conveyed after the WHR implementations.

Despite the excessive turbulent flow that undermines the already unstable exhaust gas pulse at the measuring point (after the turbine), the temperature of the exhaust gas demonstrated an adequate stability, which further benefits the WHR methods. For a clearer contrast resolution, the BSFC values from points X1 to X6 are plotted on the same graph for both the ORC and T/C systems while using the original plot as a gauge, as shown in Figure 21.

**Figure 21.** Comparison of total BSFC between ORC and T/C in relation to engines' speed.

The results indicated that the ORC system was more favorable during low engine speed operation. This is because the secondary turbine in the T/C system cannot produce enough mechanical power to compensate for the additional exhaust backpressure during low-end operation. However, it was a different story when the engine was running at high speed. At points X4, X5, and X6, the exhaust backpressure did not have a major impact on the BSFC value. In addition, the turbine thrived, due to the increased exhaust mass flow rate, and the power gain was in excess of 35 kW at the finishing point. That was not only enough to compensate for the low engine speed losses, but also to further decrease the overall BSFC by an average of 2.3%. The comparison of BSFC between the two WHR systems as well as the baseline engine is also visualized in Figure 22. As far as the physical mass and volume properties of the two WHR systems are concerned, both are at approximately the same level. A mechanical T/C system would naturally result in more weight, due to the incorporated gearbox unit.

**Figure 22.** Comparison of total BSFC between ORC and T/C in relation to six individual points.

#### **4. Conclusions**

A comparison between two WHR systems using 1-D engine code simulation in GT-Power was conducted. The main aspect considered was their ability to reduce BSFC of an experimental engine model. To recapitulate, the effects of each systems' integration are listed below, as follows:


The ORC system provided a more consistent WHR method, with progressive improvements in fuel consumption across the engines' speed range. However, the T/C system presents incomparable contributions in fuel economy during high-speed engine operation.

**Author Contributions:** A.M.A. and A.P. have written the paper context and performed the simulation works alongside the results presentation. V.E., A.S.-Z., and A.H. have carried out the design of experiment in the simulations.

**Funding:** This research received no external funding.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


© 2019 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 (http://creativecommons.org/licenses/by/4.0/).

## *Article* **Configuration Optimization and Performance Comparison of STHX-DDB and STHX-SB by A Multi-Objective Genetic Algorithm**

#### **Zhe Xu 1,\*, Yingqing Guo 1, Haotian Mao <sup>1</sup> and Fuqiang Yang <sup>2</sup>**


Received: 16 March 2019; Accepted: 7 May 2019; Published: 11 May 2019

**Abstract:** Based on the thermohydraulic calculation model verified in this study and Non-dominated Sorted Genetic Algorithm-II (NSGA-II), a multi-objective configuration optimization method is proposed, and the performances of shell-and-tube heat exchanger with disc-and-doughnut baffles (STHX-DDB) and shell-and-tube heat exchanger with segmental baffles (STHX-SB) are compared after optimization. The results show that, except in the high range of heat transfer capacity of 16.5–17 kW, the thermohydraulic performance of STHX-DDB is better. Tube bundle diameter, inside tube bundle diameter, number of baffles of STHX-DDB and tube bundle diameter, baffle cut, number of baffles of STHX-SB are chosen as design parameters, and heat transfer capacity maximization and shell-side pressure drop minimization are considered as common optimization objectives. Three optimal configurations are obtained for STHX-DDB and another three are obtained for STHX-SB. The optimal results show that all the six selected optimal configurations are better than the original configurations. For STHX-DDB and STHX-SB, compared with the original configurations, the heat transfer capacity of optimal configurations increases by 6.26% on average and 5.16%, respectively, while the shell-side pressure drop decreases by 44.33% and 19.16% on average, respectively. It indicates that the optimization method is valid and feasible and can provide a significant reference for shell-and-tube heat exchanger design.

**Keywords:** shell-and-tube heat exchanger; disc-and-doughnut baffles; segmental baffles; multi-objective configuration optimization; genetic algorithm

#### **1. Introduction**

Shell-and-tube heat exchanger (STHX) is one kind of mechanical device that can exchange heat through a tube wall between two fluids of different temperatures, which are widely used in oil cooling in aero-engines, heating, chemical and food industries, and so on [1]. In the process of heat exchanger design, different configuration parameters can lead to different performances. In order to obtain the better performances, heat exchanger optimization is increasing significantly [2].

In the pursuit of optimized designs, configuration optimization as well as optimization methods application have been analyzed from different points of view. The genetic simulated annealing algorithm [3], biogeography-based algorithm [4], firefly algorithm [5], Jaya algorithm [6] and differential evolution algorithm [7] were used to minimize the cost of the heat exchanger. However, in the actual design process of the heat exchanger, total cost minimization is not the only objective. Many other objectives need to be considered, just like heat transfer capacity maximization, pressure drop

minimization, and so on. Thus, the results of research that can consider more than one objective simultaneously could be more meaningful.

In the research of Sanaye et al. [8], the Bell-Delaware procedure and ε-NTU method were used to calculate the shell side heat transfer coefficient and pressure drop of shell-and-tube heat exchanger with segmental baffles (STHX-SB), and the fast and elitist non-dominated sorted genetic algorithm was utilized to solve the optimization problem, which was to maximize the heat transfer coefficient and minimize pressure drop simultaneously, and obtain the Pareto optimal solutions. A method of a multi-objective genetic algorithm combined with numerical simulation was used by Wen et al. [9] and Wang et al. [10,11] to optimize the configuration of STHXs with helical baffles, and the Pareto optimal points of multi-objective optimization were obtained. In addition, this method was also used in the configuration optimization of a spiral-wound heat exchanger by Wang et al. [12]. Some novel theories were also used to design, calculate or optimize heat exchangers. Both Amidpour et al. [13] and Mirzaei et al. [14] applied the constructal theory combined with a genetic algorithm for STHX optimization. An optimization strategy that combined the entransy theory, numerical simulation and a genetic algorithm was proposed by Gu et al. [15] to minimize the entransy dissipation thermal resistance of STHX with helical baffles. Similarly, Chahartaghi et al. [16] used the entransy dissipation theory and Non-dominated Sorted Genetic Algorithm-II (NSGA-II) to minimize the entransy dissipation numbers separately caused by thermal conduction and fluid friction for STHX-SB. In addition, some other novel optimization algorithms were also proposed to optimize STHXs, just like the modified version of teaching-learning-based optimization [17], particle swarm optimization [18], Taguchi method [19], bat algorithm [20], and multi-objective heat transfer search algorithm [21]. In the investigation of Rao et al. [22], an artificial immune system is used with a generalized regression neural network to solve the optimization problem of heat transfer and overall pressure drop for a STHX-SB. Saldanha et al. [23] used Predator–Prey, Multi-objective Particle Swarm Optimization, and NSGA-II evolutionary algorithms to tackle the optimization problem of STHX and chose the best evolutionary algorithms by the Preference Ranking Organization Method for the Enrichment Evaluations decision-making method.

This research had successfully used theoretical approaches or numerical simulation to get the relationships between some configuration parameters and heat exchanger performances, and the corresponding configuration parameters were optimized based on some multi-objective optimization algorithms. Unfortunately, there was only one kind of heat exchanger optimized in each research, and the heat exchanger mentioned was not compared with some different kinds of heat exchangers which are commonly used. Thus, it is worth mentioning that Wang et al. [24] compared heat transfer and pressure drop performances among the STHXs separately with segmental baffles, continuous helical baffles and staggered baffles, and used the multi-objective optimization using a genetic algorithm to obtain the optimal Pareto front of STHX with staggered baffles. However, the comparison was only between the heat exchangers before optimization. The comparison of heat exchangers after optimization under the same condition could reveal their respective advantages and disadvantages after their corresponding configuration parameters were optimal.

In this paper, the configuration parameters of tube bundle diameter, inside tube bundle diameter, number of baffles of shell-and-tube heat exchanger with disc-and-doughnut baffles (STHX-DDB) and the configuration parameters of tube bundle diameter, baffle cut, number of baffles of STHX-SB are analyzed for heat transfer capacity and shell-side pressure drop based on thermohydraulic calculation of the Slipcevie method and Bell-Delaware method [25–29]. In addition, a multi-objective optimization method based on NSGA-II is utilized to optimize these two kinds of heat exchangers and compare the performances of them after optimization.

#### **2. Physical Model and Thermohydraulic Calculation Model**

#### *2.1. Physical Model and Configuration Parameters*

Both STHX-DDB and STHX-SB primarily consist of heat transfer tubes, baffles, end plates, spacer tubes and shell. The spacer tubes are strung on four heat transfer tubes to locate the positions of baffles. All the baffles are uniformly distributed between two end plates. The differences between disc-and-doughnut baffles and segmental baffles result in the different fluid flow styles of shell-side between STHX-DDB and STHX-SB, as shown in Figure 1. The main configuration parameters of STHX-DDB and STHX-SB are shown in Table 1. The material of STHXs is an aluminium alloy.

**Figure 1.** Differences between shell-and-tube heat exchanger with disc-and-doughnut baffles (STHX-DDB) and shell-and-tube heat exchanger with segmental baffles (STHX-SB): (**a**) STHX-DDB; (**b**) STHX-SB.


**Table 1.** Configuration parameters of STHX-DDB and STHX-SB.

To simplify calculations, assumptions are demonstrated as follows: (1) the leakages between baffles and tubes and those between baffles and shell are neglected; (2) the heat transfer between the fluid inside heat exchanger and the environment outside heat exchanger is neglected; (3) both the fluid of shell-side and tube-side are Newtonian and incompressible fluid with constant properties; and (4) the heat exchangers are new and there is no fouling.

#### *2.2. Thermohydraulic Calculation Model*

Thermohydraulic calculation of STHXs includes heat transfer capacity calculation and shell-side pressure drop calculation. Basic formulas of heat transfer capacity calculation and shell-side pressure drop calculation are given in Appendices A and B, respectively. The main steps of this thermohydraulic calculation model are shown as follows:


The flowchart of thermohydraulic calculation model is shown in Figure 2.

**Figure 2.** Flowchart of the thermohydraulic calculation model.

#### **3. Optimization Method**

#### *3.1. Multi-Objective Genetic Algorithm*

It is quite difficult to get the absolute optimal solution of a multi-objective problem because of the conflict of objectives. However, a Pareto optimal front, which consists of a series of solutions, can be obtained through NSGA-II [30]. NSGA-II incorporates elitism [31], and any solution of the Pareto optimal front is generated at the price of other lower priority performances. The ranking scheme of this method adopts the non-dominated sorting method, which is faster than the traditional Pareto ranking method. The constraint handling method also uses a non-dominance principle as the objective, thus penalty functions and Lagrange multipliers are not needed, which guarantees that the feasible solutions are always ranked higher than the unfeasible solutions [32]. A genetic algorithm is a stochastic global searching probabilistic method, which emulates the genetic mechanism of Darwinian evolution and the process of natural selection. This method contains the main steps of selection, crossover and mutation and has the characteristics of self-adaptive, parallel, random [33].

The main steps of the optimization method are shown as follows:


Details of crossover and mutation are referred in [34]. The flowchart of the optimization method is shown in Figure 3.

**Figure 3.** Flowchart of the optimization method.

#### *3.2. Objective Functions and Design Parameters*

In this case, in consideration of the thermohydraulic performances of STHX-DDB and STHX-SB, the heat transfer capacity *Q* and shell-side pressure drop Δ*P* are selected as objective functions. For STHX-DDB, tube bundle diameter *Do*, inside tube bundle diameter *Di* and number of baffles *Nb* are regarded as design parameters, and the multi-objective optimization problem can be formulated as Equation (1):

$$\begin{cases} \text{Min} & -Q, \Delta P \\ \text{S.t.} & 80 \text{ mm} \le D\_{\vartheta} \le 113 \text{ mm} \\ & 10 \text{ mm} \le D\_{i} \le 70 \text{ mm} \\ & N\_{b} = 2,3,4,5,6 \end{cases} \tag{1}$$

For STHX-SB, tube bundle diameter *Do*, baffle cut *P* and number of baffles *Nb* are regarded as design parameters, and the multi-objective optimization problem can be formulated as Equation (2):

$$\begin{cases} \text{Min} & -Q, \Delta P \\ \text{S.t.} & 80 \text{ mm} \le D\_o \le 113 \text{ mm} \\ & 0.15 \le P \le 0.45 \\ & N\_b = 2, 3, 4, 5, 6 \end{cases} \tag{2}$$

Except for the design parameters, the other configuration parameters of STHX-DDB and STHX-SB are shown in Table 1. The working conditions are shown in Table 2.


**Table 2.** Experiment conditions for STHX-DDB and STHX-SB.

#### **4. Results and Discussion**

#### *4.1. Model Validation*

In order to validate the accuracy of thermohydraulic calculation model, calculation results are compared with experimental data. In this case, the tube-side fluids of these two heat exchangers, which are heated, are both fuel, and the shell-side fluids, which are cooled, are both oil. The experiment conditions and the testing equipment for STHX-DDB and STHX-SB are given in Table 2 and Figure 4, respectively, and the comparison results are shown in Figure 5.

25

25

25

**Figure 4.** Testing equipment: (**a**) testing equipment of the oil side; (**b**) one part of testing equipment of the fuel side; (**c**) the other part of testing equipment of the fuel side.

**Figure 5.** Comparison between calculation results and experimental data: (**a**) heat transfer capacity; (**b**) shell-side pressure drop.

The results illustrate that both the heat transfer capacity and the shell-side pressure drop obtained from calculation model are in agreement with the experimental results. The errors of the heat transfer capacity of STHX-DDB and STHX-SB are 1.5–3.0% and 1.1–1.7% with the average deviations of 2.4% and 1.4%, respectively, and those of shell-side pressure drop are 10.8–14.4% and 10.0–13.0% with the average deviations of 13.0% and 11.6%, respectively. Hence, it can be concluded that the thermohydraulic calculation model is reliable and the errors are acceptable. The differences may be caused by simplification of physical model, deviation of formula itself, and unavoidable measurement error in experiment.

#### *4.2. Design Parameters E*ff*ects of STHX-DDB*

#### 4.2.1. Effects of Tube Bundle Diameter

The effects of tube bundle diameter *Do* on heat transfer capacity *Q* and shell-side pressure drop Δ*P* of STHX-DDB are represented in Figure 6 while inside tube bundle diameter and number of baffles are 20 mm and 3, respectively. As depicted, while tube bundle diameter increases from 80 mm to 113 mm, heat transfer capacity increases by 11.47–16.19 kW, and shell-side pressure drop drops slowly from 14.76 kPa first, reaching its lowest point of 14.27 kPa when tube bundle diameter is 102 mm, and then rises more and more dramatically to 19.02 kPa. As shown in Figure 7, as tube bundle diameter increases, the number of tubes increases. Thus, heat transfer area increases, which enhances heat transfer capacity. In addition, the growth of tube bundle diameter means that the flow area between baffles increases and the flow area between disc baffles and shell decreases, which results in the decrease of the pressure drop between baffles and the increase of the pressure drop between disc baffles and shell. Moreover, the pressure drop between disc baffles and shell grows increasingly significantly after tube bundle diameter reaches 102 mm. Thus, shell-side pressure drop falls increasingly slowly first and then grows increasingly sharply after tube bundle diameter reaches 102 mm.

**Figure 6.** Heat transfer capacity and shell-side pressure drop separately versus tube bundle diameter.

**Figure 7.** Tube number and heat transfer area separately versus tube bundle diameter.

#### 4.2.2. Effects of Inside Tube Bundle Diameter

The effects of inside tube bundle diameter *Di* on heat transfer capacity *Q* and shell-side pressure drop Δ*P* of STHX-DDB are represented in Figure 8 while tube bundle diameter and number of baffles are 100 mm and 3, respectively. As depicted, while inside tube bundle diameter increases from 10 mm to 70 mm, heat transfer capacity decreases by 14.82–9.80 kW, and shell-side pressure drop drops significantly from 154.42 kPa first, reaching the point of 5.96 kPa when inside tube bundle diameter is 32 mm, and then falls slightly to 4.05 kPa. As shown in Figure 9, the increase of inside tube bundle diameter results in the decrease of number of tubes and the growth of flow area of the central hole of doughnut baffles. Thus, both heat transfer capacity and shell-side pressure drop decrease. As the impact on shell-side pressure drop is more serious before the point of 5.96 kPa than after, the figure decreases more seriously before this point than after.

#### 4.2.3. Effects of Number of Baffles

The effects of number of baffles *Nb* on heat transfer capacity *Q* and shell-side pressure drop Δ*P* of STHX-DDB are represented in Figure 10 while tube bundle diameter and inside tube bundle diameter are 100 mm and 20 mm, respectively. As depicted, while number of baffles increases from 2 to 6, heat transfer capacity and shell-side pressure drop increase by 14.21–14.95 kW and 8.63–25.60 kPa, respectively. As number of baffles increases, the distance between baffles decreases, which leads to the increase of fluid velocity between baffles. Thus, both heat transfer capacity and shell-side pressure drop increase. In addition, as number of baffles grows, disc baffle number and doughnut baffle number increase alternately, which results in an unsmooth figure of shell-side pressure drop.

**Figure 8.** Heat transfer capacity and shell-side pressure drop separately versus inside tube bundle diameter.

**Figure 9.** Tube number and heat transfer area separately versus inside tube bundle diameter.

**Figure 10.** Heat transfer capacity and shell-side pressure drop separately versus number of baffles.

#### *4.3. Design Parameter E*ff*ects of STHX-SB*

#### 4.3.1. Effects of Tube Bundle Diameter

The effects of tube bundle diameter *Do* on heat transfer capacity *Q* and shell-side pressure drop Δ*P* of STHX-SB are represented in Figure 11 while baffle cut and number of baffles are 0.25 and 3, respectively. As depicted, while tube bundle diameter increases from 80 mm to 113 mm, heat transfer capacity and shell-side pressure drop increase by 10.16–16.60 kW and 6.59–57.07 kPa, respectively. As shown in Figure 12, as tube bundle diameter grows, the number of tubes increases, which improves

heat transfer capacity. At the same time, both the flow areas of cross flow region and window region decrease, which leads to a higher shell-side pressure drop.

**Figure 11.** Heat transfer capacity and shell-side pressure drop separately versus tube bundle diameter.

**Figure 12.** Tube number and heat transfer area separately versus tube bundle diameter.

#### 4.3.2. Effects of Baffle Cut

The effects of baffle cut *P* on heat transfer capacity *Q* and shell-side pressure drop Δ*P* of STHX-SB are represented in Figure 13 while tube bundle diameter and number of baffles are 100 mm and 3, respectively. As depicted, while baffle cut increases from 0.15 to 0.45, heat transfer capacity and shell-side pressure drop decrease by 14.50–12.44 kW and 16.74–9.99 kPa, respectively. The increase of baffle cut means that cross flow decreases and parallel flow increases. As cross flow has a better heat transfer effect and a bigger pressure drop than parallel flow, both heat transfer capacity and shell-side pressure drop decrease.

**Figure 13.** Heat transfer capacity and shell-side pressure drop separately versus baffle cut.

#### 4.3.3. Number of Baffles

The effects of number of baffles *Nb* on heat transfer capacity *Q* and shell-side pressure drop Δ*P* of STHX-SB are represented in Figure 14 while tube bundle diameter and baffle cut are 100 mm and 0.25, respectively. As depicted, while number of baffles increases from 2 to 6, heat transfer capacity and shell-side pressure drop increase by 13.88–14.25 kW and 9.17–43.22 kPa, respectively. As number of baffles increases, turbulence of fluid is strengthened. Thus, both heat transfer capacity and shell-side pressure drop increase.

**Figure 14.** Heat transfer capacity and shell-side pressure drop separately versus number of baffles.

#### *4.4. Multi-Objective Optimization Results*

In order to consider the comprehensive thermohydraulic performances of STHX-DDB and STHX-SB, the multi-objective configuration optimization using NSGA-II is conducted. The optimization goals are to maximize heat transfer capacity and to minimize shell-side pressure drop, and two objective functions (−*Q* and Δ*P*) are conflicting. Population size, maximum iteration number, analog binary cross distribution index and polynomial mutation distribution index, which are main settings of optimization, are 100, 500, 20 and 20, respectively. Figure 15 shows the Pareto optimal points which are obtained based on the objective functions and the design parameters described in Equations (1) and (2), and the Pareto optimal points in the heat transfer capacity ranges of 7.5–16.5 kW and 16.5–17 kW are shown in Figures 16 and 17, respectively.

**Figure 16.** Pareto optimal points for the heat transfer capacity range of 7.5–16.5 kW.

**Figure 17.** Pareto optimal points for the heat transfer capacity range of 16.5–17 kW.

As shown in Figure 15, for both heat exchangers, shell-side pressure drop increases with the growth of heat transfer capacity. Especially in the range of high heat transfer capacity, shell-side pressure drop increases significantly sharply. As depicted in Figure 16, after optimization, for the same heat transfer capacity, the shell-side pressure drop of STHX-DDB is lower than this value of STHX-SB in the heat transfer capacity range of 7.5–16.5 kW. Thus, STHX-DDB has a better thermohydraulic performance than STHX-SB in this range. As shown in Figure 17, after optimization, for the same heat transfer capacity, the shell-side pressure drop of STHX-DDB is above or almost equal to this value of STHX-SB in the heat transfer capacity range of 16.5–17 kW. Thus, in this high range of heat transfer capacity, the thermohydraulic performance of STHX-SB is better.

Three Pareto optimal solutions for STHX-DDB and three Pareto optimal solutions for STHX-SB are chosen and shown in Figure 16, Tables 3 and 4. For either STHX-DDB or STHX-SB, the selection method of three Pareto optimal solutions follows the following principles:



**Table 3.** Optimization thermohydraulic performances comparison of STHX-DDB.


**Table 4.** Optimization thermohydraulic performances comparison of STHX-SB.

Compared with the original configuration of STHX-DDB shown in Table 3, the heat transfer capacity values of optimal configuration A, B and C increase separately by 10.56%, 7.11% and 1.10%, and the shell-side pressure drop values of them decrease by 15.20%, 49.30% and 68.49%, respectively. It is clearly depicted that the heat transfer capacity of optimal configurations of STHX-DDB increases by 6.26% on average, while the shell-side pressure drop decreases by 44.33% on average. Compared with the original configuration of STHX-SB shown in Table 4, the heat transfer capacity values of optimal configuration D, E and F increase separately by 8.35%, 5.14% and 2.00%, and the shell-side pressure drop values of them decrease by 5.61%, 19.51% and 32.35%, respectively. It is clearly shown that the heat transfer capacity of optimal configurations of STHX-SB increases by 5.16% on average, while the shell-side pressure drop decreases by 19.16% on average. Based on the comparison and analysis of original and optimal configurations, it indicates that the optimization method is valid and feasible, and the comprehensive thermohydraulic performances of STHX-DDB and STHX-SB can be enhanced by the optimization method.

#### **5. Conclusions**

In this paper, a multi-objective configuration optimization method for STHX-DDB and STHX-SB based on NSGA-II is proposed and used. Based on the thermohydraulic calculation model verified in this paper, the effects of design parameters are analyzed. Then, the optimal design parameters (tube bundle diameter, inside tube bundle diameter, number of baffles) for STHX-DDB and the optimal design parameters (tube bundle diameter, baffle cut, number of baffles) for STHX-SB are found separately by a trade-off between the two common objectives of heat transfer capacity maximization and shell-side pressure drop minimization.

The distribution of Pareto optimal points reveals that the two objectives are conflicted, and, for both heat exchangers, shell-side pressure drop increases significantly sharply in the range of high heat transfer capacity, which means that high heat transfer capacity will cost a high shell-side pressure drop. In addition, after optimization, except in the high range of heat transfer capacity of 16.5–17 kW, the shell-side pressure drop of STHX-DDB is lower than this value of STHX-SB when the heat transfer capacity values of them are the same. In other words, under this condition, the thermohydraulic performance of STHX-DDB is better.

Finally, three optimal solutions for STHX-DDB and three optimal solutions for STHX-SB are selected and compared with the original configurations. The results illustrate that all the selected optimal configurations are better than the original configurations. For STHX-DDB, the heat transfer capacity of optimal configurations increases by 6.26% on average compared with the original configuration, while the shell-side pressure drop decreases by 44.33% on average. For STHX-SB, the heat transfer capacity of optimal configurations increases by 5.16% on average compared with the original configuration, while the shell-side pressure drop decreases by 19.16% on average. It indicates that the optimization method is feasible and valid and can provide a significant reference for STHX design.

**Author Contributions:** Conceptualization, Z.X.; Investigation, Z.X. and H.M.; Methodology, Z.X.; Project administration, Y.G.; Resources, Z.X.; Software, Z.X.; Supervision, Y.G.; Validation, Z.X.; Visualization, Z.X.; Writing—original draft, Z.X.; Writing—review and editing, Z.X., Y.G. and F.Y.

**Funding:** This research received no external funding.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **Nomenclature**




#### **Appendix A**

Basic formulas of heat transfer capacity calculation. The basic formula of heat transfer capacity calculation of heat exchanger is given by Equation (A1) [29]:

$$
\Omega = \eta \mathcal{W}\_{\text{min}}(t\_1 - t\_2). \tag{A1}
$$

As the number of shell passes and tube passes are 1 and 2, respectively, η is given by Equations (A2)–(A6) [29,35]:

$$\eta = \frac{2}{1 + \mathbb{C}^\* + \sqrt{1 + \mathbb{C}^2}(1 + e^{-\text{NTU}\sqrt{1 + \mathbb{C}^2}})/(1 - e^{-\text{NTU}\sqrt{1 + \mathbb{C}^2}})},\tag{A2}$$

$$\mathbf{C}^\* = \mathcal{W}\_{\text{min}} / \mathcal{W}\_{\text{max}} \,\tag{A3}$$

NTU <sup>=</sup> *KA W*min , (A4)

$$\frac{1}{K} = \frac{1}{\alpha\_i} (\frac{d\_o}{d\_i}) + r\_{s,i}(\frac{d\_o}{d\_i}) + \frac{\delta\_{w}}{\lambda\_{w}}(\frac{d\_o}{d\_{\text{m}}}) + r\_{s,o} + \frac{1}{\alpha\_o} \tag{A5}$$

$$d\_m = \frac{d\_o - d\_i}{\ln(d\_o/d\_i)}.\tag{A6}$$

In general, the thermal-conduction resistance of a metal tube wall is much smaller than the thermal resistance of fluids' heat convection, and, for the new heat exchanger, the fouling resistance can be neglected. In addition, the tube wall thickness used in this case is only 0.305 mm, which is too small to be ignored. Thus, Equation (A5) can be simplified to Equation (A7) [29]:

$$K = \frac{\alpha\_o \alpha\_i}{\alpha\_o + \alpha\_i}.\tag{A7}$$

For STHXs, α*<sup>i</sup>* is given by Equations (A8)–(A10) [26]:

$$
\mu\_i = 0.027 \frac{\lambda\_i}{d\_i} Re\_i^{0.8} \Pr\_i^{1/3} \left(\frac{\mu\_i}{\mu\_{inv}}\right)^{0.14} \tag{A8}
$$

(*Rei* > 10,000, 0.7 < *Pri* < 16,700, *<sup>L</sup> di* ≥ 60),

$$\alpha\_{i} = 0.027(1 - \frac{6 \times 10^{5}}{Re\_{i}^{1.8}}) \frac{\lambda\_{i}}{d\_{i}} Re\_{i}^{0.8} Pr^{1/3} \left(\frac{\mu\_{i}}{\mu\_{inv}}\right)^{0.14} \tag{A9}$$

(2300 <sup>≤</sup> *Rei* <sup>≤</sup> 10,000, 0.7 <sup>&</sup>lt; *Pri* <sup>&</sup>lt; 16,700, *<sup>L</sup> di* ≥ 60),

$$\alpha\_{i} = 1.86 \frac{\lambda\_{i}}{d\_{i}} (\text{Re}\_{i} \Pr\_{i} \frac{d\_{i}}{L})^{1/3} \left(\frac{\mu\_{i}}{\mu\_{iw}}\right)^{0.14} \tag{A10}$$

(*Rei* <2300, 0.6 < *Pri* < 6700, *Rei* Pr*<sup>i</sup> <sup>L</sup> di* ≥ 100).

In Equations (A8)–(A10), μ*iw* is the viscosity of tube-side fluid under the temperature of tube wall. Because of the different shell-side structures between STHX-DDB and STHX-SB, α*<sup>o</sup>* calculation formulas are different. For STHX-DDB, the Slipcevie method is used, and α*<sup>o</sup>* is given by Equations (A11) and (A12) [26]:

$$\alpha\_{\ o} = \sum\_{j=1}^{n} \frac{A(j)}{A} \alpha\_{\ o}(j),\tag{A11}$$

$$\alpha\_{\vartheta}(j) = 0.33 \frac{\lambda\_o}{d\_o} \text{Re}\_o^{0.6}(j) \stackrel{0.33}{\text{Pr}}\_o \left(\frac{t\_o}{t\_w}\right)^{0.14}.\tag{A12}$$

In Equations (A11) and (A12), *A*(*j*), α*o*(*j*) and *Reo*(*j*) represent heat transfer area, heat transfer coefficient and Reynolds number of the fluid over the tube row of *j*, respectively.

For STHX-SB, the Bell–Delaware method is used and α*<sup>o</sup>* is given by Equation (A13) [26–29]:

$$
\alpha\_{\vartheta} = \alpha\_{\mathrm{id}} \|\_{\mathrm{c}} \|\_{\mathrm{l}} \|\_{\mathrm{b}} \|\_{\mathrm{s}} \mathrm{l}\_{\mathrm{r}}.\tag{A13}
$$

In Equation (A13), α*id* is the ideal heat transfer coefficient for pure cross flow over a tube bundle. *Jc*, *Jl*, *Jb*, *Js* and *Jr* represent correction factors for baffle cut, baffle leakage, bundle bypassing, variable baffle spacing at the inlet and outlet of shell, and adverse temperature gradient in laminar flow, respectively. *Re* and *Pr* used in the calculation above can be calculated by Equations (A14)–(A16) [29]:

$$\mathcal{Re} = \frac{\rho u d\_c}{\mu},\tag{A14}$$

$$\text{Pr} = \frac{\mu \mathcal{C}\_p}{\lambda},\tag{A15}$$

$$
\mu = \frac{q\_m}{\rho S}.\tag{A16}
$$

#### **Appendix B**

Basic formulas of shell-side pressure drop calculation. Because of the different shell-side structures between STHX-DDB and STHX-SB, formulas of Δ*P* calculation are different.

For STHX-DDB, Δ*P* is given by Equations (A17)–(A24) [36]:

$$
\Delta P = \Delta P\_N + \Delta P\_B + \Delta P\_{s1} + \Delta P\_{s2\prime} \tag{A17}
$$

*Energies* **2019**, *12*, 1794

$$
\Delta P\_N = 1.5 \frac{\rho\_0 u\_N^2}{2},
\tag{A18}
$$

$$
\Delta P\_B = \varepsilon\_B (N\_B + 1) \frac{\rho\_\vartheta u\_B^2}{2},
\tag{A19}
$$

$$
\Delta P\_{s1} = N\_{B1} \varepsilon\_{s1} \frac{\rho\_o u\_{s1}^2}{2},
\tag{A20}
$$

$$
\Delta P\_{s2} = N\_{B2} \varepsilon\_{s2} \frac{\rho\_o u\_{s2}^2}{2},
\tag{A21}
$$

$$\varepsilon\_B = \begin{cases} \frac{\frac{60}{L - d\_B}}{\frac{L - d\_B}{d\_B} \, R c\_B} & \left( R c\_B < \frac{42.3 d\_B}{L - d\_o} \right) \\\frac{3}{\left( \frac{L - d\_B}{d\_o} \, R c\_B \right)^{0.2}} & \left( R c\_B \ge \frac{42.3 d\_B}{L - d\_o} \right) \end{cases} \tag{A.22}$$

$$
\kappa\_{s1} = 2.2 + 286 R \varepsilon\_{s1}^{-0.845},
\tag{A23}
$$

$$
\varepsilon\_{\circ2} = 2.2 + 286 \text{Re}\_{\circ2}^{-0.845}. \tag{A24}
$$

For STHX-SB, the Bell-Delaware method is used and Δ*P* is given by Equations (A25)–(A28) [26–28].

$$
\Delta P = 2\Delta P\_{bk} R\_b (1 + \frac{N\_{\rm GW}}{N\_{\rm \varepsilon}}) R\_{\rm \varepsilon} + [(N\_B - 1)\Delta P\_{bk} R\_b + N\_B \Delta P\_{\rm wk}] R\_l + \Delta P\_{\rm N\prime} \tag{A25}
$$

$$
\Delta P\_{\rm lk} = 4f\_{\rm o} N\_{\rm c} \frac{\rho\_{\rm o} \mu\_{\rm lk}^2}{2} (\frac{\mu\_{\rm o}}{\mu\_{\rm ow}}) \stackrel{-0.14}{\ }\tag{A26}
$$

$$f\_o = \begin{cases} 47.1 Re\_{bk}^{-0.965} & (Re\_{bk} < 100) \\ 13.0 Re\_{bk}^{-0.685} & (100 \le Re\_{bk} \le 300) \\ 3.2 Re\_{bk}^{-0.44} & (300 < Re\_{bk} \le 1000) \\ 0.505 Re\_{bk}^{-0.176} & (Re\_{bk} > 1000) \end{cases} \tag{A27}$$

$$\Delta P\_{\rm uk} = \begin{cases} (2 + 0.6 \text{N}\_{\rm cw}) \frac{\rho\_0 u\_{\rm wk}^2}{2} & (\text{Re}\_{\rm wk} \ge 100) \\ 26 \mu\_o u\_{\rm wk} (\frac{\text{N}\_{\rm cw}}{p\_l - d\_o} + \frac{l\_b}{d\_{\rm wv}^2}) + 2 (\frac{\rho\_o u\_{\rm wk}^2}{2}) & (\text{Re}\_{\rm wk} < 100) \end{cases} . \tag{A28}$$

In Equations (A25)–(A28), *Rb*, *Rl*, *Rs* represent correction factors for bundle bypassing, baffle leakage, and variable baffle spacing at the inlet and outlet of shell, respectively. μ*ow* is the viscosity of shell-side fluid under the temperature of tube wall. *dew* is the equivalent diameter of tube bundle in the window region.

#### **References**


© 2019 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 (http://creativecommons.org/licenses/by/4.0/).

*Article*

## **Improving Performance of Simplified Computational Fluid Dynamics Models via Symmetric Successive Overrelaxation**

#### **Vojt ˇech Turek**

Sustainable Process Integration Laboratory–SPIL, NETME Centre, Faculty of Mechanical Engineering, Brno University of Technology–VUT Brno, Technická 2, 616 00 Brno, Czech Republic; turek@fme.vutbr.cz

Received: 23 May 2019; Accepted: 21 June 2019; Published: 25 June 2019

**Abstract:** The ability to model fluid flow and heat transfer in process equipment (e.g., shell-and-tube heat exchangers) is often critical. What is more, many different geometric variants may need to be evaluated during the design process. Although this can be done using detailed computational fluid dynamics (CFD) models, the time needed to evaluate a single variant can easily reach tens of hours on powerful computing hardware. Simplified CFD models providing solutions in much shorter time frames may, therefore, be employed instead. Still, even these models can prove to be too slow or not robust enough when used in optimization algorithms. Effort is thus devoted to further improving their performance by applying the symmetric successive overrelaxation (SSOR) preconditioning technique in which, in contrast to, e.g., incomplete lower–upper factorization (ILU), the respective preconditioning matrix can always be constructed. Because the efficacy of SSOR is influenced by the selection of forward and backward relaxation factors, whose direct calculation is prohibitively expensive, their combinations are experimentally investigated using several representative meshes. Performance is then compared in terms of the single-core computational time needed to reach a converged steady-state solution, and recommendations are made regarding relaxation factor combinations generally suitable for the discussed purpose. It is shown that SSOR can be used as a suitable fallback preconditioner for the fast-performing, but numerically sensitive, incomplete lower–upper factorization.

**Keywords:** computational fluid dynamics; symmetric successive overrelaxation; preconditioning; performance

#### **1. Introduction**

In engineering practice, it is often the case that process equipment is designed according to various rules of thumb. No optimization is generally done and, at best, a single computational fluid dynamics (CFD) simulation is carried out to verify that the design meets the key requirements of the future operator of the apparatus. This means that suboptimal designs or solutions, potentially leading to operating problems, are not uncommon.

One of the ways to remedy the situation is to use simplified CFD models. In spite of them not being as accurate as the standard CFD models, it has been shown [1] that they can provide useful quantitative information. What is more, these models feature significantly shorter computational times and their application in optimization algorithms is therefore much less cumbersome. To obtain solutions even faster, however, the numerical methods used to solve the underlying linear systems of equations can also be preconditioned. This means that instead of solving the original linear system

$$\mathbf{Ax} = \mathbf{b},\tag{1}$$

where **A** is the coefficient matrix, **x** the solution vector, and **b** the right-hand side vector, one considers the system

$$\mathbf{M}^{-1}\mathbf{A}\mathbf{x} = \mathbf{M}^{-1}\mathbf{b}.\tag{2}$$

here, **M** denotes the preconditioning matrix such that **M**−1**A** has a smaller condition number than **A** and, therefore, linear system (2) features better convergence. It should also be noted that **M**−1**A** often is not formed explicitly but, instead, **Mu** = **v** is solved for various auxiliary vectors **u** and **v** within the numerical solution method itself.

There are many different preconditioning techniques (i.e., ways to choose **M**) available, and the selection of the best one for a particular purpose depends mainly on the type of equation that is being solved and the employed ordering of the variables. However, the most commonly used techniques are likely various flavors of incomplete lower–upper factorization (ILU) [2] (these were numerically investigated by Chapman et al. [3]) and symmetric successive overrelaxation (SSOR) [4]. Although SSOR was originally intended for symmetric matrices, it was shown [5] to also work when the matrices are not symmetric. Assuming the splitting **A** = **D** + **L** + **U**, where **D** is the diagonal of **A** and **L** and **U** its strictly lower and upper triangular parts, respectively, the SSOR preconditioning technique is applied within the numerical solution method as two SOR sweeps using different values of ω. The iteration process for auxiliary vectors **u**, **v** (depend on the actual solution method used) can then be written as

$$\begin{aligned} \mathbf{u}^{(k+1/2)} &= (\mathbf{D} + \omega\_{\mathbb{F}} \mathbf{L})^{-1} [(1 - \omega\_{\mathbb{F}}) \mathbf{D} - \omega\_{\mathbb{F}} \mathbf{U}] \mathbf{u}^{(k)} + \omega\_{\mathbb{F}} (\mathbf{D} + \omega\_{\mathbb{F}} \mathbf{L})^{-1} \mathbf{v} \text{ (forward sweep)}\\ \mathbf{u}^{(k+1)} &= (\mathbf{D} + \omega\_{\mathbb{R}} \mathbf{U})^{-1} [(1 - \omega\_{\mathbb{R}}) \mathbf{D} - \omega\_{\mathbb{R}} \mathbf{L}] \mathbf{u}^{(k+1/2)} + \omega\_{\mathbb{R}} (\mathbf{D} + \omega\_{\mathbb{R}} \mathbf{U})^{-1} \mathbf{v} \text{ (backward scorep)}\end{aligned} \tag{3}$$

where ω<sup>F</sup> and ω<sup>R</sup> denote the corresponding forward and backward relaxation factors. In other words, direct application of the inverse of the preconditioning matrix, **M**<sup>−</sup>1, is replaced by preconditioned fixed point iteration.

The advantages of SSOR are evidenced by the existence of a multitude of papers discussing improved versions of this technique or its extensions to various specific applications. Bai [6] studied SSOR-like preconditioners for non-Hermitian positive definite linear systems, for which the respective matrix was either Hermitian-dominant or skew–Hermitian-dominant. The paper also discussed the results of numerical implementations, showing that Krylov subspace iteration methods, when accelerated using SSOR-like preconditioners, are efficient solvers for classes of non-Hermitian positive definite linear systems. A "shifted" version of SSOR for non-Hermitian positive definite linear systems with a dominant Hermitian part was proposed by Tan [7]. Zhang [8], on the other hand, introduced an SSOR-like preconditioner for saddle point problems with a dominant skew-Hermitian part. A class of hybrid preconditioning methods for accelerated solution of saddle point problems was discussed by Wang [9], while Chen et al. [10] proposed a version of SSOR suitable for preconditioning of large dense complex linear systems arising from three-dimensional electromagnetic scattering. Wu and Li [11] introduced a modified SSOR technique for the solution of Helmholtz equations.

Preconditioning can also be done block-wise. This was discussed, e.g., by Zhang and Cheng [12] in terms of large sparse saddle point problems, and by Huang and Lu [13], who focused on SSOR block preconditioners applied in image restoration. Because, in fact, preconditioning means obtaining an easily invertible approximation of the original matrix, one can also use SSOR for just this purpose as shown, for example, by Meng et al. [14] in the context of fast recovery of density 3D data from gravity data. Similarly, a massively-parallel GPU implementation of the conjugate gradient method, which uses the approximate inverse matrix derived from SSOR as the preconditioning matrix, was proposed by Helfenstein and Koko [15].

Performance comparisons of SSOR and other, simpler preconditioning techniques were presented, e.g., by Meyer [16], who focused on genomic evaluation and by Sanjuan et al. [17], who used SSOR to accelerate parallel wind field calculations. In the latter paper, the authors also evaluated a new, reordered sparse matrix storage format and showed that this format can markedly shorten computational time. This confirms the earlier findings of Duff and Meurant [18], who investigated the effect of ordering

on the convergence of the conjugate gradient method preconditioned, among others, using SSOR, or DeLong and Ortega [19], who focused on parallel implementations of SOR in terms of natural and multicolor orderings. Chen et al. [20], on the other hand, proposed a novel reordering technique for SSOR approximate inverse preconditioner, used together with a GPU-accelerated conjugate gradient solver, which should maximize the coalescing of global memory accesses.

Many SSOR-like, a priori preconditioned numerical solution methods using different splittings of the coefficient matrix have also been proposed for various types of problems. A three-parameter extension of the SSOR method intended for singular saddle point problems, which commonly arise, e.g., in fluid dynamics, was proposed by Li and Zhang [21]. A differently accelerated generalized three-parameter method for both singular and non-singular problems was introduced by Pan [22]. Similarly, a three-parameter unsymmetric SOR method for such saddle point problems was proposed, for example, by Liang and Zhang [23], who also discussed the choice of optimal values of the parameters. Many different SSOR-like methods are available for augmented systems, as well. Wang and Huang [24] introduced a four-parameter method, Louka and Missirlis [25] introduced a five-parameter extrapolated form of SSOR, and Najafi and Edalatpanah [26] introduced an improved version of the modified SSOR method for large sparse augmented systems, proposed earlier by Darvishi and Hessari [27]. Another improved SSOR method intended for the solution of augmented systems was proposed by Salkuyeh et al. [28]. In the case of complex systems, one can use, for instance, the accelerated method by Huang et al. [29], that is, an accelerated version of the method by Edalatpour et al. [30], in which the solution vector is split into two subvectors and different relaxation factors are used when solving for each of them. Likewise, one can employ the preconditioned variant of the generalized SSOR method by Hezari et al. [31] or the method by Salkuyeh et al. [32], which solves a real system obtained from the original, complex one. Block linear systems can be solved, e.g., using the block-preconditioned SSOR method by Pu and Wang [33].

The majority of the papers mentioned above discuss convergence (or at least semi-convergence) of the proposed methods, and many also include some information on the optimal selection of the relaxation factors. Kushida [34] focused on the estimation of convergence of the original SSOR preconditioner via a condition number, while general discussion related to SSOR-like methods for non-Hermitian positive definite linear systems was published in [35]. Augmented systems were addressed, for example, by Wang and Huang [36]. Similarly, there are papers focusing on convergence and optimal selection of the relaxation factors in the case of methods for block 2 × 2 linear systems [37], saddle point problems [38], parallel SSOR implementations [39], the Poisson equation [40], etc. In all these cases, however, convergence was investigated via spectral analysis, which is often prohibitively expensive [41]. The present paper, focusing on fast estimation of suitable SSOR relaxation factors in engineering practice, therefore, investigates the convergence experimentally using several different simplified 3D CFD flow models. The suitability of specific combinations of relaxation factors is assessed on the basis of mean computational times needed to reach converged steady-state solutions. The best-performing combinations of relaxation factors are then given together with the obtained relaxation factor trends.

#### **2. Materials and Methods**

Three different flow systems, with both the "U" and the "Z" flow arrangements ("U": outlet on the same side of the flow system as the inlet, "Z": outlet on the opposite side), were used to generate test cases. Moreover, in two of these three flow systems, the mesh fineness was also varied (coarser and finer mesh). This yielded ten flow system configurations in total, with simplified, cuboid cell-only meshes of different sizes ranging from ~6000 cells to ~41,000 cells (see Table 1). The meshes were generated automatically by the employed benchmarking software (see further) using the key parameters listed in Table 1. Due to the cuboid nature of the meshes, cell sizes were, in all computational domains, governed primarily by how many cell faces comprised a tube cross section (coarser mesh: 1 face only, finer mesh: 2 × 2 faces) and by the utilized cell growth factor. Sample meshes are shown in Figure 1.


**Table 1.** Flow system parameters: *W*, *L*, and *H* denote width, length, and height, respectively, *p*r row pitch, and *p*t tube pitch.

As for boundary conditions, 0.5 kg/s of water at 300 K was fed into the inlet while pressure at the outlet was set to 101,325 Pa. All walls were adiabatic except for tube walls, where a specific heat flux of 15 kW/m2 was set when the energy equation was enabled. Steady state simulations were carried out using the same CFD setup as in [42], that is, the SIMPLEC pressure-velocity coupling [43] was employed together with the Power Law discretization scheme [44]. Standard scaled residual limits, i.e., 10−<sup>3</sup> for continuity and momentum and 10−<sup>6</sup> for energy, were used. Only the natural ordering of the variables was considered in this study.

**Figure 1.** Sample simplified (cuboid cell-only) meshes generated automatically by the benchmarking software using the growth factor of 1.15; clockwise from top left: (**a**) layout of the "Z"-arranged flow system 1 from Table 1; (**b**) the corresponding finer mesh (~25,000 cells, groups of red cells in the top view denote locations of tube cross-sections); (**c**) top view of the finer mesh of the distribution header from flow system 2 (inlet region being on the left); other parts of this mesh, as well as the remaining meshes, were generated analogously.

The SSOR preconditioning technique was paired with two widely used numerical solution methods, which were shown by an earlier study [42] to perform very well in simplified 3D CFD models. The conjugate gradient (CG) numerical solution method [45] was employed in the pressure correction equation. The momentum and energy equations, on the other hand, were solved using the bi-conjugate gradient stabilized numerical method with the minimization of residuals over *L*-dimensional subspaces (BiCGstab(*L*)) [46], with *L* = 1, 2, or 3. Performance of the ILU preconditioner (which is efficient, but the construction of the respective preconditioning matrix may not always be possible) was taken as the baseline.

The SSOR-preconditioned numerical solution methods were tested with various tuples of ωF, ω<sup>R</sup> = 0.1, 0.2, ... , 1.8, 1.9 in successive steps, depending on the obtained results. Promising combinations of ω<sup>F</sup> and ω<sup>R</sup> were then taken as pivots, and their square neighborhoods were tested further—that is, all combinations of ωF, ω<sup>R</sup> with the respective values being ω – 0.04, ω – 0.02, ω, ω + 0.02, and ω + 0.04 were evaluated except for the original pivot point (ωF, ωR). In order to be able to compare SSOR to the baseline (ILU), all the ten meshes were also evaluated using the ILU-preconditioned combinations of solution methods. In total, 50,620 CFD model setups were tested.

The same benchmarking simplified 3D CFD Java software application was used as in [42], and, therefore, the reader is kindly referred to this paper for details (please note that the software is not publicly available). The benchmarking procedure itself was almost the same, as well, with the only difference being that the numbers of warm-up and test runs were smaller for the larger meshes (see Table 2). Such a measure was necessary to keep the times required to complete the respective benchmarks within reasonable bounds. This did not introduce any problems, because with larger cell counts all Java initializations and compilations had been finished within much less warm-up runs, and, therefore, it was not needed to carry out many of them before the timing phase. Mean test-run computational time was then taken as the final performance metric. Unlike in [42], however, only one machine (Intel Xeon E5 2698 v4 CPU, 128 GB RAM) was used instead of two largely disparate ones. The reason for this simplification was that, as shown in the respective paper, single-core computational times proved to be virtually identical, irrespective of whether the machine was a high-performance server or a regular laptop with an ultra-low voltage CPU. Please see the paragraph titled Supplementary Materials on how to obtain the data set containing all the mean computational times together with other relevant information.


**Table 2.** Numbers of warm-up and test runs and test case limits.

Because this study targeted fast computation, two kinds of limits were set in the solution process as detailed in Table 2. The first one concerned the number of CFD solver iterations, while the second one applied to the actual computational time. Any combination of numerical solution methods and preconditioning techniques which exceeded at least one of these two limits was marked as failing to reach a solution. Additionally, since robustness is one of the factors that must be considered when evaluating the suitability of numerical solution methods, no user interventions (e.g., changes to the internal residual limits of the numerical methods, CFD relaxation factors, etc.) were allowed during a solution process.

#### **3. Results**

This section is split into four parts. The first part presents the results pertaining to flow-only simulations. The second part discusses flow and energy transport simulations, in which only the energy equation was solved using SSOR-preconditioned numerical methods. The third part, on the other hand, summarizes data obtained via benchmarks, where SSOR was used for both the momentum and the energy equations. The last part then compares the SSOR-related data with results corresponding to the cases where solely the ILU preconditioning technique was used.

#### *3.1. Flow-Only Simulations with SSOR-Preconditioned Momentum Equations*

As mentioned above, the pressure correction equation was always solved using CG. This solution method was first coupled with SSOR, while the momentum equations were solved using BiCGstab(3):ILU. Within the several initial benchmarks, however, it became clear that CG:SSOR is wholly unsuitable. No matter the actual relaxation factors ω<sup>F</sup> and ωR, the majority of test cases either failed (mostly due to divergence) or resulted in very long computational times. Those setups which did not fail featured ω<sup>F</sup> ≈ ω<sup>R</sup> and, what is more, their amount decreased with an increasing mesh size, as shown in Figure 2. This was most probably a consequence of the fact that SSOR is sensitive to the ordering of variables. Only the ILU preconditioning technique was, therefore, used for the pressure correction equation, from then on.

**Figure 2.** Mean computational times for three different meshes and various combinations of ω<sup>F</sup> and ωR, the pressure correction equation was solved using CG:SSOR and the momentum equations using BiCGstab(3):ILU; diagonal crosses indicate failing combinations of ω<sup>F</sup> and ωR. Please note that the colorbars have been adjusted so that the best (i.e., the most relevant) combinations of ωF, ω<sup>R</sup> are easily identifiable. Further, due to the necessary colorbar ranges, the respective lower limits have been set to zero even though the data start at higher values.

The shortest mean computational times obtained with BiCGstab(2) and BiCGstab(1) instead of BiCGstab(3) were, on average, ~80% and ~260% longer, respectively. Conversely, stability of the solution process was greater with these two methods, and fewer combinations of ω<sup>F</sup> and ω<sup>R</sup> resulted in failures (again, mostly because of divergence; see Figure 3). It is also evident from the figure that the best-performing combinations were around ω<sup>F</sup> = ω<sup>R</sup> = 1.0, while with ω<sup>F</sup> below 0.5 or above 1.5 the computational times were much longer, or solution failures occurred. As for the CFD setup involving BiCGstab(3) in particular, mean computational times started at 2.56 s for the smallest mesh and 117.33 s when the largest mesh was used.

The solution behaviormentioned above was also observed for thelargermeshes. Only BiCGstab(3):SSOR was, therefore, used to solve the momentum equations in the following flow-only benchmarks because of its superior performance. To generate these, the combinations of ω<sup>F</sup> and ω<sup>R</sup> obtained for the ten flow systems were sorted by mean computational time and the respective ordered sets were then used to find the most common combinations yielding the most favorable computational times. In other words, the best tuples (ωF, ωR) were taken as pivot points whose square neighborhoods were evaluated further.

**Figure 3.** Mean computational times for the smallest mesh and the cases when the pressure correction equation was solved using CG:ILU and the momentum equations using BiCGstab(3):SSOR, BiCGstab(2):SSOR, and BiCGstab(1):SSOR; diagonal crosses indicate failing combinations of ω<sup>F</sup> and ωR. Please note that the colorbars have been adjusted so that the best (i.e., the most relevant) combinations of ωF, ω<sup>R</sup> are easily identifiable.

Example plots resulting from such evaluations are shown in Figure 4. Both pertain to the same mesh—the smallest one in this particular case. The plot on the left (Figure 4a) shows mean computational times for all the pivots and their respective neighborhoods. The plot on the right (Figure 4b) displays a cropped area corresponding to ωF, ω<sup>R</sup> ∈ [0.5, 1.5], where the best-performing combinations of relaxation factors are located. The dotted lines in both these plots represent the trend obtained using the standard weighted least squares method. Because the goal was to minimize computational time, the weights for individual combinations (ωF, ωR) were calculated as *wi* = (*t*min/*t*i) 4, where *t*min denotes the minimum computational time observed with a specific mesh and *ti* denotes the mean computational time corresponding to the respective (*i*-th) combination of relaxation factors evaluated using this mesh. The fourth power of the computational time ratio instead of just the ratio itself was used to adequately limit the influence of combinations that yielded solutions in longer time frames. Weights for combinations leading to solution failures were set to zero.

**Figure 4.** Mean computational times for the neighborhoods of pivot points corresponding to the smallest mesh; the pressure correction equation was solved using CG:ILU and the momentum equations using BiCGstab(3):SSOR; please note that, for the sake of clarity, the time ranges have been severely limited in both plots. Again, the colorbars have been adjusted so that the best (i.e., the most relevant) combinations of ωF, ω<sup>R</sup> are easily identifiable.

Trends for all the mesh sizes as well as the overall relaxation factor trends are listed in Table 3 and shown in Figure 5. Although the *R*<sup>2</sup> values are relatively low, this is caused by the fact that the data featured many less-relevant points scattered over the (0, 2) × (0, 2) relaxation factor domain, which were assigned small, but still non-zero weights. In any case, the standard errors for the trend coefficients are quite reasonable, and it is obvious that all the trends are very similar. Considering the actual values of the coefficients *a* and *b* and the fact that the best combinations of relaxation factors

featured ω<sup>F</sup> ≈ 0.9, it follows that, when solving the momentum equations in flow-only scenarios, both SSOR sweeps should be slightly underrelaxed to gain the shortest computational time.


**Table 3.** Relaxation factor trends, ω<sup>R</sup> = *a*ω<sup>F</sup> + *b*, for individual mesh sizes and the overall relaxation factor trend for the solution of momentum equations using BiCGstab(3):SSOR; *R*<sup>2</sup> denotes the coefficient of determination and *SE*(*a*) and *SE*(*b*) the standard error values for the coefficients *a*, *b.*

**Figure 5.** Relaxation factor trends for the solution of momentum equations using BiCGstab(3):SSOR; the overall trend was obtained using the merged data set and identical weights.

#### *3.2. Flow & Energy Transport Simulations with SSOR-Preconditioned Energy Equation*

Based on the solution behavior observed in the flow-only scenarios, only CG:ILU was used to solve the pressure correction equation in the flow and energy transport simulations. Momentum equations were also preconditioned only with ILU, while SSOR was used just for the energy equation. Various combinations of BiCGstab(*L*) for the momentum and energy equations were tested first, and the two most suitable combinations were then evaluated in detail using square neighborhoods of promising relaxation factor tuples (i.e., the pivot points).

The best results overall were obtained using BiCGstab(3) and BiCGstab(1) for the momentum equations and the energy equation, respectively (see Figure 6). This setup was relatively robust, most probably because of the better stability and smoothness of convergence resulting from the use of BiCGstab(1). Figure 7 shows the second-best combination, featuring only BiCGstab(2). On average, this setup was ~38% slower, and more combinations of ω<sup>F</sup> and ω<sup>R</sup> resulted in solution failures. It can also be seen that all the suitable relaxation factor tuples were in a relatively small neighborhood of (1, 1). Furthermore, with the SSOR-preconditioned energy equation, a significantly larger percentage of failures than before was due to slow convergence (that is, the respective iteration limits were exceeded).

The data sets mentioned above were then combined with data sets obtained by evaluating square neighborhoods of the promising tuples of ω<sup>F</sup> and ω<sup>R</sup> to get the corresponding relaxation factor trends. Again, the standard weighted least squares method was used with the weights being calculated in the same manner as before. Because the best setup and the second-best one (featured in Figures 6 and 7) were, at least in some cases, on par, the trends were calculated for both of them (see Table 4 and Figure 8). It is of note here that all benchmarks involving the largest mesh and BiCGstab(1) had failed. This suggests that the respective solution method simply is too slow when combined with SSOR. In any case, the best results were obtained with ω<sup>F</sup> around 1.2 or 1.1 when BiCGstab(1) or BiCGstab(2) were used, respectively. This means that, given the calculated relaxation factor trends, ω<sup>R</sup> should also

be a little above 1.0 (i.e., it is best to slightly overrelax both SSOR sweeps when solving the energy equation).

**Figure 6.** Mean computational times for three different meshes and various combinations of ω<sup>F</sup> and ωR, CG:ILU was used to solve the pressure correction equation, BiCGstab(3):ILU was used for the momentum equations, and BiCGstab(1):SSOR was used for the energy equation; diagonal crosses indicate failing combinations of ω<sup>F</sup> and ωR. Please note that the colorbars have been adjusted so that the best (i.e., the most relevant) combinations of ωF, ω<sup>R</sup> are easily identifiable.

**Figure 7.** Mean computational times corresponding to the meshes from Figure 6 and the cases when the momentum equations were solved using BiCGstab(2):ILU, and the energy equation was solved using BiCGstab(2):SSOR; diagonal crosses indicate failing combinations of ω<sup>F</sup> and ωR. Please note that the colorbars have been adjusted so that the best (i.e., the most relevant) combinations of ωF, ω<sup>R</sup> are easily identifiable.

**Table 4.** Relaxation factor trends, ω<sup>R</sup> = *a*ω<sup>F</sup> + *b*, for individual mesh sizes and the overall relaxation factor trends for the two discussed setups using either BiCGstab(1):SSOR or BiCGstab(2):SSOR to solve the energy equation; *R*<sup>2</sup> denotes the coefficient of determination, *SE*(*a*) and *SE*(*b*) denote the standard error values for the coefficients *a*, *b*, and "n/a" denotes the fact that all the respective benchmarks had failed, and thus the trend could not be obtained.


**Figure 8.** Relaxation factor trends for the two discussed setups using either BiCGstab(1):SSOR or BiCGstab(2):SSOR to solve the energy equation; the overall trends were obtained using the merged data sets and identical weights.

#### *3.3. Flow & Energy Transport Simulations with SSOR-Preconditioned Momentum and Energy Equations*

The last set of SSOR benchmarks involved both the momentum and the energy equations being preconditioned using this technique. However, because of a significantly larger relaxation factor domain (four factors had to be chosen instead of two), only ωF, ω<sup>R</sup> = 0.7, 0.8, ... , 1.3 were evaluated, i.e., only the subdomain where the best-performing relaxation factor quadruples were expected to lie. Additionally, only the "Z"-arranged flow system meshes, and the two setups identified in Section 3.2 as the most promising ones, were considered. Such a reduction led to a decrease in the number of combinations to be evaluated from more than 2.6 million (10 meshes × 2 setups × 130,321 factor quadruples) to ~24 thousand (5 meshes × 2 setups × 2,401 factor quadruples). This still provided enough information to get a general sense of how solution processes would likely behave.

The respective benchmarks generally resulted in computational times and solution failure percentages comparable to those reached when just the energy equation was preconditioned using SSOR. As before, the failures mostly occurred due to slow convergence or—less often—because of divergence. Only with rare combinations of SSOR relaxation factors were the solution processes so slow that the respective time limit was exceeded.

The best-performing combinations of relaxation factors are listed in Table 5. It can be seen that with almost all meshes, the setup involving only BiCGstab(2) resulted in markedly longer computational times. Additionally, it should be noted that there were other combinations of factors providing similar numerical performance, but all of them were clustered around the values mentioned in the table.


**Table 5.** Combinations of relaxation factors resulting in the shortest mean computational times when both the momentum and the energy equations were preconditioned using SSOR; setup "B3/B1" denotes the case when BiCGstab(3) was used to solve the momentum equations and BiCGstab(1) the energy equation, while setup "B2/B2" corresponds to only BiCGstab(2) being used for both these equation types.

Visualizing the obtained data in one plot per data set is not possible because that would require four-dimensional plots. One could, however, fix the momentum relaxation factor tuple to, e.g., the values from the best-performing quadruple mentioned in Table 5, and then plot the corresponding two-dimensional energy relaxation factor map (or vice versa). Examples of such plots are shown in Figures 9 and 10.

**Figure 9.** Two-dimensional plots of mean computational times obtained for the smallest "Z"-arranged mesh with (**a**) the energy relaxation factor tuple fixed to (ωF, ωR) = (1.2, 0.7) and (**b**) the momentum relaxation factor tuple fixed to (ωF, ωR) = (1.3, 1.1); BiCGstab(3):SSOR was used to solve the momentum equations and BiCGstab(1):SSOR the energy equation. Please note that the colorbars have been adjusted so that the best (i.e., the most relevant) combinations of ωF, ω<sup>R</sup> are easily identifiable.

**Figure 10.** Two-dimensional plots of mean computational times obtained for the smallest "Z"-arranged mesh with (**a**) the energy relaxation factor tuple fixed to (ωF, ωR) = (0.9, 1.2) and (**b**) the momentum relaxation factor tuple fixed to (ωF, ωR) = (1.0, 0.9); BiCGstab(2):SSOR was used to solve both the momentum equations and the energy equation; diagonal crosses indicate failing combinations of ω<sup>F</sup> and ωR. Please note that the colorbars have been adjusted so that the best (i.e., the most relevant) combinations of ωF, ω<sup>R</sup> are easily identifiable.

Because, here, one must choose two relatively independent relaxation factor tuples, it is best to generate two trends for each combination of numerical solution methods. This can be done by "flattening" the four-dimensional data to two dimensions (while still considering all the data points). In other words, if one sought, e.g., the momentum relaxation factor trend, one would disregard the energy-related part of the relaxation factor quadruple and thus have multiple data points with different weights (calculated just as before) for each momentum relaxation factor tuple. The respective trends would then, again, be calculated via the standard weighted least squares method (see Tables 6 and 7 and Figures 11 and 12). From the results, it follows that for the momentum equations, the SSOR forward sweeps should generally be carried out with ω<sup>F</sup> between ca. 1.1 and 1.3, while the backward sweeps should use ω<sup>R</sup> ≤ ωF. As for the energy equation and BiCGstab(1), the forward sweep should, again, be slightly overrelaxed (ω<sup>F</sup> up to ca. 1.2) and ω<sup>R</sup> ≤ ωF, while with BiCGstab(2) ω<sup>F</sup> should be around 1.0 (i.e., without any forward sweep relaxation) and ω<sup>R</sup> ≥ ωF.


**Table 6.** Relaxation factor trends, ω<sup>R</sup> = *a*ω<sup>F</sup> + *b*, for individual mesh sizes and the overall relaxation factor trends for the setup where BiCGstab(3):SSOR was used to solve the momentum equations and BiCGstab(1):SSOR the energy equation; *R*<sup>2</sup> denotes the coefficient of determination and *SE*(*a*) and *SE*(*b*) the standard error values for the coefficients *a*, *b*.

**Table 7.** Relaxation factor trends, ω<sup>R</sup> = *a*ω<sup>F</sup> + *b*, for individual mesh sizes and the overall relaxation factor trends for the setup where BiCGstab(2):SSOR was used to solve both the momentum equations and the energy equation; *R*<sup>2</sup> denotes the coefficient of determination and *SE*(*a*) and *SE*(*b*) the standard error values for the coefficients *a*, *b*.


**Figure 11.** Relaxation factor trends for (**a**) momentum and (**b**) energy obtained with BiCGstab(3):SSOR and BiCGstab(1):SSOR as the solution methods for the momentum equations and the energy equation, respectively; the overall trends were obtained using the merged data sets and identical weights.

#### *3.4. Comparison to ILU-Only Simulations*

In order to be able to assess the potential benefit of using SSOR, the meshes listed in Table 1 were also evaluated via the two combinations of numerical solution methods from Table 5, but only with the ILU preconditioning technique being employed. The results, summarized in Table 8 and visually compared in Figure 13, suggest that utilizing SSOR leads to at least a ~23% increase (on average) in computational time.

**Figure 12.** Relaxation factor trends for (**a**) momentum and (**b**) energy obtained with BiCGstab(2):SSOR as the solution method for both the momentum equations and the energy equation; the overall trends were obtained using the merged data sets and identical weights.

**Table 8.** Comparison of the mean computational times (s) obtained when solely the ILU preconditioning technique was used, and the overall best-case mean computational times reached with the momentum equations and/or the energy equation being preconditioned with SSOR.



**Figure 13.** Comparison of the mean computational times obtained when solely the ILU preconditioning technique was used, and the overall best-case mean computational times reached with the momentum equations and/or the energy equation being preconditioned with SSOR.

#### **4. Discussion**

The aim of this study was to establish whether, in the case of simplified CFD models, the SSOR preconditioning technique can be a viable replacement for ILU. From the obtained data, it follows that SSOR should not be used in conjunction with CG to solve the pressure correction equation. When applied to the momentum and/or energy equations, computational times tend to be significantly longer even when the relaxation factors are chosen favorably (for the cases evaluated in this study, the increase was at least ~23% on average). However, because the SSOR preconditioning matrix can always be constructed, the respective techniques could be used in conjunction with ILU as a fallback option. From an engineering point of view, this would mean that ILU would be employed by default, and only in case of numerical issues would the CFD solver try to reach a converged solution using

SSOR. Such an approach would capitalize on the efficiency of ILU while maintaining reasonable numerical robustness due to the possibility of falling back to a technique with guaranteed existence of the preconditioning matrix. The resulting models would, ultimately, be much more suitable for implementation in optimization algorithms or for other use cases where large batches of simulations must be carried out without user intervention.

The best-performing combinations of numerical solution methods and SSOR forward and backward relaxation factors differ according to whether energy transport is included in the model or not. In the flow-only scenario, the momentum equations should preferably be solved using BiCGstab(3), with ω<sup>F</sup> ≈ 0.9 and ω<sup>R</sup> slightly less than ωF, that is, both SSOR sweeps should be a little underrelaxed. Computational times obtained using other variants of BiCGstab(*L*) proved to be at least 80% longer. If also the energy transport is included and only the energy equation is preconditioned using SSOR, it is best to solve the momentum equations using BiCGstab(3) and the energy equation using BiCGstab(1). Here, both SSOR sweeps should be slightly overrelaxed, with ω<sup>F</sup> ≈ 1.2 and ω<sup>R</sup> ≈ 1.1. Similar performance can in some cases be obtained by employing BiCGstab(2) for both types of equations with the energy SSOR sweeps being overrelaxed using ω<sup>F</sup> ≈ ω<sup>R</sup> ≈ 1.1; however, a much greater solution failure probability can then be expected. If SSOR is utilized for both the momentum and the energy equations, then it is, again, preferable to use the combination of BiCGstab(3) and BiCGstab(1). The respective forward sweeps should be a little overrelaxed (ω<sup>F</sup> between ca. 1.1 and 1.3 for the momentum, and up to ca. 1.2 for the energy), while the backward sweeps should feature ω<sup>R</sup> slightly lower than ωF. The best-case computational times obtained with BiCGstab(2) proved to be up to ~67% longer and, therefore, the use of this numerical solution method is discouraged in this scenario.

**Supplementary Materials:** The mean computational times together with other relevant information are available online at http://www.mdpi.com/1996-1073/12/12/2438/s1.

**Funding:** This research was funded by the Czech Republic Operational Programme Research, Development, and Education, Priority 1: Strengthening capacity for quality research, grant No. CZ.02.1.01/0.0/0.0/15\_003/ 0000456 "Sustainable Process Integration Laboratory—SPIL".

**Conflicts of Interest:** The author declares no conflict of interest.

#### **References**


© 2019 by the author. 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 (http://creativecommons.org/licenses/by/4.0/).

*Article*

## **One Convenient Method to Calculate Performance and Optimize Configuration for Annular Radiator Using Heat Transfer Unit Simulation**

## **Zhe Xu 1,2,\*, Yingqing Guo 1, Huarui Yang 2, Haotian Mao 1, Zongling Yu <sup>2</sup> and Rui Li <sup>2</sup>**


Received: 4 December 2019; Accepted: 3 January 2020; Published: 5 January 2020

**Abstract:** In order to calculate heat transfer capacity and air-side pressure drop of an annular radiator (AR), one performance calculation method was proposed combining heat transfer unit (HTU) simulation and plate-and-fin heat exchanger (PFHX) performance calculation formulas. This method can obtain performance data with no need for meshing AR as a whole, which can be convenient and time-saving, as grid number is reduced in this way. It demonstrates the feasibility of this performance calculation method for engineering applications. In addition, based on the performance calculation method, one configuration optimization method for AR using nondominated sorted genetic algorithm-II (NSGA-II) was also proposed. Fin height (FH) and number of fins in circumferential direction (NFCD) were optimized to maximize heat transfer capacity and minimize air-side pressure drop. Three optimal configurations were obtained from the Pareto optimal points. The heat transfer capacity of the optimal configurations increased by 22.65% on average compared with the original configuration, while the air-side pressure drop decreased by 33.99% on average. It indicates that this configuration optimization method is valid and can provide a significant guidance for AR design.

**Keywords:** annular radiator; performance calculation; configuration optimization; heat transfer unit; plate-and-fin heat exchanger; nondominated sorted genetic algorithm-II

#### **1. Introduction**

Heat exchangers, which can transfer heat between two fluids of different temperatures, are widely used in several industrial applications, just like aerospace engineering, petrochemical progress, nuclear power plants, oil refining, and so on [1]. As performance calculation and configuration optimization are main steps in the process of heat exchanger design, numerous up-to-date technologies have been applied in this research area from various points of view.

Theoretical analysis is a common method to calculate heat exchanger performance. As one kind of well-known theoretical analysis method, the Bell–Delaware method is widely used in performance calculation of heat exchangers. In the research of L. Yi et al. [2], one cooling water heat exchanger was designed preliminarily using the Bell–Delaware method, and the design results were verified by HTRI 6.0 commercial software. Some mathematical relationship formulas for performance calculation of shell-and-tube heat exchanger (STHX) with helical baffles were proposed by B. A. Abdelkader et al. [3] based on Bell–Delaware method. In addition, in the research of B. A. Abdelkader et al. [4], Kern method, Bell–Delaware method, and flow stream analysis (Wills Johnston) methods were applied, respectively, to predict both heat-transfer coefficient and pressure drop on the shell side of a heat exchanger. The calculation results were compared with the experimental data, and the comparison results showed that the Bell–Delaware method was the most accurate method.

Computational fluid dynamics (CFD) technology is also widely used in performance analysis of heat exchangers. The performances of STHXs with and without baffles were compared through CFD calculation program OpenFOAM-2.2.0 by E. Pal et al. [5]. In the research of J. Du et al. [6], a midtemperature gravity heat pipe exchanger was taken as the research object, and the effects of different operating parameters and fin parameters on heat transfer performance were studied using Fluent software. A multitube tank was proposed by M. Ramadan et al. [7] as a heat exchanger. The performance of it was analyzed using CFD simulation code, and the best scenario among three different configurations was obtained by an optimization procedure. Y. Yang et al. [8] studied heat transfer and flow characteristics in a type of plate heat exchanger by numerical simulation, and the correlations of single-phase heat transfer coefficient and friction coefficient were presented. In the research of plate-and-fin heat exchangers (PFHXs), the effects of inlet header configuration on fluid flow maldistribution and the effects of top bypass flow of fins on thermal performance were studied by A. Raul et al. [9] and H. Cai et al. [10] using CFD simulation, respectively. In addition, based on CFD technology, the performances of STHXs separately with segmental baffles [11,12], trefoil-hole baffles [12,13], and helical baffles [11,12,14] were analyzed. Especially in the research of A. El Maakoul et al. [12], the simulation results of these three kinds of heat exchangers were compared, and STHX with helical baffles was found to be the one that had the best comprehensive performance.

Combining theoretical analysis with CFD simulation can be another useful way to research heat exchanger performance. R. Amini et al. [15] compared the CFD simulation results with the calculation results by the Bell–Delaware method to validate the accuracy of the simulation method. D.M. Godino et al. [16] calculated the heat transfer coefficient of the preheater separately by Bell–Delaware method and Kern method, and the calculation results were both compared with the CFD simulation data to analyze the accuracy. X. Gu et al. [17] came up with periodic whole cross-section computation models to obtain performance data of segmental baffle heat exchanger, shutter baffle heat exchanger, and trapezoid-like tilted baffle heat exchanger, and the reliability of the method was verified by comparing the simulation data with the calculation results using the Bell–Delaware method. I. Milcheva et al. [18] improved the traditional Bell–Delaware method and introduced an enhancement factor to calculate the performance of a STHX with double-segmental baffles. The calculation results were compared with the CFD simulation data to validate the effectiveness of this method.

As for configuration optimization for heat exchanger, nondominated sorted genetic algorithm-II (NSGA-II) is commonly used in some research. The Bell–Delaware procedure and the ε-NTU method were applied in STHX performance estimation by S. Sanaye et al. [19] and NSGA-II was used to maximize heat transfer coefficient and minimize pressure drop simultaneously. M. Chahartaghi et al. [20] combined NSGA-II with entransy dissipation theory to minimize the entransy dissipation numbers separately caused by thermal conduction and fluid friction for STHX. Z. Xu et al. [21] calculated the performances of two kinds of STHXs by theoretical formulas and optimized their structural parameters using NSGA-II, respectively. In order to optimize the configuration of a PFHX with offset strip fins, NSGA-II combined with the ε-NTU method was adopted in the research of R. Song et al. [22]. In addition, as a popular optimization method, NSGA-II combined with response surface method was also used in the configuration optimizations of STHX with helical baffles [23–25], spiral-wound heat exchanger [26], helically coiled tube heat exchanger [27], torsional flow heat exchanger [28], and triple concentric-tube heat exchanger [29]. Except for NSGA-II, some other novel optimization algorithms were also proposed to optimize configuration parameters of different heat exchangers, such as firefly algorithm [30], Tsallis differential evolution algorithm [31], bat algorithm [32], Taguchi method [33], particle swarm optimization (PSO) [34], cohort intelligence algorithm [35], tree traversal method [36], surrogate-based optimization algorithm [37], wale optimization [38], topology optimization [39], Jaya algorithm [40], and so on.

As annular radiator (AR) is a neoteric heat exchanger, the research about its thermal-hydraulic performance calculation and configuration optimization, which are badly in need of design processes, are rare. In this paper, a feasible and reliable method for performance calculation and configuration optimization based on heat transfer unit (HTU) simulation and NSGA-II was proposed, which can conveniently obtain AR heat transfer capacity and air-side pressure drop, while avoiding the problem of the huge amount of grids generated by meshing AR as a whole directly, and getting the optimized fin height (FH) and number of fins in circumferential direction (NFCD), which are a trade-off on maximizing heat transfer capacity and minimizing air-side pressure drop.

#### **2. Method**

#### *2.1. Physical Model and Normal Operating Conditions*

AR is an air–liquid heat exchanger, which mainly consists of liquid channels, annular substrate, and fins. As depicted in Figure 1, liquid flows through the channels inside the annular substrate, and air passes by the fins, which are averagely distributed on the circular inner side of the annular substrate. So, the heat can be transferred from the liquid of higher temperature to the air of lower temperature through fins and annular substrates.

The material of AR is aluminum alloy, and the normal operating conditions are shown in Table 1. In this case, the working fluid of liquid is lubricating oil.

**Figure 1.** Schematic diagram of annular radiator (AR).



#### *2.2. Performance Calculation Method*

#### 2.2.1. Structural Equivalence

As AR is structurally similar to PFHX with fins of rectangular straight wave, the performance calculation formulas of this kind of heat exchanger were applied to the calculation of AR heat transfer capacity in this method. Basic formulas of PFHX heat transfer capacity calculation are given in Appendix A. The structure of the equivalent PFHX is depicted in Figure 2. The main configuration parameters of AR and the corresponding PFHX are shown in Figure 3 and Table 2.

**Figure 2.** Structure of the equivalent plate-and-fin heat exchanger (PFHX).

**Figure 3.** Configuration parameter correspondence between AR and PFHX: (**a**) configuration parameters of AR; (**b**) corresponding configuration parameters of PFHX.


**Table 2.** Configuration parameter values.

#### 2.2.2. HTU Simulation Method

A traditional way to obtain the performance data of a heat exchanger is through emulating it as a whole directly. While, for AR simulation, a huge number of grids will be generated by this traditional method in the process of meshing. That is extremely time-consuming for solving and is beyond the available computing resources. In order to avoid this problem, HTU simulation, rather than overall simulation of AR, was used in this method. When an AR is working, the air above the fins of the AR does not pass through the clearance between fins, which means that it is unscientific if the air inlet flow rate of AR is used in the PFHX performance calculation formulas directly. So, the percentage of the air calculated in the formulas needs to be obtained through HTU simulation in this method. This percentage is defined by Equation (1):

$$k = q\_{fin} / q\_{HTU} \,\text{.}\tag{1}$$

In Equation (1), *qfin* and *qHTU* are the volume flow rates of the air passing through the clearance between fins of HTU and the air of HTU air-side, respectively, and *k* reflects the percentage of the air that can be calculated in PFHX performance calculation formulas.

HTUs are obtained through dividing AR, and two kinds of HTUs can be obtained according to the symmetry of AR, as shown in Figure 4. As HTU A and HTU B are structurally similar, and air flows through HTU A first, HTU A was chosen to be the only one kind of HTU to be emulated in this method. In this case, the central angle of HTU was 2 degrees.

**Figure 4.** Two kinds of heat transfer units (HTUs): (**a**) HTU A; (**b**) HTU B.

In the process of CFD simulation, ANSYS commercial software was adopted. There were two fluid domains, including air domain and liquid domain, and one solid domain in the CFD simulation of HTU, as shown in Figure 5. Unstructured grids were adopted, and the normal operating conditions of AR were used in HTU simulation. In order to improve simulation accuracy, grid independence validation was carried out. In Figure 6, Δ*P* represents air-side pressure drop of HTU. As shown in this figure, the simulation deviations of *k* and Δ*P* were both in the acceptable range when the grid number increases above the point of 24,020,699. So, the mesh generation settings of this point were used.

**Figure 5.** HTU computational domains.

**Figure 6.** Grid independence validation.

In this case, the air flow rate was less than Ma0.3, and the temperature deviation between air inlet and outlet was small. So, like liquid, air was treated as a Newtonian and incompressible fluid with constant physical property parameters. In addition, in the process of simulation, fluid flow and heat transfer process were considered as turbulent and in steady-state, and fouling resistance was neglected.

In this method, a realizable *k-*ε turbulence model and the default constant values were adopted. Velocity-inlet and pressure-outlet were separately used as inlets and outlets of both air and liquid. In the process of solving this kind of fluid flow heat transfer problem, the SIMPLE (semi-implicit method for pressure-linked equations) algorithm is the most widely used. The core of the algorithm is to use continuous equations and momentum equations to construct an approximate pressure correction equation on staggered grids to calculate pressure field and correct velocity. The SIMPLE algorithm is very useful in engineering applications. However, the convergence velocity of the SIMPLE algorithm is not fast. Compared with the SIMPLE algorithm, the SIMPLEC (SIMPLE-consistent) algorithm, which is the improved version of SIMPLE algorithm, can achieve a higher rate of convergence as it synchronizes speed field improvement process with pressure field improvement process [41]. Thus, in this study, SIMPLEC algorithm, as one built-in solution algorithm in ANSYS Fluent commercial software, was chosen to solve the simulation problem. The second order upwind was used for the momentum, turbulent kinetic energy, turbulent dissipation rate, and energy. The default convergence criterion was adopted, which is that the normalized residuals are less than 1 <sup>×</sup> <sup>10</sup>−<sup>6</sup> for energy equation and 1 <sup>×</sup> <sup>10</sup>−<sup>3</sup> for the other equations. Iteration number was set as 2000.

#### 2.2.3. Method Procedure

AR performance calculation in this method included heat transfer capacity calculation and air-side pressure drop calculation. The main procedures are shown as follows:


$$q\_{\upsilon} = kq\_{\upsilon}^{\prime} \tag{2}$$

In Equation (2), *qv* is the air inlet volume flow rate of AR used in PFHX performance calculation formulas and *qv*' is the air inlet volume flow rate of AR.

5. Calculate air-side pressure drop of AR. As HTU is obtained by cutting AR along its symmetry plane, air-side pressure drop of AR can be computed by Equation (3):

$$
\Delta P = \text{2} \Delta P'\tag{3}
$$

In Equation (3), Δ*P* is air-side pressure drop of AR.

	- (a) Assume one heat transfer capacity for AR.
	- (b) Calculate the mean temperature of liquid/air-side under the assumptive heat transfer capacity.
	- (c) Calculate the physical property parameters (such as fluid density, kinematic viscosity, etc.) of liquid/air-side under the mean temperature of the corresponding side.
	- (d) Calculate the effective heat transfer area values and the heat transfer coefficients of liquid-side and air-side.
	- (e) Calculate NTU (number of transfer units), heat transfer efficiency, and heat transfer capacity.
	- (f) Give the value of the calculated heat transfer capacity to the assumptive heat transfer capacity, and go to step (b) until the deviation of them is within acceptable limits.
	- (g) Output the calculated heat transfer capacity.

The flowchart of this method is shown in Figure 7.

**Figure 7.** Flowchart of AR performance calculation method.

#### *2.3. Configuration Optimization Method*

In this research, in consideration of the performances of AR, heat transfer capacity *Q* and air-side pressure drop Δ*P* were chosen as two objective functions, and FH and NFCD were regarded as two design parameters. Thus, the multiobjective optimization problem was be formulated as Equation (4):

$$\begin{cases} \text{Min} & -Q, \Delta P \\ \text{S.t.} & \text{7mm} \le FH \le 16 \text{mm} \& FH \in Z \\ & 238 \le NFCD \le 544 \& NFCD \in Z \end{cases} \tag{4}$$

For a multiobjective optimization problem, it is always hard to find a solution that is absolutely optimal. However, through some evolutional algorithms, a Pareto optimal front that contains a series of optimal solutions can be obtained. NSGA-II is one kind of evolutional algorithm that is based on a genetic algorithm for multiobjective optimization. One or more satisfactory solutions can be selected by some criteria from the Pareto optimal front obtained by NSGA-II [24]. It incorporates elitism [42], and no sharing parameter needs to be chosen a priori [43]. As the nondominated sorting method is used as the ranking scheme in this algorithm, the convergence velocity of NSGA-II is faster than the traditional Pareto ranking method. Besides, as the constraint handling method also uses a nondominance principle as the objective, which guarantees that the feasible solutions are always ranked higher than the unfeasible solutions, penalty functions and Lagrange multipliers are not needed in NSGA-II [44]. So, due to these advantages, and in order to optimize these two configuration parameters shown in Equation (4) under the normal operating conditions shown in Table 1, the multiobjective genetic algorithm NSGA-II is adopted in this configuration optimization method to obtain the optimal solutions. The main procedures are shown as follows:


$$k = f\_1(FH, NFCD),\tag{5}$$

$$
\Delta P' = f\_2(FH, NFCD) \tag{6}
$$

	- (a) Initialize NSGA-II parameters, including population size, generation number, and so on, and generate a random population in the constraints of design parameters.
	- (b) Calculate objective function values for each chromosome of the population using AR performance calculation method and the functional relationships obtained in step 3.
	- (c) Sort chromosomes based on non-domination and crowding distance. In this method, the crowding distance is compared only if the ranks for both chromosomes are the same.
	- (d) Choose the chromosomes that are fit for reproduction as the parents of the next generation using tournament algorithm.
	- (e) Generate children by crossover and mutation, and calculate their objective function values using AR performance calculation method and the functional relationships obtained in step 3.
	- (f) Combine parents and children, and sort them based on nondomination and crowding distance.
	- (g) A new generation is extracted based on ranking.

Other details of NSGA-II are referred in [43]. The flowchart of this method is shown in Figure 8.

**Figure 8.** Flowchart of AR configuration optimization method.

#### **3. Results and Discussion**

#### *3.1. Performance Calculation Method Validation*

HTU simulation results are shown in Figure 9. As depicted, when air passes by fins, only a small part of air passes through the clearance between fins, and the velocity of the air goes above fins is higher than that of the other. So, the calculation of *k* is significant. Based on the plane perpendicular to the air flow direction and located in the middle place of the finned area, the calculated *k* was 11.5%. That means only 11.5% of the air passes through the clearance between fins and can be used in PFHX performance calculation formulas.

In order to validate the accuracy of the AR performance calculation method, the calculation results were compared with the experimental data. The experimental conditions, the experimental schematic diagram, and the testing equipment for AR are shown in Table 1 and Figure 10; Figure 11, respectively. The experimental equipment included an air supply system, oil circulation system, and measuring system. An air supply system can provide AR the airflow of the required flow rate, temperature, and pressure. The maximum air supply capacity of it is 8600 kg/h, which can meet the requirement of this test. An oil circulation system can heat the lubricating oil to the required temperature and then supply it to AR. The measuring system includes some instruments distributed in the inlets and outlets of both air and oil. As the radius of AR is large, five pressure sensors and five temperature

sensors were averagely distributed in the air inlet of AR, and the same is true for the air outlet. The experimental temperature or pressure value of the air inlet or outlet of AR was computed by averaging the corresponding five measured values. The results of the comparison between calculation results and experimental data are given in Figure 12.

**Figure 9.** Velocity streamline of HTU air-side.

**Figure 10.** Experimental schematic diagram.

**Figure 11.** Testing equipment for AR.

**Figure 12.** Comparison between calculation results and experimental data.

As depicted in Figure 12, Δ*P* obtained by this method was consistent with the experimental data, while the error of *Q* was higher than Δ*P*. The deviations of *Q* and Δ*P* were 20.7–22.4% and 5.5–6.9% with the average errors of 21.5% and 6.2%, respectively. The reason for this is that, as shown in Equation (3), half of Δ*P* is directly obtained by CFD simulation, which is accurate. However, the calculation of *Q* using PFHX performance calculation formulas only can consider the heat transfer of the air that flows through the clearances of fins, while the air that flows above fins also participates in the heat transfer process. Even though the heat transfer effect of the air flowing above fins is very small, it can still influence the precision of *Q* calculation and make the calculated values by this method lower than the experimental data. Thus, if a highly accurate *Q* result is needed, this part of air also needs to be considered, which is not realizable through using PFHX theoretical formulae directly. AR needs to be emulated as a whole directly to meet the high accuracy requirement, which is too time-consuming to realize in engineering applications. In the design process of engineering applications, the precision of this method is enough for AR performance prediction, and it is more convenient and time saving. Hence, it can be concluded that this method is feasible for engineering applications, and can be used in the configuration optimization process in this research.

Ignoring factors of this method itself, there are also some other reasons that can cause the differences between the calculation results and the experimental data, just like deviations of formulas themselves, simplification of simulation model, and unavoidable experimental errors.

#### *3.2. Design Parameter E*ff*ects and Configuration Optimization Results*

#### 3.2.1. Functional Relationships Fitting

In this case, 16 sets of HTU simulation data under different values of FH and NFCD were averagely obtained within design parameter constraints, and the surfaces of fitting were obtained using least square method and polynomial fitting based on these data. The fitting surfaces are shown in Figures 13 and 14, and the functional relationships can be represented by Equations (7) and (8).

$$\begin{array}{l}k = 0.03211 + 0.03009\mathbf{x} - 0.0001297\mathbf{y} + 0.0002819\mathbf{x}^2 - 3.732 \times 10^{-5}\mathbf{x}\mathbf{y} \\ + 4.445 \times 10^{-8}\mathbf{y}^2 - 8.684 \times 10^{-6}\mathbf{x}^3 + 1.589 \times 10^{-6}\mathbf{x}^2\mathbf{y} - 2.16 \times 10^{-8}\mathbf{x}\mathbf{y}^2 \end{array} \tag{7}$$

$$\begin{aligned} \Delta P &= -0.6039 + 0.1434\mathbf{x} + 0.004469y - 0.009944\mathbf{x}^2 - 0.000316\mathbf{x}y\\ -1.103 \times 10^{-5}y^2 &+ 0.000389\mathbf{x}^3 - 8.93 \times 10^{-6}\mathbf{x}^2y + 1.378 \times 10^{-6}\mathbf{x}y^2 \end{aligned} \tag{8}$$

Sum of squares due to error (SSE) and the coefficient of determination *R*<sup>2</sup> are significant goodness of fit criteria to evaluate the accuracy of fitting function [45]. The closer the SSE and *R*<sup>2</sup> are to 0 and 1, respectively, the better is the fitting function. The SSEs of Equations (7) and (8) are 1.9739 <sup>×</sup> 10−<sup>5</sup> and

0.0033, and the *R*<sup>2</sup> values of that are 0.9999 and 0.9993, respectively. Hence, the functional relationships can be obtained accurately.

**Figure 13.** Fitting surface for the percentage of the air that can be calculated in PFHX performance calculation formulas (*k*).

**Figure 14.** Fitting surface for air-side pressure drop of AR (Δ*P*).

#### 3.2.2. Design Parameter Effects

The effects of design parameters on *Q* and Δ*P* were analyzed based on Equations (7) and (8) and the AR performance calculation method. The effects of FH are represented in Figure 15, while NFCD was 544. As depicted, while FH increased from 7 mm to 16 mm, *Q* and Δ*P* increased by 9.66–26.26 kW and 1.25–4.87 kPa, respectively. The growth of FH led to the increase of heat transfer area, which can lead to the increase of *Q*. At the same time, as FH rose, the influence area of fins rose, and the resistance influence of fins was enhanced. Thus, Δ*P* increased continuously.

The effects of NFCD are shown in Figure 16, while FH was 10 mm. As depicted, while NFCD increased from 238 to 544, Δ*P* increased by 0.96–2.53 kPa, and *Q* increased from 16.71 kW first, reaching its highest point of 19.77 kW when NFCD was 408, and then decreased to 17.25 kW. The increase of NFCD indicates the growth of the heat transfer area, which results in the increase of *Q* at first. Meanwhile the increase of NFCD can also lead to the decrease of the clearance space between two adjacent fins, which means that air is increasingly harder to flow through these clearances and the

velocity of air that flows through the air channel without fins increases. So, Δ*P* increased continuously in this process, and *k* decreased continuously. The decrease of *k* will weaken the heat exchange capability of AR. After NFCD increases to a certain level, the decrease of *k* can lead to the decline of *Q,* even though heat transfer area is still increasing. Thus, the curve of *Q* rose first and then fell after the point.

**Figure 15.** Heat transfer capacity (*Q*) and Δ*P* separately versus fin height (FH).

**Figure 16.** *Q* and Δ*P* separately versus number of fins in circumferential direction (NFCD).

#### 3.2.3. Optimization Results

In order to consider the comprehensive performance of AR, the multiobjective configuration optimization method driven by NSGA-II was conducted. The conflicting optimization objectives were set as the minimization of −*Q* and Δ*P* both. Population size, maximum iteration number, analog binary cross distribution index, and polynomial mutation distribution index were set as 100, 500, 20, and 20, respectively. The obtained Pareto optimal points are shown in Figure 17.

**Figure 17.** Pareto optimal points.

As depicted in Figure 17, both the values of −*Q* and Δ*P* of some obtained optimal points were lower than that of the original point, which means that these obtained configurations had better comprehensive performances than the original one. In addition, three Pareto optimal solutions were chosen and are shown in Figure 17 and Table 3. The selection method of these three Pareto optimal solutions follows the following principles:



**Table 3.** Optimization performances comparison of AR.

Compared with the original point shown in Table 3, *Q* values of the point A, B, and C separately increased by 41.77%, 24.51%, and 1.66%, and Δ*P* values of them decreased by 4.35%, 38.74%, and 58.89%, respectively. It is clearly shown that *Q* of the optimal configurations increased by 22.65% on average, while Δ*P* decreased by 33.99% on average. Based on the above comparisons, it indicates that the proposed configuration optimization method is valid and feasible, and the comprehensive performance of AR can be enhanced by this method.

#### **4. Conclusions**

In this research, a performance calculation method for AR was proposed and verified by experiment. Heat transfer capacity and air-side pressure drop were calculated through combining HTU simulation and PFHX performance calculation formulas, rather than through emulating AR as a whole directly. So, the problem of the huge amount of grids generated by meshing AR can be

effectively avoided, which means that this method is convenient to use and can save calculation time and computing resources. It demonstrates the feasibility of this performance calculation method for engineering applications.

Based on this performance calculation method, a configuration optimization method for AR was also come up with using NSGA-II in this research. Heat transfer capacity maximization and air-side pressure drop minimization were regarded as two conflicting objectives, and FH and NFCD were set as two design parameters. A set of Pareto optimal solutions were obtained and some of them had better comprehensive performances than the original configurations. Three optimal solutions were chosen and compared with the original configuration. The comparison results illustrate that the heat transfer capacity of the optimal configurations increased by 22.65% on average compared with the original configuration, while the air-side pressure drop decreased by 33.99% on average. It indicates that this configuration optimization method is valid and can provide a significant guidance for AR design.

**Author Contributions:** Conceptualization, Z.X.; investigation, Z.X. and H.M.; methodology, Z.X.; project administration, Y.G. and H.Y.; resources, Z.X. and Z.Y.; software, Z.X.; supervision, Y.G. and H.Y.; validation, Z.X.; visualization, Z.X.; writing—original draft, Z.X.; writing—review and editing, Z.X., Y.G., Z.Y., and R.L. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received no external funding.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **Nomenclature**


*w* wall

#### **Appendix A**

The basic formula of heat transfer capacity calculation of PFHX is shown by Equation (A1) [46]:

$$Q = \eta \mathcal{W}\_{\text{min}} (T\_L - T\_A). \tag{A1}$$

As the numbers of air passes and liquid passes are 1 and 2, respectively, η is given by Equations (A2)–(A7) [47]:

$$\eta = \frac{\left(\frac{1 - C^\* \eta\_{\parallel}}{1 - \eta\_{\parallel}}\right)^2 - 1}{\left(\frac{1 - C^\* \eta\_{\parallel}}{1 - \eta\_{\parallel}}\right)^2 - C^\*} \tag{A2}$$

$$\eta\_{l} = 1 - \exp\left\{ \frac{NTlI^{0.22}}{\mathcal{C}^\*} [\exp\left(-\mathcal{C}^\* NTlI^{0.78}\right) - 1] \right\} \tag{A3}$$

$$NTU = \frac{UL}{W\_{\text{min}}} \tag{A4}$$

$$\mathcal{C}^\* = \frac{\mathcal{W}\_{\text{min}}}{\mathcal{W}\_{\text{max}}} \tag{A5}$$

$$\frac{1}{\text{LIA}} = \frac{1}{\eta\_{\text{L}}\alpha\_{\text{L}}A\_{L}} + R\_{\text{w}} + \frac{1}{\eta\_{\text{A}}\alpha\_{\text{A}}A\_{A}}\tag{A6}$$

$$R\_{\nu} = \frac{\delta\_p}{\lambda\_{\nu} A\_p} \tag{A7}$$

In Equation (A6), η*<sup>L</sup>* and η*<sup>A</sup>* reflect surface efficiencies of liquid-side and air-side, respectively, and *AL* and *AA* separately represent heat transfer areas of liquid-side and air-side. These variables can be calculated by Equations (A8)–(A15) [47,48]:

$$\eta\_L = 1 - \frac{A\_{fL}}{A\_L} (1 - \eta\_{fL}),\tag{A8}$$

$$\eta\_A = 1 - \frac{A\_{fA}}{A\_A} (1 - \eta\_{fA}) \tag{A9}$$

$$A\iota\_{\!\!\!\!\!\!\!\/} = A\_{\!\!\!\!\!\!\/} + A\_{\!\!\!\!\/\!\!\!\/} \tag{A10}$$

$$A\_A = A\_p + A\_{fA} \tag{A11}$$

$$\eta\_{fL} = \frac{\tanh(m\_L h\_L)}{m\_L h\_L} \tag{A12}$$

$$\eta\_{fA} = \frac{\tanh(m\_A h\_A)}{m\_A h\_A} \tag{A13}$$

$$m\_{L} = \sqrt{\frac{2\alpha\_{L}}{\lambda\_{fL}\delta\_{fL}}}\tag{A14}$$

$$m\_A = \sqrt{\frac{2\alpha\_A}{\lambda\_{fA}\delta\_{fA}}}\tag{A15}$$

*W* and α refer to thermal capacity rate and heat transfer coefficient, respectively, and they are given by Equations (A16) and (A17) [46,48]:

$$\mathcal{W} = q\_{\rm w} c\_{\rm p} \tag{A16}$$

$$\alpha = \frac{j\omega c\_p}{Pr^{2/3}}\tag{A17}$$

#### **References**


© 2020 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 (http://creativecommons.org/licenses/by/4.0/).

MDPI St. Alban-Anlage 66 4052 Basel Switzerland Tel. +41 61 683 77 34 Fax +41 61 302 89 18 www.mdpi.com

*Energies* Editorial Office E-mail: energies@mdpi.com www.mdpi.com/journal/energies

MDPI St. Alban-Anlage 66 4052 Basel Switzerland

Tel: +41 61 683 77 34 Fax: +41 61 302 89 18