An Integrated Three-Level Synergetic and Reliable Optimization Method Considering Heat Transfer Process, Component, and System

: Optimization of heat transfer systems (HTSs) beneﬁts energy e ﬃ ciency. However, current optimization studies mainly focus on the improvement of system design, component design, and local process intensiﬁcation separately, which may miss the optimal results and lack reliability. This work proposes a synergetic optimization method integrating levels of the local process, component to system, which could guarantee the reliability of results. The system-level optimization employs the heat current method and hydraulic analysis, the component level optimization adopts heuristic optimization algorithm, and the process level optimization applies the ﬁeld synergy principle. The introduction of numerical simulation and iteration provides the self-consistency and credibility of results. Optimization results of a multi-loop heat transfer system present that the proposed method can save 16.3% pumping power consumption comparing to results only considering system and process level optimization. Moreover, the optimal parameters of component originate from the trade-o ﬀ relation between two competing mechanisms of performance enhancement, i.e., the mass ﬂow rate increase and shape variation. Finally, the proposed method is not limited to heat transfer systems but also applicable to other thermal systems.


Introduction
Thermal energy is still of great importance in modern society. A heat transfer system (HTS) is widely used to transport thermal energy, its performance improvement benefits the energy efficiency improvement, and hence efficient and reliable optimization methods and strategies for HTS are necessary [1]. Current related studies can be divided into three levels from top to bottom, i.e., the system parameter level, component design level, and local process intensification level. In the system level, many studies apply various algorithms to optimize the structure and operating parameters such as heat transfer area and operation cost. Ponce-Ortega et al. [2] optimized a heat transfer network including streams with phase change using an mixed integer nonlinear programming MINLP model. Adjiman et al. [3] used the branch and bound method to minimize the total heat transfer area of a heat transfer network. Ravaganani et al. [4] combined the pinch-point analysis and genetic algorithm to minimize the total cost of a heat exchange system. Sameti and Haghighat [5] recently presented a detailed review for optimization researches focusing on district heating and cooling network. On the other hand, there are also contributions of system optimization employing different optimization principles, where entropy generation minimization principle (EGMP) and exergy destruction rate minimization are two obtained from simulation results can be used as boundary conditions. Moreover, optimization results often deviate from simulation conditions, and hence the simulation is required once more to determine these parameters. For the self-consistency of calculation, an iterating and updating strategy is required.
In this study, a synergetic optimization method for HTSs integrating all three levels of optimization is proposed. A typical HTS is taken as an example and optimized to minimize its pumping power consumption. The system-level optimization combines the heat current method and hydraulic analysis, the component design level optimization uses the heuristic algorithm, and the local heat transfer process optimization is performed numerically based on the FSP. An iterating and updating strategy of system parameters is realized using a multi-disciplinary optimization platform. Finally, optimization results are analyzed to show the advantages of the proposed synergetic method, and the mechanism of optimal parameters is briefly discussed.

Physical Model and Optimization Problem of the HTS
A typical multi-loop HTS is investigated in this work to present the proposed optimization method. As Figure 1 presents, the HTS contains a hot water reservoir, an evaporator, two counter-flow heat exchangers, and three variable speed pumps (VSPs). The tubes in the evaporator are elliptical tubes due to their heat transfer enhancement comparing to round tubes [36]. Temperatures of fluids are denoted by 1-3 according to the fluid loop, and subscripts h and c refer to the hot and cold fluid, respectively. The temperatures of hot water in the tank and refrigerant remain as T 1,h, and T e , respectively. In the system operation, VSPs drive working fluids, by which the heat duty is transferred from the hot water reservoir to the refrigerant. updating strategy is required.
In this study, a synergetic optimization method for HTSs integrating all three levels of optimization is proposed. A typical HTS is taken as an example and optimized to minimize its pumping power consumption. The system-level optimization combines the heat current method and hydraulic analysis, the component design level optimization uses the heuristic algorithm, and the local heat transfer process optimization is performed numerically based on the FSP. An iterating and updating strategy of system parameters is realized using a multi-disciplinary optimization platform. Finally, optimization results are analyzed to show the advantages of the proposed synergetic method, and the mechanism of optimal parameters is briefly discussed.

Physical Model and Optimization Problem of the HTS
A typical multi-loop HTS is investigated in this work to present the proposed optimization method. As Figure 1 presents, the HTS contains a hot water reservoir, an evaporator, two counterflow heat exchangers, and three variable speed pumps (VSPs). The tubes in the evaporator are elliptical tubes due to their heat transfer enhancement comparing to round tubes [36]. Temperatures of fluids are denoted by 1-3 according to the fluid loop, and subscripts h and c refer to the hot and cold fluid, respectively. The temperatures of hot water in the tank and refrigerant remain as T1,h, and Te, respectively. In the system operation, VSPs drive working fluids, by which the heat duty is transferred from the hot water reservoir to the refrigerant.
The design of this HTS could be optimized by choosing design parameters, such as the geometry design of elliptical tubes, heat transfer areas of heat exchangers, and other operating parameters. In this study, the optimization objective is chosen as minimizing the pumping power consumption P under a given heat transfer duty Q. Besides, heat transfer areas of components are given in advance. Decision variables of this optimization problem can be attributed to three levels: (1) operation parameters, i.e., operating frequencies of three VSPs, ωi (i = 1, 2, 3), (2) geometry design parameter, i.e., the semi-minor axis rb of the elliptical tube under a constant semi-major axis ra, and (3) local heat transfer process parameters, i.e., the flow and temperature fields in the elliptical tubes, U(x,y,z) and T(x,y,z). Figure 1. Sketch of a multi-loop heat transfer system [37]. The design of this HTS could be optimized by choosing design parameters, such as the geometry design of elliptical tubes, heat transfer areas of heat exchangers, and other operating parameters. In this study, the optimization objective is chosen as minimizing the pumping power consumption P under a given heat transfer duty Q. Besides, heat transfer areas of components are given in advance. Decision variables of this optimization problem can be attributed to three levels: (1) operation parameters, i.e., operating frequencies of three VSPs, ω i (i = 1, 2, 3), (2) geometry design parameter, i.e., the semi-minor axis r b of the elliptical tube under a constant semi-major axis r a , and (3) local heat transfer process parameters, i.e., the flow and temperature fields in the elliptical tubes, U(x,y,z) and T(x,y,z).

Heat Transfer and Fluid Flow Constraints of the HTS
The heat current model of the HTS shown in Figure 1 is given in Figure 2, where the heat flow direction in the HTS is presented. Applying the Kirchhoff's Voltage Law (KVL) on the heat current model yields the governing equation.
where Q is the heat transfer rate in the system. R 1 , R 2, and R 3 are inlet temperature difference-based thermal resistances of heat exchanger HX1, HX2 and the evaporator, respectively. Their expressions are [14]: where k i (i = 1, 2, 3) is the overall heat transfer coefficient of each heat exchanger. In the calculation, these parameters are initialized using estimated values and hence are marked by overline. Their values will be updated using simulation results. A i stands for the heat transfer area, and G i stands for the heat capacity rate of working fluids, i.e., the product of mass flow rate and specific heat (mc p ) i .  Figure 1 is given in Figure 2, where the heat flow direction in the HTS is presented. Applying the Kirchhoff's Voltage Law (KVL) on the heat current model yields the governing equation.
( ) where Q is the heat transfer rate in the system. R1, R2, and R3 are inlet temperature difference-based thermal resistances of heat exchanger HX1, HX2 and the evaporator, respectively. Their expressions are [14]: where i k (i = 1, 2, 3) is the overall heat transfer coefficient of each heat exchanger. In the calculation, these parameters are initialized using estimated values and hence are marked by overline. Their values will be updated using simulation results. Ai stands for the heat transfer area, and Gi stands for the heat capacity rate of working fluids, i.e., the product of mass flow rate and specific heat (mcp)i. The hydraulic constraints of the system are also required, which can be established by the balance relation of driving force and resistance of the fluid flow. The flow resistance in the system comes from fluid pressure drops in heat exchangers and pipelines. The total flow resistance in the ith loop is where Hi is the total head, Hs,i is the static head, Hd,i is the dynamic head of the pipeline network, and He,i is the head from pressure drops in heat exchangers. Static heads Hs,i are controlled by the structure of the pipeline network and can be treated as constants, while the dynamic head Hd,i is [19]: where bi are lumped parameters determined by the pipeline. The hydraulic constraints of the system are also required, which can be established by the balance relation of driving force and resistance of the fluid flow. The flow resistance in the system comes from fluid pressure drops in heat exchangers and pipelines. The total flow resistance in the i-th loop is where H i is the total head, H s,i is the static head, H d,i is the dynamic head of the pipeline network, and H e,i is the head from pressure drops in heat exchangers. Static heads H s,i are controlled by the structure of the pipeline network and can be treated as constants, while the dynamic head H d,i is [19]: where b i are lumped parameters determined by the pipeline.
where ρ i is the density of the working fluid, S i is the sectional area of the pipeline, L i stands for the length of the pipeline, F i refers to Darcy's coefficient, d i stands for the diameter of the pipeline, and K i is the local resistance coefficient [37]. Besides, heads related heat exchangers H e,i are determined by internal pressure drops: where ∆p i,h and ∆p i,c are pressure drops of both hot and cold fluids in heat exchangers. Like the overall heat transfer coefficients, here pressure drops are also initialized by estimated values to start the calculation and to be updated with simulation results. The governing equation of VSPs is required to determine the driving force of working fluids as follows: where a j,i (j, i = 0, 1, 2) are characteristic parameters of VSPs, and they remain constants in the optimization procedure. Therefore, the hydraulic constraint equation of the HTS can be obtained: The optimization objective, in this case, is the total pump power consumption of VSPs, which could be expressed as:

Optimization Equations of the HTS
The Lagrange multiplier method is applied to derive the optimization equations of the problem, and the Lagrange function is constructed as follows: where α and β i are Lagrange multipliers. Taking partial differentials of the Lagrange function with respect to each variable (m i and ω i ) and letting them equal to zero yields the Lagrange equations: Energies 2020, 13, 4112 6 of 19 Partial derivatives related to thermal resistances are provided in Appendix A. Solving these Lagrange equations together with heat current Equation (1) as well as the hydraulic constraint Equation (12) gives optimal results. However, the obtained solution is the optimal parameters with estimated overall heat transfer coefficients k i and pressure drops, ∆p i,h and ∆p i,c . Therefore, these optimization results are to be updated using fluid flow simulation results.

Simulation of Heat Exchangers HX1 and HX2
CFD simulation is used to determine overall heat transfer coefficients k i and pressure drops ∆p i,h and ∆p i,c under different operating conditions. Heat exchangers 1 and 2 are counter-flow plate exchangers with periodical structures, and hence they can be divided into heat transfer units consisting of two fluid channels and a heat conduction plate between them. Figure 3 3  2 3 3  2  3  3  3 3  3  3 3  2  3   3  2   3   3 3 Partial derivatives related to thermal resistances are provided in Appendix A. Solving these Lagrange equations together with heat current Equation (1) as well as the hydraulic constraint Equation (12) gives optimal results. However, the obtained solution is the optimal parameters with estimated overall heat transfer coefficients i k and pressure drops, these optimization results are to be updated using fluid flow simulation results.

Simulation of Heat Exchangers HX1 and HX2
CFD simulation is used to determine overall heat transfer coefficients ki and pressure drops Δpi,h and Δpi,c under different operating conditions. Heat exchangers 1 and 2 are counter-flow plate exchangers with periodical structures, and hence they can be divided into he   Boundary conditions for the simulation are given as follows. The entrance of the unit is treated with the velocity-inlet condition, where the velocities are set by mass flow rates m i : Energies 2020, 13, 4112 where A 1,inlet and A 2,inlet are inlet sectional areas of heat transfer units to calculate the average flow velocity. N u1 and n u2 are numbers of heat transfer units contained in two heat exchangers, respectively. The outlets adopt the standard pressure-outlet condition. The plate in the heat transfer unit is a copper plate with a thickness of 1.5 mm, and no-slip wall condition with conjugate heat conduction is assigned to the plate. Besides, the upper and lower surfaces are assigned to the periodic boundary condition, and the lateral sides and other zones contacted with the environment are treated with insulating surfaces. The gravity force, thermal radiation, and viscous heat generation are neglected in the simulation. The working fluids in three loops are water, which is regarded as incompressible, and its properties are treated as constant. Therefore, the governing equations are [38]: where U is the velocity, p stands for the pressure, and λ is the thermal conductivity. The simulation is implemented by ANSYS Fluent, and the Semi-Implicit Method for Pressure Linked Equations-Consistent (SIMPLEC) algorithm is used. The renormalization group (RNG) k-ε model is additionally applied when the flow is turbulent. The overall heat transfer coefficients k i and pressure drops ∆p i,h and ∆p i,c are extracted from simulation results when the computation reaches convergence. Specifically, the overall heat transfer coefficient is determined employing heat transfer rate q i , the heat transfer area A u,I and logarithm mean temperature difference (LMTD) ∆T lm,I of heat transfer units.

Convective Heat Transfer Optimization in the Evaporator
Both design parameters and convective heat transfer processes in the evaporator are to be optimized. The optimization method of the convective heat transfer process under a certain geometry condition is presented first. The flow in the evaporator remains laminar flow due to the small heat transfer duty Q and mass flow rate m 3 . The corresponding governing equations for the optimal flow field are [32,33]: where F is an additional volume force determined by: where C Φ is a parameter relating to the viscous dissipation rate: Here, C 0 , C 1 , and C 2 are Lagrange multipliers introduced, and C 0 is a constant [39], while C 1 and C 2 are field functions of temperature T, velocity U, and spatial coordinates (x, y, z). Equation (27) is the momentum equation of the working fluid with an additional volume force F. The value of C Φ represents the strength of volume force, and different values of C Φ lead to different longitudinal vortex structures. Equation (28) is the governing equation of Lagrange multiplier C 1 . Solving Equations (27) and (28) as well as the mass and energy conservation equations offers the optimal fluid velocity field with the largest overall heat transfer coefficient for a laminar forced convection process with a prescribed inlet velocity and a specific viscous dissipation rate [33].
The geometry and boundary conditions used in the evaporator optimization are as follows. The evaporator consists of 15 copper elliptical tubes with a length of 150 mm and a wall thickness of 0.7 mm. The structure sketch of the evaporator is presented in Figure 4. Besides, the semi-major axis of the elliptical tube r a remains 5 mm while the semi-minor axis of the elliptical tube r b varies between 2 mm to 5 mm. Since the length-diameter ratio of the tube is quite large, the elliptical tubes in the evaporator can also be divided into heat transfer units to simplify the numerical calculation. The heat transfer unit here is chosen as a segment of the elliptical tube with a length of 20 mm, and the computation grid for r a = 5 mm and r b = 2.5 mm is shown in Figure 5. The computation grid is the combination of the structured hexahedron grid in the boundary layer and the unstructured hexahedron grid in the inner domain to improve the calculation precision, and the grid number is around 116,000. Comparing with earlier works on the convective heat transfer optimization in tubes [31,32], the grid number is large enough to ensure the grid independency.
where CФ is a parameter relating to the viscous dissipation rate: Here, C0, C1, and C2 are Lagrange multipliers introduced, and C0 is a constant [39], while C1 and C2 are field functions of temperature T, velocity U, and spatial coordinates (x, y, z).
Equation (27) is the momentum equation of the working fluid with an additional volume force F. The value of CФ represents the strength of volume force, and different values of CФ lead to different longitudinal vortex structures. Equation (28) is the governing equation of Lagrange multiplier C1. Solving Equations (27) and (28) as well as the mass and energy conservation equations offers the optimal fluid velocity field with the largest overall heat transfer coefficient for a laminar forced convection process with a prescribed inlet velocity and a specific viscous dissipation rate [33].
The geometry and boundary conditions used in the evaporator optimization are as follows. The evaporator consists of 15 copper elliptical tubes with a length of 150 mm and a wall thickness of 0.7 mm. The structure sketch of the evaporator is presented in Figure 4. Besides, the semi-major axis of the elliptical tube ra remains 5 mm while the semi-minor axis of the elliptical tube rb varies between 2 mm to 5 mm. Since the length-diameter ratio of the tube is quite large, the elliptical tubes in the evaporator can also be divided into heat transfer units to simplify the numerical calculation. The heat transfer unit here is chosen as a segment of the elliptical tube with a length of 20 mm, and the computation grid for ra = 5 mm and rb = 2.5 mm is shown in Figure 5. The computation grid is the combination of the structured hexahedron grid in the boundary layer and the unstructured hexahedron grid in the inner domain to improve the calculation precision, and the grid number is around 116,000. Comparing with earlier works on the convective heat transfer optimization in tubes [31,32], the grid number is large enough to ensure the grid independency.   Boundary conditions for the optimization calculation are as follows. Since the refrigerant evaporates from the liquid to vapor, its temperature keeps at the boiling point Te. Besides, in the evaporator, the thermal resistance of heat conduction in the tube wall can be neglected due to the high conductivity of the tube wall. Therefore, the temperature of the tube wall is regarded as Te. That  Boundary conditions for the optimization calculation are as follows. Since the refrigerant evaporates from the liquid to vapor, its temperature keeps at the boiling point T e . Besides, in the evaporator, the thermal resistance of heat conduction in the tube wall can be neglected due to the high conductivity of the tube wall. Therefore, the temperature of the tube wall is regarded as T e . That is, the tube wall adopts the condition of the no-slip wall with a constant temperature. Moreover, periodical conditions are applied to the inlet and outlet of the heat transfer unit. Finally, the mass flow rate in each tube m u can be calculated by the total mass flow rate m 3 and the number of tubes in the evaporator n t : The heat transfer area of the heat transfer unit A u is given as: where L u is the length of the heat transfer unit of elliptical tubes (20 mm), D is the perimeter of the tube cross-section, and can be calculated by Ramanujan's approximation [40]: The total heat transfer area of the evaporator A 3 is then derived as: where L is the length of the tube (150 mm). The overall heat transfer coefficient of the evaporator, k 3 , is also determined by the LMTD approach. The discretization scheme and algorithms in the simulation are as follows: the SIMPLEC algorithm is applied for the coupling of pressure and velocity, the second order upwind scheme is adopted for the momentum and energy conservation equation, and the QUICK (Quadratic Upstream Interpolation for Convective Kinematics) scheme is used for the governing equation of Lagrange multiplier C 1 . Considering the robustness of the computation, the value of C Φ is chosen to be −4.8, which leads to a 4-vortex structure in the flow field [37,41].

Synergetic Optimization Procedure
The system-level optimization, heat exchanger simulation, and local heat transfer process level optimization in elliptical tubes of the evaporator have been presented, and they are required to be integrated by a synergetic optimization strategy. Besides, the local heat transfer optimization above is performed. Moreover, with fixed geometry configuration of elliptical tubes, while the design optimization of elliptical tubes is absent.
First, the heat transfer system with a given semi-minor axis of elliptical tubes r b is considered. The system-level optimization of the HTS is first performed with estimated parameters, i.e., estimated overall heat transfer coefficients k i and estimated pressure drops, ∆p i,h and ∆p i,c . Next, the CFD simulation is conducted for plate heat exchangers HX1, HX2. The convective heat transfer optimization in the evaporator is also performed using the numerical simulation. Overall heat transfer coefficients k i and pressure drops ∆p i,h and ∆p i,c derived from computation results generally differ from estimated values, and if deviations between actual parameters and estimated parameters, ε i , κ i,h , and κ i,c are out of the tolerance range, they will be then updated: where φ i , Ψ i,h , and Ψ i,c are relaxation factors, and they are empirically chosen to be 0.7, 0.5, and 0.5, respectively. After the update, the system optimization will be performed again with updated k i , ∆p i,h , and ∆p i,c . This iteration continues until deviations are within the tolerance range. Once the iteration converges, the optimal pump power consumption P with the current semi-minor axis r b can be derived from optimization results. The last step of the synergetic optimization is to search the optimum semi-minor axis r b within the given range to obtain the minimum pumping power consumption, where the downhill simplex (DS) algorithm is employed. When the iteration process of DS algorithm converges, the minimized pumping power consumption and corresponding optimal operation parameters including the frequencies of VSPs ω i , the semi-minor axis of elliptical tubes r b , and the optimal flow field U in the elliptical tubes are obtained. That is, the optimization of the elliptical tubes includes both component level and local heat transfer process level optimization. The component level optimization is to find the optimal semi-minor axis r b , and the local heat transfer process level optimization is to find the optimal flow field. Figure 6 gives a graphical summary of the proposed synergetic optimization method, and Figure 7 presents a detailed flowchart of the nested iteration process. Energies 2020, 13, x FOR PEER REVIEW 11 of 23

Optimization Results and Discussion
The main conditions used in the HTS optimization case are presented in Table 1. Besides, the geometry parameters of heat exchangers 1 and 2, characteristic parameters of the pipeline network, and VSPs remain the same in the optimization computation. The characteristic parameters of the pipeline network and VSPs are presented in Tables 2 and 3, which are identical with them in References [37,42] for comparison.

Optimization Results and Discussion
The main conditions used in the HTS optimization case are presented in Table 1. Besides, the geometry parameters of heat exchangers 1 and 2, characteristic parameters of the pipeline network, and VSPs remain the same in the optimization computation. The characteristic parameters of the pipeline network and VSPs are presented in Tables 2 and 3, which are identical with them in References [37,42] for comparison.    Table 4 compares the results of two cases, i.e., with and without component geometry optimization. Both two cases include system-level optimization and local convective heat transfer process optimization. In the case with component geometry optimization, the Reynolds number of the elliptical tube flow in the evaporator is 425.27, which verifies the working condition of laminar flow. The comparison shows that the optimized pumping power consumption P, decreases from 19.35 to 16.2 W with the semi-minor axis length decreases from 5.0 to 3.75 mm. Besides, the heat transfer area of the evaporator decreases from 0.071 to 0.0623 m 2 , and the operating frequencies of VSPs and the mass flow rates in loops also decrease. Therefore, the proposed synergetic optimization method combining the local heat transfer process, geometrical design, and system parameter optimization is more efficient for both energy and cost-saving. Table 4. Comparison of synergetic optimization results of the HTS.

Geometry Optimization
Included?  Figure 8 presents the sectional velocity ( Figure 8a) and the temperature field (Figure 8b) of the elliptical tubes after the optimization procedure converges. The results present that the optimal flow field for the elliptical tube has vertical vortexes to enhance the internal convective heat transfer. However, compared to the original flow field, the temperature field does not change significantly, which is attributed to the limited strength of volume force F. Besides, the convective heat transfer in the elliptical tube is indeed enhanced since the velocity field is changed, and the average field synergy angle of the entire field is increased [43]. Figure 9 presents the iteration process of searching the optimal semi-minor axis length r b within the given range, and iteration points are not distributed uniformly due to the application of the heuristic algorithm. Each point stands for an optimal total pumping power consumption P under a certain r b , and among them, the minimal total pump consumption point is the optimal one considering the geometry design optimization. The results read that the optimal semi-minor axis r b is 3.75 mm, and the optimal total pumping power consumption first decreases and then increases with an increasing r b . Besides, the optimal frequencies of VSPs ω i and mass flow rates m i in three loops have the same variation trend as shown in Figure 10a,b, respectively.  Figure 9 presents the iteration process of searching the optimal semi-minor axis length rb within the given range, and iteration points are not distributed uniformly due to the application of the heuristic algorithm. Each point stands for an optimal total pumping power consumption P under a certain rb, and among them, the minimal total pump consumption point is the optimal one considering the geometry design optimization. The results read that the optimal semi-minor axis rb is 3.75 mm, and the optimal total pumping power consumption first decreases and then increases with an increasing rb. Besides, the optimal frequencies of VSPs ωi and mass flow rates mi in three loops have the same variation trend as shown in Figure 10a and 10b, respectively.  The variations of the heat transfer area of the evaporator and overall heat transfer coefficients in the system are also investigated. The heat transfer area of evaporator A3, versus the semi-minor axis of elliptical tubes rb is presented in Figure 11. It is found that A3 increases linearly versus rb approximately. Therefore, it can be predicted that the overall heat transfer coefficient of the evaporator k3, will also be affected by the geometry optimization, and has a different variation trend comparing to the overall heat transfer coefficients of heat exchanger HX1 and HX2. The variations of the heat transfer area of the evaporator and overall heat transfer coefficients in the system are also investigated. The heat transfer area of evaporator A 3 , versus the semi-minor axis of elliptical tubes r b is presented in Figure 11. It is found that A 3 increases linearly versus r b approximately. Therefore, it can be predicted that the overall heat transfer coefficient of the evaporator k 3 , will also be affected by the geometry optimization, and has a different variation trend comparing to the overall heat transfer coefficients of heat exchanger HX1 and HX2.  Figure 12 shows three overall heat transfer coefficients k1, k2 and k3 under different rb, and the variation trends do match the prediction. Overall heat transfer coefficients of HX1 and HX2, i.e., k1 and k2, are merely determined by operation parameters via heat transfer correlations [34], since their geometry and structure parameters remain unchanged. Therefore, the variation trend of k1 and k2 corresponds to the trend of pump frequencies ω1 and ω2 as well as mass flow rates m1 and m2. However, the overall heat transfer coefficient of the evaporator k3, presents a different variation mode. When the semi-minor axis rb > 3.75 mm, the mass flow rate in the tube decreases with rb decreases, which weakens the heat transfer performance of the evaporator. On the other hand, the decrease of rb would strengthen the heat transfer in the tube comparing to round tubes. The competition of two factors together derives a slowly increasing overall heat transfer coefficient k3 with a decreasing rb. Meanwhile, when the semi-minor axis rb ≤ 3.75 mm, the mass flow rate m3 increases with a decreasing rb. In this case the two factors do not compete but result in a strengthened heat transfer performance, and the overall heat transfer coefficient k3 increases significantly with a decreasing rb.  Figure 12 shows three overall heat transfer coefficients k 1 , k 2 and k 3 under different r b , and the variation trends do match the prediction. Overall heat transfer coefficients of HX1 and HX2, i.e., k 1 and k 2 , are merely determined by operation parameters via heat transfer correlations [34], since their geometry and structure parameters remain unchanged. Therefore, the variation trend of k 1 and k 2 corresponds to the trend of pump frequencies ω 1 and ω 2 as well as mass flow rates m 1 and m 2 . However, the overall heat transfer coefficient of the evaporator k 3 , presents a different variation mode. When the semi-minor axis r b > 3.75 mm, the mass flow rate in the tube decreases with r b decreases, which weakens the heat transfer performance of the evaporator. On the other hand, the decrease of r b would strengthen the heat transfer in the tube comparing to round tubes. The competition of two factors together derives a slowly increasing overall heat transfer coefficient k 3 with a decreasing r b . Meanwhile, when the semi-minor axis r b ≤ 3.75 mm, the mass flow rate m 3 increases with a decreasing r b . In this case the two factors do not compete but result in a strengthened heat transfer performance, and the overall heat transfer coefficient k 3 increases significantly with a decreasing r b . Recall that in the physical picture presented by heat current method, heat flows from the hot end to cold end through thermal resistances, which describe heat transfer ability of components. In this system, the pressure drop in the evaporator is far smaller comparing to pressure drops in HX1 and HX2, and hence it only has a negligible effect on the system performance optimization. Therefore, the Recall that in the physical picture presented by heat current method, heat flows from the hot end to cold end through thermal resistances, which describe heat transfer ability of components. In this system, the pressure drop in the evaporator is far smaller comparing to pressure drops in HX1 and HX2, and hence it only has a negligible effect on the system performance optimization. Therefore, the variation trend of the system optimization result is mainly determined by the trade-off relation between the heat transfer area and the overall heat transfer coefficient of the evaporator. More specifically, the optimal result is generated from the trade-off between the heat transfer enhancement due to the increase of mass flow rate and due to the shape variation of tubes.
Although a rather simple multi-loop heat transfer system is adopted in this study as an example, the proposed synergetic optimization framework is not limited to this case. Recent studies have shown that heat current method can be used for heat transfer systems with different complicated topology structures and even for complex thermodynamic systems. Besides, component design optimization and local heat transfer process optimization also work for other heat exchangers. Therefore, the proposed optimization method here is not limited by specific topology structures of heat transfer systems but has a wide application scope.

Conclusions
Optimization of heat transfer systems benefits energy efficiency. In this study, a synergetic optimization method combining system-level optimization, component design optimization, and local heat transfer process optimization is proposed. The system-level optimization is implemented by the heat current method and hydraulic analysis, the component design optimization of the evaporator is performed using the downhill simplex algorithm, and the simulation of heat exchangers uses CFD computation. The local convective heat transfer process in the evaporator is optimized by the field synergy principle, where an additional volume force is introduced to obtain the optimal flow field. These three levels are integrated by using the proposed iterating and updating strategy to provide self-consistent and credible results such as heat transfer coefficients and pressure drops in the system.
A typical multi-loop HTS is optimized to minimize pumping power consumption using the proposed method as an example. The optimization results present that there exists an optimal value for the semi-minor axis r b for elliptical tubes in the evaporator, and the optimal value of r b is 3.75 mm. The optimal flow field in the elliptical tubes can be obtained by designing of the inner surface shape of tubes in practical applications. Besides, the synergetic optimization gives a 16.2 W total pump consumption, which is a 16.3% saving comparing the results without component design optimization. Therefore, the synergetic optimization including component design optimization is more efficient to optimize the system performance. The coupled relation between the system optimization, geometry design optimization, and local heat transfer process optimization can be properly handled by the iterating and updating process. On the other hand, optimization results under different semi-minor axis r b presents that the optimum r b is the result of two competing mechanisms of heat transfer performance enhancement. That is, the optimal r b is determined by the trade-off relation between heat transfer enhancement due to the mass flow rate increase and due to the shape variation of the tube. Finally, since the proposed method relies on the heat current method and CFD simulation as well as the iteration strategy, it is not limited in this pure heat transfer case and could be applied to various other thermal systems even integrated energy systems. Therefore, the proposed method has a wide application scope and would benefit the thermal system optimization.