Topology design of two-fluid heat exchange

Heat exchangers are devices that typically transfer heat between two fluids. The performance of a heat exchanger such as heat transfer rate and pressure loss strongly depends on the flow regime in the heat transfer system. In this paper, we present a density-based topology optimization method for a two-fluid heat exchange system, which achieves a maximum heat transfer rate under fixed pressure loss. We propose a representation model accounting for three states, i.e., two fluids and a solid wall between the two fluids, by using a single design variable field. The key aspect of the proposed model is that mixing of the two fluids can be essentially prevented without any penalty scheme. This is because the solid constantly exists between the two fluids due to the use of the single design variable field. We demonstrate the effectiveness of the proposed approach through three-dimensional numerical examples in which an optimized design is compared with a simple reference design, and the effects of design conditions (i.e., Reynolds number, Prandtl number, design domain size, and flow arrangements) are investigated.


Introduction
Heat exchangers are devices that transfer heat between two or more fluids. Recently, designing high-performance heat exchangers has been prioritized owing to an increased demand for such systems with low energy consumption. Most heat exchangers involve two-fluid heat exchange with indirect contact, i.e., there are two fluids separated by a wall. In such heat exchangers, the performances depend on the flow regime. For instance, the curving and branching of the flow channels of each fluid and the adjacency of the flow channels to each other have a significant effect on the heat transfer rate and the pressure loss. Therefore, consideration of such characteristics is significant in heat exchanger design.
Generally, in heat exchanger design, a type is first selected from the existing ones, and the details, e.g., size, shape, and fin type, are determined (Shah and Sekulic, 2003). Optimization methods for the existing types of heat exchangers have been well researched in previous works. For instance, Hilbert et al. (2006) studied a multiobjective optimization method concerning the blade shape of a tubular heat exchanger via a genetic algorithm. Kanaris et al. (2009) applied a response surface method to optimization problems for a plate heat exchanger with undulated surfaces. Guo et al. (2018) proposed a design method to simultaneously determine the fin types and their optimal sizes in plate fin heat exchangers.
The existing types of heat exchangers are typically manufactured through traditional processes such as pressing of plates and bending of tubes. This implies that the typical heat exchangers are designed under strict constraints pertaining to their manufacturability. However, there is an increasing demand to produce heat exchangers more efficient and compact, as the traditional manufacturing techniques are usually too inefficient to satisfy this demand. To produce innovative heat exchangers, these constraints pertaining to their manufacturability should be eliminated. In this regard, using additive manufacturing is a promising option, and has been employed for manufacturing innovative heat exchangers (e.g. Scheithauer et al., 2018) in instances where traditional techniques cannot be easily be applied. According to such progression of additive manufacturing technologies, revealing and understanding the ideal flow regime of the two fluids in heat exchange systems become increasingly important. However, the ideal flow regime is not obvious in most cases. For such a problem, topology optimization is a promising methodology, as it enables shape and topological changes of the structure.
Topology optimization was first proposed by Bendsøe and Kikuchi (1988). It essentially comprises introduction of a fixed design domain and replacement of the original optimization problem with a material distribution problem. A characteristic function is introduced to represent either solid material or void at the appropriate points in the design domain. Although topology optimization was originally applied to solid mechanics problems, it has been developed for application in various physical problems (Bendsøe and Sigmund, 2003;Sigmund and Maute, 2013;Deaton and Grandhi, 2014).
For fluid flow problems, Borrvall and Petersson (2003) proposed a topology optimization method for solid/fluid distribution problems in Stokes flow. The critical feature in their method is the introduction of a fictitious body force term, which can discriminate fluid and solid at points in the design domain, according to each design variable value. Their method was extended for Navier-Stokes equations (Gersborg-Hansen et al., 2005;Olesen et al., 2006), and design problems of fluidic devices, such as microreactors (Okkels and Bruus, 2007), micromixers (Andreasen et al., 2009;Deng et al., 2018), Tesla valves (Lin et al., 2015;Sato et al., 2017), and flow batteries (Yaji et al., 2018b;Chen et al., 2019). The thermal-fluid problems that have been studied include forced convection problems (Matsumori et al., 2013;Koga et al., 2013;Yaji et al., 2015), natural convection problems (Alexandersen et al., 2014;Coffin and Maute, 2016), and large-scale three-dimensional problems involving cluster computing (Alexandersen et al., 2016;Yaji et al., 2018a). In addition, Dilgen et al. (2018) applied density-based topology optimization to turbulent heat transfer problems. For practical approaches in utilizing topology optimization in complex fluid problems, simplified models have recently been studied for solving turbulent heat transfer problems (Zhao et al., 2018;Yaji et al., 2020;Kambampati and Kim, 2020), natural convection problems (Asmussen et al., 2019;Pollini et al., 2020), and heat exchanger design problems (Haertel and Nellis, 2017;Kobayashi et al., 2019). The history and the state-of-the-art of topology optimization for fluid problems including heat transfer problems are summarized in the recently published review paper .
Many studies regarding topology optimization of thermal-fluid problems dealt with solid/fluid distribution, i.e., only one fluid is considered. A topology optimization method for two-fluid and one-solid distribution problems is essential to generate an innovative design of a two-fluid heat exchanger. Tawk et al. (2019) has recently proposed a topology optimization method for heat transfer problems involving two-fluid and one-solid domains. In their method, two types of design variable fields are introduced. The first design variable field expresses either solid or fluid, and the second one expresses the ratio of the two fluids. The interpolation function in Borrvall and Petersson (2003) is extended for two-fluid problems using an interpolation function for multi-material topology optimization presented in the literature (Bendsøe and Sigmund, 1999). In addition, a penalty function is introduced to prevent the mixing of the two fluids, as an essential constraint. The proposed approach was applied to two-dimensional examples, and it was revealed that the optimized designs realized an increasing heat transfer rate and a decreasing pressure loss. The authors (Tawk et al., 2019), however, stated that the solution search becomes unstable in certain conditions because of the penalty scheme.
In the general framework of multi-material topology optimization, three states (two types of materials and void) are accounted for by using multiple design variable fields. Conversely, in topology optimization of two-fluid heat exchange, the two fluids should be separated by a wall. For two-fluid problems, solving the multi-material topology optimization typically requires a penalty scheme as with the previous work (Tawk et al., 2019), since the representation model for the multiple states has an excessive degree of freedom. There arises a possibility that a simpler representation model can be obtained by focusing on the characteristics of the two-fluid problems. For instance, de Souza and Silva (2020) proposed a representation model for design-dependent pressure load problems, in which three states (interior and exterior fluids, and a wall) are represented using a single design variable field. They demonstrated that the leakage between the two fluids can be circumvented by only using the single design variable field, whose middle value corresponds to the wall between the fluids.
In this study, we propose a representation model for topology optimization of two-fluid heat exchange. The novel feature of the proposed model is the representation of two different types of fluids and a wall by using a single design variable field, whereas, in the previous work (de Souza and Silva, 2020), the two fluids used were required to be of the same type. Therefore, we introduce two types of fictitious force terms in each flow field for representing a wall between the two types of fluid, which essentially cannot mix. The topology optimization problem is formulated as a maximization problem of the heat transfer rate under a fixed pressure loss. We emphasize that the two-fluid heat exchange problems should be treated as three-dimensional problems since the topological change is significantly restricted in the case of two-dimensional problems. Thus, we apply the proposed approach to three-dimensional heat exchange problems and examine the effects of several design conditions, such as Reynolds number, Prandtl number, design domain size, and flow arrangements. The remaining sections of this paper are organized as follows. In Section 2, the design problem of two-fluid heat exchange is presented, as well as the method to represent two fluids and one solid by the single design variable field is briefly discussed. In Section 3, the analysis model and the proposed representation model is presented, and the topology optimization problem is formulated. Furthermore, the numerical implementation is presented. In Section 4, several three-dimensional numerical examples are provided to validate the proposed method and to discuss the ideal flow regime of the two-fluid heat exchange. Finally, Section 5 concludes this study.
2 Topology optimization for two-fluid heat exchange 2.1 Design problem A simple heat exchanger, in which flow channels of the two fluids are separated by a wall as presented in Fig. 1, is considered. An inlet and an outlet of each fluid are respectively set, and each fluid flows independently. The heat is transferred between the two fluids in the design domain D. The hot and cold fluids are denoted by "fluid 1" and "fluid 2," respectively.
In this study, we consider a topology optimization problem for designing the shape and topology of the flow channels and the separating wall in the design domain D. The objectives of the heat exchanger design are the maximization of the heat transfer rate and minimization of the pressure loss. This multi-objective design problem is transformed into a single-objective design problem by fixing the pressure loss. The detailed formulation is provided in Section 3.

Basic concept of topology optimization
Topology optimization essentially involves the replacement of the original structural optimization problem with a material distribution problem in a fixed design domain D that contains the original design domain Ω. Bendsøe and Kikuchi (1988) introduced a characteristic function χ Ω that assumes a discrete value 0 and 1, as follows: where x represents the position in the fixed design domain D. Since the characteristic function χ Ω is a discontinuous function, a relaxation technique is required for gradient-based optimization. In the density-based method (Bendsøe and Sigmund, 2003), the characteristic function χ Ω is replaced with a continuous scalar function 0 ≤ ψ(x) ≤ 1, which corresponds to the raw design variable field in the topology optimization problem.
Solid Mixture of solid and fluid 1 Fluid 1 0 < γ 2 < 1 Solid Mixture of solid, fluid 1, and fluid 2 Mixture of fluid 1 and fluid 2 γ 2 = 1 Solid Mixture of solid and fluid 2 Fluid 2 Table 2: Representation model of two fluids and one solid in the proposed method To ensure the spatial smoothness of the material distribution in D, the filtering technique is generally applied in topology optimization. In this study, Helmholtz-type filter (Kawamoto et al., 2011;Lazarov and Sigmund, 2011) is employed, as follows: where γ is a smoothed design variable field, and R is a filtering radius. By solving (2), ψ is mapped to γ, thus ensuring a spatial smoothness in D.

Representation model
For the topology optimization of two-fluid heat exchange, the design variable field should represent the two fluids and the one solid, as shown in Fig. 1. Here, we discuss the method to represent the two fluids and the one solid by the design variable field. In the previous work regarding topology optimization of two fluids (Tawk et al., 2019), the method presented was based on multi-material topology optimization. The representation of the two fluids by design variable fields are presented in Table 1. Here, design variable 1, corresponding to γ 1 , represents the ratio of the solid and the fluids, whereas design variable 2, corresponding to γ 2 , represents the ratio of fluid 1 and 2.
Multi-material topology optimization is extensively studied regarding distribution problems of materials and void (Bendsøe and Sigmund, 1999). In the typical multi-material topology optimization, two types of design variable fields are utilized to represent two types of materials and void. In this regard, the distribution problem of two fluids and one solid is similar in that the three states are accounted for. However, it differs in the adjacency of the three states. In the two-fluid heat exchange problems, fluids 1 and 2 should not be adjacent and mixed. As revealed in Table 1, fluid 1 and 2 can be adjacent and mixed when γ 1 = 1 and 0 < γ 2 < 1. Since the penalty scheme is required to prevent mixing and adjoining of the two fluids, the method based on the multi-material topology optimization is complex and may lead to an unstable solution search.
The state representation should be simple, and the penalty scheme should not be required for a stable solution search. In this study, we propose a new method that uses a single design variable field. Table 2 demonstrates the state represented by the design variable field. The design variable field γ assumes 1 in fluid 1 and 0 in fluid 2. The significant difference to the previous work (Tawk et al., 2019) is to represent the solid in the intermediate value of γ. This is because the two fluids are separated by the wall in the heat exchange problem. The two fluids can consistently be separated by the wall when the spatial smoothness of the design variable field is preserved in the design domain D. It should be noted that a filtering technique (2), is essential in the proposed representation model. In this setting, the mixing of fluids 1 and 2 can be eliminated without using the penalty function. Thus, the main feature of this method is that the two fluids and the one solid can be represented simply by the single design variable field without any penalty schemes.

Governing equations and boundary conditions
The governing equations of an incompressible steady flow are given as follows: where u is the velocity, ρ is the density, p is the pressure, µ is the viscosity, and F is a body force. In the fluid topology optimization, F is a fictitious body force caused by the solid objects in the flow (Borrvall and Petersson, 2003). The energy equation without the heat source is defined as follows: where C p is the specific heat at constant pressure, T is the temperature, and k is the thermal conductivity.
The boundary conditions are defined as follows: where subscripts 1 and 2 represent the fluid 1 and fluid 2, respectively. As mentioned in Section 2.1, the constant pressure loss is assumed by fixing the pressure at the inlet and outlet.

Proposed interpolation scheme
In this study, the different fictitious forces in the governing equations for fluid 1 and fluid 2 are imposed separately to realize the state represented by the design variable field indicated in Table 2. Therefore, the governing equations are separated for fluid 1 and fluid 2, as follows: where subscripts 1 and 2 represent fluid 1 and fluid 2, respectively. The last terms in (12) and (14) are crucial in the proposed method. These fictitious body forces are defined as follows: where α 1 (x) and α 2 (x) are the inverse permeability of the porous media at position x (Borrvall and Petersson, 2003). The fluid hardly flows when α 1 or α 2 assumes a large value. Consequently, the fluid can flow freely when α 1 or α 2 assumes a zero value. The inverse permeabilities α 1 and α 2 are defined as follows: where α max is the maximum inverse permeability and q is the parameter that controls the convexity of α. Both inverse permeabilities are defined based on the interpolation function by Borrvall and Petersson (2003). However, these two  Figure 2 displays the graph of the interpolation functions at α max = 10 4 and q = 0.01. As revealed in Fig. 2, if γ = 0, α 1 assumes a large value α max and α 2 assumes a zero value. For this reason, only fluid 2 can flow whereas fluid 1 cannot flow. Conversely, if γ = 1, only fluid 1 can flow since α 1 assumes a zero value and α 2 assumes α max . In the intermediate value of γ, both fluids hardly flow because both α 1 and α 2 are not zero. Therefore, the flow characteristics of Table 2 can be achieved via these interpolation functions of the inverse permeability.
To simplify the numerical implementation, the governing equations are replaced with dimensionless form, and the dimensionless quantities, i.e., the Reynolds number and the Péclet number, are introduced.
The dimensionless equations of (11)-(14) are defined as follows: where the asterisk symbol represents dimensionless variables defined by a characteristic length L, characteristic speeds U 1 and U 2 as follows: The characteristic speeds U 1 and U 2 cannot be defined using the magnitude of the inlet velocity because the inlet pressure is fixed, as mentioned in Section 3.1. Therefore, U 1 and U 2 are defined using the pressure loss (Yaji et al., 2015), as follows:  As a result, the proposed interpolation functions α * 1 (γ) and α * 2 (γ) in the dimensionless form are given as: where α * max is the maximum inverse permeability. The dimensionless equation of (5) is defined as follows: Here, Pe is the Péclet number Pe = PrRe, in which Pr is the Prandtl number. The dimensionless temperature T * and the Prandtl number Pr are respectively defined as: In (29), since the dimensional parameters of thermal-fluid properties are integrated into the dimensionless parameter, Pe, an interpolation function should be introduced so that Pe depends on the states, i.e., fluid 1, fluid 2 or the solid. Therefore, we propose an interpolation function based on the Gaussian function, as follows: where s is the parameter that controls the shape of the Gaussian function. Pe f1 , Pe f2 , and Pe s are the Péclet number of fluid 1, fluid 2, and the solid, respectively. As shown in Fig. 3, Pr in the intermediate region corresponding to the solid can be smaller than that of the two fluids. This indicates that the solid has larger thermal conductivity than that of the two fluids, under the condition that the remaining parameters of thermal-fluid properties are fixed.

Optimization formulation
As discussed in Section 2.1, the objective of the design problem is to maximize the total heat transfer with a fixed pressure loss. The total heat transfer is given as follows: whereṁ is the mass flow, and T is the mean temperature. Introducing the boundary integral form, the equation for total heat transfer can be rewritten as: Since incompressible flow is assumed in this study, the density is constant and the volume flows at the inlet and outlet are equal. In addition, the inlet temperature T is constant as defined in (6) and (8). Thus, the equation can be further rewritten as: To replace the objective function with dimensionless form, a total heat transfer Q is nondimensionalized as follows: Based on (35), we define the objective function in this study as the sum of the total heat transfer in fluid 1 and fluid 2, as follows: For brevity, the asterisk symbol of the dimensionless variables is dropped henceforth. Finally, the topology optimization problem is formulated as follows: where the thermal-fluid field is governed by (19)-(22), and (29).

Numerical implementation
The numerical solution of the governing equations (19)- (22) and (29) is obtained using COMSOL Multiphysics 5.4. The analysis domain is discretized using tetrahedral finite elements. The scalar function ψ(x) is transformed into a discrete form ψ i (i = 1, . . . , n), where n is the total number of nodal points in the design domain D. The initial values of ψ i are set to 0.5 uniformly. For the sensitivity analysis, we employ the adjoint method based on the discrete formulation described in the literature (Bendsøe and Sigmund, 2003). The optimization procedure in this study is described as follows: Step 1. The initial values of ψ i are set in D.
Step 2. The smoothed design variables γ i are derived using the Helmholtz-type filter in (2).
Step 4. If the value of J is converged, the optimization process ends. Otherwise, the sensitivity is calculated using the adjoint method.
Step 5. The design variables are updated using sequential linear programming (SLP), and the optimization procedure then returns to Step 2.
4 Numerical examples 4.1 Optimized flow fields Figure 4 shows the three-dimensional analysis model and its dimensions. The dimensions L x , L y , and Lz are set to 4L, 2L, and 4L, respectively. A symmetrical boundary condition is imposed in the ZX plane. The analysis domain is discretized using 5.3 × 10 5 tetrahedral elements. The inlets and outlets of the two fluids are set as shown in Fig. 5, which is referred to as counter flow. The dimensionless quantities are set to Re 1 = 100, Re 2 = 100, Pr f1 = 7, Pr f2 = 7, Pr s = 3.5. The parameters regarding the interpolation functions are α max = 1 × 10 4 , q = 0.01, s = 0.1. The filtering radius R is set to L/12. Figure 6 demonstrates the optimized result in which the isosurface of γ = 0.5 and the velocity vectors are depicted. The separating wall is colored with orange on fluid 1 side, and light blue on fluid 2 side. The color of the arrows represents the temperature of the fluids. The optimized result of the half sectional view is shown in Fig. 7. A validation of the optimized result is presented in the Appendix, where the optimization model is compared with the high-fidelity model that uses a body-fitted mesh.
The optimization history of γ is illustrated in Fig. 8, and the convergence history of J is shown in Fig. 9. From  Fig. 8, fluid 1 and fluid 2 domains appear as the optimization proceeds. Since the intermediate value of γ exists between fluid 1 and fluid 2, the separating wall can be represented. The objective function value converges by approximately the 100th step, as depicted in Fig. 9.
Geometrically complex flow channels are obtained, as shown in Figs. 6 and 7. The flow channels are gradually curved and branched to circulate the fluids in the whole heat exchange domain, and fluid 1 and fluid 2 are placed alternately. This is because increasing the heat transfer area without the sharp bend is effective for improving the total heat transfer with low pressure loss. In addition, it is found that the outlet temperature of fluid 1 is lower than that of fluid 2. As illustrated in Fig. 7, the downstream region of fluid 1 adjoins the upstream region of fluid 2 and vice versa. Thus, the fluids can exchange heat in the downstream region, whereas the heat is already transferred in the upstream branched region. This feature is unique to the counter flow and can enhance the heat transfer efficiency.

Performance comparison with reference design
We now verify the performance of the optimized design shown in Fig. 6, by comparison with the simple reference design that is composed of straight flow channels. Figure 10 illustrates the reference design and its velocity vector. It should be noted that the scale of the vector is not the same as that of the optimized design.    Table 3 shows the comparison of the total heat transfer J, the flow rate of fluid 1, i.e.,V 1 , and the mean outlet temperature of fluid 1, i.e.,T out,1 . The flow rate of the reference design is over ten times higher than that of the optimized design. However, the mean outlet temperature is approximately the same as the inlet temperature in the reference design. Consequently, the mean outlet temperature is low in the optimized design. Thus, the total heat transfer is over four times higher than that of the reference design. Therefore, it is observed that the optimized design enhances heat transfer under the investigated operating condition.

Effects of design conditions 4.3.1 Effect of Reynolds number
Here, we investigate the effect of Reynolds number, Re, settings on the optimized configuration. In this study, the Reynolds number represents the pressure loss, because the characteristic speed is based on the pressure loss as defined in (25). Figure   From Figs. 7 and 11, the number of branches increases as the Reynolds number increases. In Fig. 11, a relatively simple structure is obtained in low Re. The heat transfer surface is approximately planar instead of the circular channels. Conversely, narrow branched channels are obtained for high Re. This is because the heat transfer area can be increased when the high pressure loss is allowed. In Fig. 11(b), it can be observed that the outlet temperature of fluid 1 is significantly low for high Re, owing to the complex flow channels. These results indicate that the ideal topology of the two fluids varies by the Reynolds number setting. Figure 12 shows the optimized results in a case of Pr f1 = 14, Pr f2 = 14. The parameters except for the Prandtl number are similar to the case of Pr f1 = 7, Pr f2 = 7 presented in Section 4.1. The Prandtl number shown in (30) is defined as the ratio of momentum diffusivity to thermal diffusivity. Hence, high Prandtl number results in low total heat transfer.

Effect of Prandtl number
In Fig. 12, complex flow channels are obtained, which is similar to the case of Re 1 = 200, Re 2 = 200 shown in Fig. 11(b). The reason for this similarity is related to the trade-off between the total heat transfer and the pressure loss. Increasing the Reynolds number corresponds with the eased pressure loss limit. In contrast, increasing the Prandtl number corresponds with the strict requirement for total heat transfer. This can be explained by the fact that Péclet number of both cases in Fig. 11(b) and Fig. 12 are identical; Pe f1 = 1400, Pe f2 = 1400.

Effect of design domain size
To investigate the effect of the design domain size, the optimization is performed in various L x , L y , and L z . Figure 13 shows the optimized results for L x = 2L, L y = 2L, and L z = 4L, as well as for L x = 4L, L y = 4L, and L z = 4L. The parameters are similar to the case of L x = 4L, L y = 2L, and L z = 4L presented in Section 4.1.
Although the flow channels are placed alternately regardless of the design domain size, different features are obtained in each case. For Fig. 13(a), the branched channels are in the z direction instead of the x direction as found in the cases of Figs. 7 and 13(b). Since the length in x direction is short and the heat exchange domain is small, a relatively sharp bend is obtained. From Fig. 13(b), thick and flat flow channels are realized when the length in y direction is large. This is because the narrow flow channels are not necessary in case of a large heat exchange domain.

Effect of flow arrangement
Finally, the effect of the flow arrangement, i.e., the arrangement of the inlet and the outlet in each fluid, is investigated. Figure 14 displays the schematic diagram and the optimized result. Therefore, we denote the arrangement of (b) and (c) in Fig. 14 by "U-counter" and "U-parallel", respectively. From Fig. 14(a), the flow channels of each fluid are placed alternately at right angles in parallel flow. It is discovered that the heat is mainly exchanged in the upstream region, and the temperature is nearly unchanged in the downstream region. This is because both fluids have approximately the same temperature after passing the branched section, and downstream fluid 1 cannot be adjacent to the upstream fluid 2, and vice versa, unlike in the case of counter flow. Similar features are found in Fig. 14(c), U-parallel flow. In Fig. 14(b), the high temperature difference can be maintained in the whole heat exchange surface, which is similar to counter flow. Based on these results, it is revealed that the outlet of the fluid 1 should be close to the inlet of the fluid 2, and vice versa, for efficient heat exchange.

Conclusions
In this paper, we proposed a density-based topology optimization method for two-fluid heat exchange problems. The novel aspect of the paper is that a representation model involving three states (two types of fluids and a solid wall) is defined using a single design variable field without any penalty schemes. To this end, we introduced two types of fictitious force terms in each flow field so that the two fluids cannot essentially mix.
We formulated the optimization problem as a maximization problem of the heat transfer rate under the fixed pressure loss. Based on the formulation, we investigated the effectiveness of the proposed approach through threedimensional numerical examples.
In the numerical examples, we first showed the performance of the optimized design in comparison with a simple reference design composed of straight flow channels. It was found that the optimized design achieves a heat transfer performance over four times higher than that of the reference design. We demonstrated the dependency of the optimized design with respect to the design conditions, i.e., Reynolds number, Prandtl number, design domain size, and flow arrangements. Accordingly, it was found that the proposed approach can generate effective optimized designs for optimization problem concerning the two-fluid heat exchange system, at least under the operating conditions investigated in this paper.

Replication of results
The necessary information for replication of the results is presented in this paper. The interested reader may contact the corresponding author for further implementation details.