Topology Optimization of Turbulent Flow Cooling Structures Based on the k-ε Model

Topology optimization (TO) is an effective approach to designing novel and efficient heat transfer devices. However, the TO of conjugate heat transfer has been essentially limited to laminar flow conditions only. The present study proposes a framework for TO involving turbulent conjugate heat transfer based on the variable density method. Different from the commonly used and oversimplified Darcy model, this approach is based on the more accurate and widely accepted k-ε model to optimize turbulent flow channels. We add penalty terms to the Navier–Stokes equation, turbulent kinetic energy equation, and turbulent energy dissipation equation, and use interpolation models for the thermal properties of materials. A multi-objective optimization function, aiming to minimize the pressure drop and the average temperature, is set up to balance the thermal and hydraulic performance. A case study is conducted to compare various optimization methods in the turbulent regime, and the results show that the present method has substantially higher optimization effectiveness while remaining computationally inexpensive.


Introduction
Many engineered products are subject to various degrees of thermal loads during operation and often must be actively cooled with a fluid.For many of these products, their performance relies heavily on cooling effectiveness.For example, gas turbine vanes and blades draw cold air to flow through their internal channels, thus allowing higher working temperatures and boosting power and efficiency [1,2].Electronic devices continue to be miniaturized and are generating increasingly more waste heat that is difficult to remove, which has become a bottleneck for future device development [3][4][5].Lithium-ion batteries also produce large amounts of heat during very fast charging and accidental thermal runaways and require more effective liquid cooling than is currently available [6][7][8].
The design of the cooling channels has a significant impact on the cooling performance, and the goal of cooling channel design is typically to minimize the flow resistance while maximizing the convective heat transfer rate [9][10][11][12].Conventional fluid channel design practice is largely based on intuition and a trial-and-error process, which is both inefficient and ineffective in generating the best structure.Topology optimization as an algorithmbased design method has been successfully applied to the optimization of mechanical structures since it was proposed by Bendsøe and Kikuchi [13].Popular topology optimization methods include homogenization methods [13], evolutionary structural optimization (ESO) methods [14], variable density methods [15], level set methods [16], moving morphable components (MMC) methods [17], etc.As for topology optimization of fluid flow, Borrvall and Petersson [18] pioneered a method of adding a penalty term to the Stokes equation to penalize the fluid velocity and realize the optimal design of flow channels.Among the following studies, A. Gersborg-Hansen et al. [19] extended the fluid topology optimization to Navier-Stokes flow with low Reynolds numbers based on Borrvall's research.Pietropaoli M et al. [20] proposed a design method based on topology optimization for the internal cooling channels of the gas turbine.Han et al. [21] used topology optimization to design spider-web-shaped heat sinks and experimented with different objective functions.The resulting structure has a 57.35% lower spatial temperature difference (T max −T min ) and 8% lower thermal resistance compared with the traditional design.Chen et al. [22] designed a topology-optimized cold plate for a lithium-ion battery thermal management system and achieved a higher heat transfer coefficient and significantly lower spatial temperature difference.Hu et al. [23] applied topology optimization to design microchannel heat sinks and demonstrated superior cooling performance compared with a conventional design of straight fin microchannel.
Turbulent flow is ubiquitous in advanced cooling structures, such as the active cooling of gas turbine blades, the regenerative cooling of rocket engines, and even the liquid cooling of high-power density electronics.However, topology optimization of conjugate heat transfer involving turbulent flow has been challenging because of the combined strong nonlinearity of turbulence, thermal-fluid coupling, and topology optimization itself.For pure turbulent fluid flow problems, a number of topology optimization schemes have been developed.Yoon [24] proposed a topology optimization method based on the Spalart-Allmaras (SA) turbulence model using the finite element method.Luís F.N. Sáet al. [25] extended the Spalart-Allmaras turbulence topology optimization method to flow channel design in a rotating domain.Yoon [26] developed a new topology optimization method based on the k-ε turbulent model by adding penalty terms to the Navier-Stokes equation, the turbulent kinetic energy equation, and the turbulent energy dissipation equation, respectively.Dilgen et al. [27] proposed a topology optimization method applicable to the k-ω turbulence model and the topology optimization flow channel design significantly reduced the vortex.For conjugate heat transfer design problems in the turbulent regime, the development of topology optimization methods is still at a much earlier stage.The majority of existing studies adopt the Darcy flow model to achieve stable numerical convergence and reduce computation time.Zhao et al. [28] used a low-cost Darcy flow model as an approximate substitute for turbulence flow in a topology optimization problem considering conjugate heat transfer.An-Li [29] studied the method of permeability setting and pressure drop constraint setting based on the Darcy flow topology optimization.In the topology optimization of the two-dimensional flow channel, Li et al. [30] used a simplified Darcy flow model to approximate turbulent flow.They applied the convective heat transfer model and Darcy flow model to solve the microchannel region's temperature and pressure fields, and the microchannel layout in a heat sink is designed by topology optimization.While the computational cost is low, these Darcy-based methods do not capture the essence of turbulence due to a lack of inertia terms, and therefore greatly sacrifice the optimization quality.Dilgen et al. [31] presented a novel topology optimization method based on the k-ω turbulence model.However, the objective function did not consider flow resistance, which must be taken into account in typical cooling designs.In addition, the k-ω model is mainly applied to fluid flow at a low Reynolds number and shear flow, and the applicability is limited.
The k-ε turbulent RANS model is one of the most effective and widely used turbulent flow models in engineering.However, an effective approach is lacking currently to apply the k-ε model in topology optimization for conjugate heat transfer.In this study, a conjugate heat transfer topology optimization method based on the k-ε turbulent model is developed to design cooling channels.The variable density method is used for optimization.The novelty is mainly reflected in the penalty terms which are added to the momentum equation, turbulent kinetic energy equation, and dissipation rate equation, respectively, and the material's thermal properties are interpolated to distinguish the solid and fluid phases.Then, the turbulent conjugate heat transfer topology optimization equations are obtained.The minimum inlet and outlet pressure drop and the minimum average temperature are set as multi-objective functions to meet the requirement of hydraulic and thermal performance.The turbulent conjugate heat transfer topology optimization structure is developed by carrying out the sensitivity analysis with the MMA algorithm in the topology optimization process.The optimization results using the k-ε model are compared with those using the Darcy flow model as well as the conventional straight flow channel, and the topology optimization method based on the k-ε turbulent model showed the lowest flow resistance and average temperature in the design domain.

Model of the Topology Optimization Problem
In this paper, the k-ε turbulent model and the conjugate heat transfer model are used to realize the topology optimization design of cooling channels.The turbulent conjugate heat transfer topology optimization conceptual design model is shown in Figure 1.The physical design domain consists of the solid phase domain Ω s and the fluid domain Ω f .Materials in the design domain are regarded as porous media, and the design variable γ, which affects the material's physical properties, varies from zero to unity, where γ = 1 represents fluid and γ = 0 represents solid.The solid isotropic material with penalization (SIMP) interpolation model is used to speed up the conversion of porous media to solids or fluids during the iterative process.A uniform heat input Q is applied to the entire design domain and different boundary conditions are set, including inlet Γ in , outlet Γ out , and wall boundaries Γ wall according to the actual flow and heat transfer conditions.
Entropy 2023, 25, x FOR PEER REVIEW 3 of 20 phases.Then, the turbulent conjugate heat transfer topology optimization equations are obtained.The minimum inlet and outlet pressure drop and the minimum average temperature are set as multi-objective functions to meet the requirement of hydraulic and thermal performance.The turbulent conjugate heat transfer topology optimization structure is developed by carrying out the sensitivity analysis with the MMA algorithm in the topology optimization process.The optimization results using the k-ε model are compared with those using the Darcy flow model as well as the conventional straight flow channel, and the topology optimization method based on the k-ε turbulent model showed the lowest flow resistance and average temperature in the design domain.

Model of the Topology Optimization Problem
In this paper, the k-ε turbulent model and the conjugate heat transfer model are used to realize the topology optimization design of cooling channels.The turbulent conjugate heat transfer topology optimization conceptual design model is shown in Figure 1.The physical design domain consists of the solid phase domain Ωs and the fluid domain Ωf.Materials in the design domain are regarded as porous media, and the design variable γ, which affects the material's physical properties, varies from zero to unity, where γ = 1 represents fluid and γ = 0 represents solid.The solid isotropic material with penalization (SIMP) interpolation model is used to speed up the conversion of porous media to solids or fluids during the iterative process.A uniform heat input Q is applied to the entire design domain and different boundary conditions are set, including inlet Γin, outlet Γout, and wall boundaries Γwall according to the actual flow and heat transfer conditions.

Governing Equations
The pseudo-variable γ is introduced into the flow and heat transfer control equations based on the variable density method.The flow and heat transfer equations are solved to compute the temperature, velocity, and pressure for the entire design domain first, and then the pseudo-variables are iteratively updated through sensitivity analysis and optimization solutions.Subsequently, the updated pseudo-variables are used to optimize the temperature and velocity-pressure distributions in the design domain, which enables the coupling of topology optimization and the turbulent conjugate heat transfer model.

Turbulent Flow Model
The widely used two-equation k-ε RANS model is directly implemented in the topology optimization process.The Reynolds stress tensor is regarded as a function of the turbulent viscosity in the standard k-ε turbulent model.The turbulent viscosity is solved Topology optimization conceptual design model.

Governing Equations
The pseudo-variable γ is introduced into the flow and heat transfer control equations based on the variable density method.The flow and heat transfer equations are solved to compute the temperature, velocity, and pressure for the entire design domain first, and then the pseudo-variables are iteratively updated through sensitivity analysis and optimization solutions.Subsequently, the updated pseudo-variables are used to optimize the temperature and velocity-pressure distributions in the design domain, which enables the coupling of topology optimization and the turbulent conjugate heat transfer model.

Turbulent Flow Model
The widely used two-equation k-ε RANS model is directly implemented in the topology optimization process.The Reynolds stress tensor is regarded as a function of the turbulent viscosity in the standard k-ε turbulent model.The turbulent viscosity is solved by the turbulent kinetic energy equation (k-equation) and turbulent energy dissipation equation (ε-equation).
In laminar flow topology optimization, a penalty term f related to the flow velocity is added to the momentum equation to distinguish solid and fluid regions in the design domain.In the fluid region, the penalty term is approximately equal to 0, while in the solid region, the penalty term is infinitely large, and the flow velocity tends to reach 0. The penalty term f is obtained by interpolation of the design variables.This method is also known as the Brinkman penalty model, and it is currently the most common method for solving laminar flow topology optimization.
As for solving topology optimization involving turbulent flow using the k-ε turbulence model, in addition to the penalty term in the momentum equation, the influence of the turbulent kinetic energy and the turbulent energy dissipation rate on the topology optimization design needs to be considered.Therefore, additional penalty terms f k and f ε are added to the turbulent kinetic energy equation and turbulent energy dissipation equation, respectively.The fluid has turbulent kinetic energy and energy dissipation, so f k and f ε tend to be zero in the fluid region.Whereas there is no turbulent kinetic energy and energy dissipation in the solid region, so f k and f ε tend to be infinity in solid region.
For incompressible steady-state fluid, the conservation equations for turbulent topology optimization are given as follows: The modified (added penalty terms f k and f ε ) turbulent kinetic energy equation and energy dissipation equation are expressed as: where In the turbulent kinetic energy equation (Equation (3)) and the energy dissipation equation (Equation (4)), µ T is the turbulent viscosity and P k is the production term.The parameters in Equations ( 3) and (4) recommended by the most literature are σ k = 1, C 1ε = 1.44,C 2ε = 1.92, σ ε = 1.3, C µ = 0.09.Penalty terms f , f k , and f ε are applied to the momentum equation, the turbulent kinetic energy equation, and the energy dissipation equation, respectively.α f , α k , α ε are the local inverse permeability of the porous medium, which can be calculated by the following interpolation functions.
where γ is the design variable and p α is the resistance coefficient penalty factor which is set to 0.01.α f min , α kmin , and α εmin are the inverse permeability in the fluid region and are equal to 0 in this paper.α f max , α kmax , and α εmax are the inverse permeability in the solid region and they should be given as infinity theoretically.However, limited by the computer's storage space, the values of α f max , α kmax , α εmax are set as 10 7 , 10 8 , and 10 8 , respectively.

Conjugate Heat Transfer Model
The SIMP interpolation model is adopted for the thermal properties of materials.
where k, ρ, c are the thermal conductivity, density, and specific heat capacity of the materials.
The physical parameters of the solid and fluid materials are shown in Table 1.In order to achieve the 0-1 structure, it is necessary to introduce a penalty factor p to penalize the intermediate variables so as to suppress the influence of the intermediate variables on the results.The value of the penalty factor needs to be selected according to actual problems and numerical results.In this paper, the penalty factors p 1 , p 2 , p 3 are set to 3, 3, 1.
The conjugate heat transfer equation can be given as Table 1.The values of fluid and solid physical parameters used in this study [21].

Topology Optimization Model
When designing cooling channels, both hydraulic and thermal performance need to be considered, which can be quantitatively evaluated by the pressure drop and average temperature, respectively.In this paper, the optimization objectives are set to reduce the pressure drop of the fluid from the inlet to the outlet as well as lower the average temperature in the design domain.The pressure drop and average temperature can be calculated by the following equations, respectively.
where Γ 2 and Γ 1 are the inlet and outlet boundaries, and the Ω is the design domain.
The hydraulic and thermal objectives are normalized to avoid numerical instabilities in the solution process.Then, objective functions are weighted and added to obtain the multi-objective function, which is given as follows: where w i and w j are the weighting factors of the multi-objective function.ψ 0 and φ 0 are the initial values of pressure drop and average temperature, respectively.
The topology optimization problem can be described as where v f is the volume fraction of fluid and is set to 0.5.γ i is the design variable in each domain element.The temperature field, pressure field, and velocity field are discretized using Lagrange linear interpolation.The GCMMA algorithm is used for topology optimization [23].

Density Filtering
Checkerboard is a common problem in topology optimization results, which manifests as a periodic distribution of solid and fluid elements and is thought to be caused by poor numerical modeling [32].The appearance of the checkerboard is usually inhibited by density filtering.Moreover, it should be noted that for the conjugate heat transfer problem, density filtering can reduce the ill-posedness of the topology optimization results [21].A density filter based on the Helmholtz equation is used to avoid checkerboard and meshdependency problems in topology optimization [33] and is defined as follows: where, γ is the filtered design variable and r is the filter radius.

Numerical Example
The numerical example of turbulent conjugate heat transfer topology optimization is shown in Figure 2. Due to the symmetry of the domain, the upper half of the region is selected as the design domain and a symmetrical boundary condition is imposed on the symmetrical line.The characteristic length of the model L is 0.02 m.A Reynolds number of 10,000 (Re = ρ in U in L µ in = 1000×0.5×0.020.001 = 10, 000) is used in this example.For the temperature field, the inlet temperature is prescribed (T in = 293.15K), and the other boundary conditions are adiabatic (− → n • q = 0).For the velocity and pressure field, the velocity magnitude, turbulent kinetic energy, and turbulent energy dissipation at the inlet are given (U in = 0.5 m/s, k 0 = 0.005 m 2 /s 2 , ε 0 = 0.005 m 2 /s 3 ), and the pressure at the outlet is set as p out = 0.The wall condition is non-slip, and the turbulent flow in the near-wall region is simulated by the wall function method.The design domain is discretized into 13,427 four-node quadratic elements and the initial value of the design variable γ is set to 0.5.
The hydraulic and heat transfer performance of topology optimization channels is compared with the traditional straight channel in Section 3.2.The traditional straight channel structure is shown in Figure 3, where the length and width of the fins are 80 mm and 7.8125 mm, respectively.The volume fraction of the fins is the same as the solid region of topology optimization models.The hydraulic and heat transfer performance of topology optimization channels is compared with the traditional straight channel in Section 3.2.The traditional straight channel structure is shown in Figure 3, where the length and width of the fins are 80 mm and 7.8125 mm, respectively.The volume fraction of the fins is the same as the solid region of topology optimization models.

Comparison of Several Conjugate Heat Transfer Topology Optimization Designs
We carry out topology optimization using three different methods, namely the Navier-Stokes flow model, the Darcy flow model, and our k-ε-based model, under the same problem settings, and compare their outcomes.In the Navier-Stokes flow model, the penalty term is only added to the momentum equation but not to the turbulent kinetic energy equation and turbulent energy dissipation equation.Convergence histories of the objective function for all three methods are shown in Figure 4, and the total number of iterations and the computation time are summarized in Table 2.As can be seen in these results, the Navier-Stokes flow model requires the fewest iterative steps, but each step is time-consuming, and the total computing time is the longest.Using the Darcy flow model largely  The hydraulic and heat transfer performance of topology optimization channels is compared with the traditional straight channel in Section 3.2.The traditional straight channel structure is shown in Figure 3, where the length and width of the fins are 80 mm and 7.8125 mm, respectively.The volume fraction of the fins is the same as the solid region of topology optimization models.

Comparison of Several Conjugate Heat Transfer Topology Optimization Designs
We carry out topology optimization using three different methods, namely the Navier-Stokes flow model, the Darcy flow model, and our k-ε-based model, under the same problem settings, and compare their outcomes.In the Navier-Stokes flow model, the penalty term is only added to the momentum equation but not to the turbulent kinetic energy equation and turbulent energy dissipation equation.Convergence histories of the objective function for all three methods are shown in Figure 4, and the total number of iterations and the computation time are summarized in Table 2.As can be seen in these results, the Navier-Stokes flow model requires the fewest iterative steps, but each step is time-consuming, and the total computing time is the longest.Using the Darcy flow model largely

Comparison of Several Conjugate Heat Transfer Topology Optimization Designs
We carry out topology optimization using three different methods, namely the Navier-Stokes flow model, the Darcy flow model, and our k-ε-based model, under the same problem settings, and compare their outcomes.In the Navier-Stokes flow model, the penalty term is only added to the momentum equation but not to the turbulent kinetic energy equation and turbulent energy dissipation equation.Convergence histories of the objective function for all three methods are shown in Figure 4, and the total number of iterations and the computation time are summarized in Table 2.As can be seen in these results, the Navier-Stokes flow model requires the fewest iterative steps, but each step is time-consuming, and the total computing time is the longest.Using the Darcy flow model largely reduces the computation time.However, its final objective function value is larger than that of other flow models.When the k-ε turbulent flow model is used, the computation time is longer than using the Darcy model but still within an acceptable range.The penalty terms imposed on the turbulent kinetic energy equation and turbulent energy dissipation equation accelerate the transition of the design variables to the solid or liquid phase in the topology optimization process.Moreover, the k-ε model yields the lowest objective function value thanks to the additional penalty terms.Considering both the      The evolution of design variables and the velocity distribution during the optimization process are shown in Figure 5.There are only two flow channels in the design domain in the Navier-Stokes flow model.Such a simple flow channel configuration cannot effectively improve the heat transfer performance of the conjugated heat transfer structure.When the Darcy flow model is used, the shape of the solid region looks like multiple fins.The fluid is evenly distributed by the fins at the inlet and flows into different flow channels which are nearly parallel.In the k-ε turbulent flow model, the solids in the design domain are mostly distributed in blocks, so the flow channel has more branches to enhance the heat transfer effect between the fluid and the solid.

Validation of Topology Optimization Results
The k-ε-optimized flow channel and Darcy-optimized flow channel are compared with the traditional straight flow channel.A uniform heat source (Q = 10 7 W/m 2 ) is applied to the design domain, and the boundary conditions are the same across different models.The k-ε turbulent flow physical model and conjugate heat transfer physical model are used to evaluate the hydraulic and heat transfer performance of the resulting designs.

Grid Independence Test
Firstly, a grid independence test is conducted, and the k-ε-optimized channel is used for the test.The model is considered grid independent when the deviation of the numerical results is within 2%.The test results are shown in Table 3.It can be concluded that when the number of the grid is about 300,000, the balance between calculation accuracy and computing time can be satisfied.Therefore, Grid 2 is used for the k-ε-optimized flow channel generated.The number of grids in the Darcy-optimized flow channel is 382,437, and that in the traditional straight channel is 330,378.

Flow Characteristics Analysis
Figure 6 shows the pressure drop of the different types of flow channels with different inlet velocities.It is clear that, regardless of the design of the flow channel, the pressure drop increases with increasing inlet velocity.With the same inlet velocity, the pressure drop of the topology-optimized flow channel is obviously lower than that of the straight flow channel.Moreover, the k-ε-optimized flow channel has a lower pressure drop than the Darcy-optimized flow channel when the inlet velocity is the same.The higher the inlet velocity, the larger the difference in pressure drop between the three types of channels.With the highest inlet velocity of 1.5 m/s, the pressure drop of the straight flow channel is 2760 Pa, while the pressure drop is only 1110 Pa for the Darcy design and 837 Pa for the k-ε design.Therefore, the flow channel designed by topology optimization can significantly decrease the pressure drop and the energy dissipation during fluid flow.
The relative difference in pressure drop δ is used to demonstrate the hydraulic performance improved by different types of topology-optimized flow channels compared with the conventional channel, as is shown in Figure 6.δ is defined as: where the ∆P Str and ∆P To are the pressure drop of the conventional straight flow channel and topology-optimized flow channel, respectively.It can be seen that using the Darcyoptimized flow channel reduced pressure drop by 31.79-59.79%,depending on the inlet velocity, compared with the straight flow channel.Using the k-ε-optimized flow channel reduced the pressure drop more effectively by 56.13-69.68%.Notably, although the inlet velocity selected for the topology optimization design is 0.5 m/s, the k-ε-optimized flow channel has a better hydraulic performance at all inlet velocities, and especially at higher inlet velocities.Therefore, it can be considered that topology-optimized designs using the k-ε turbulent model have applicability beyond their design condition.Figure 7 shows the streamline profile at a low inlet fluid velocity (0.5 m/s) and a high inlet velocity (1.5 m/s) in different types of flow channels.For the straight flow channel, vortices with different sizes can be found in the inlet and outlet manifolds of the parallel flow channels.This flow pattern indicates significant pressure loss and energy dissipation in the fluid flow process whatever the inlet velocity.For the Darcy-optimized flow channel, in the case of low inlet velocity, it appears that the abrupt expansion between the inlet and the heat transfer region generates large vortexes and hinders the fluid flow, which indicates higher local energy loss and an increase in the overall pressure drop.When the inlet velocity is 1.5 m/s, the fluid is concentrated in the centerline of the structure.When the fluid flows near the outlet, the sudden shrinkage of the channel also causes vortices and some local pressure loss.In addition, when the high-speed fluid hits the outlet wall, it flows back toward the inlet along the uppermost channel, which obviously hinders effective heat exchange.As for the k-ε-optimized channel design, it can be seen that the bends with fillets are formed at the bifurcation and confluence of channels.The flow is slightly disturbed, and the velocity distribution is uniform and stable.Only a few small vortices are generated when the velocity increases.These streamline results very well explain the pressure drop differences in Figure 6.The k-ε-optimized flow channel needs a lower pressure drop under the same inlet velocity largely due to less energy dissipation in vortices.Figure 7 shows the streamline profile at a low inlet fluid velocity (0.5 m/s) and a high inlet velocity (1.5 m/s) in different types of flow channels.For the straight flow channel, vortices with different sizes can be found in the inlet and outlet manifolds of the parallel flow channels.This flow pattern indicates significant pressure loss and energy dissipation in the fluid flow process whatever the inlet velocity.For the Darcy-optimized flow channel, in the case of low inlet velocity, it appears that the abrupt expansion between the inlet and the heat transfer region generates large vortexes and hinders the fluid flow, which indicates higher local energy loss and an increase in the overall pressure drop.When the inlet velocity is 1.5 m/s, the fluid is concentrated in the centerline of the structure.When the fluid flows near the outlet, the sudden shrinkage of the channel also causes vortices and some local pressure loss.In addition, when the high-speed fluid hits the outlet wall, it flows back toward the inlet along the uppermost channel, which obviously hinders effective heat exchange.As for the k-ε-optimized channel design, it can be seen that the bends with fillets are formed at the bifurcation and confluence of channels.The flow is slightly disturbed, and the velocity distribution is uniform and stable.Only a few small vortices are generated when the velocity increases.These streamline results very well explain the pressure drop differences in Figure 6.The k-ε-optimized flow channel needs a lower pressure drop under the same inlet velocity largely due to less energy dissipation in vortices.

Heat Transfer Characteristics Analysis
The heat transfer performance of different flow channels is evaluated by comparing the temperatures of the designed structure.Figure 8 shows the average temperature and maximum temperatures with different inlet velocities for the three cooling structures.For the k-ε-optimized flow channel and the straight flow channel, the average temperature and maximum temperature decrease as the inlet velocity increases.In contrast, the temperature in the Darcy-optimized channel initially rises and then decreases with an increasing inlet velocity, possibly due to the occurrence of the reversed flow in the uppermost channel as mentioned above.When the inlet velocity is lower than 0.5 m/s, the coolant mainly flows through the channels that are close to the inlet, so the heat transfer of the flow channels farther away from the inlet is poor, and local hot spots are generated.This disadvantage is exacerbated at the velocity of 0.5 m/s.Generically, when the inlet flow velocity increases, the convection between the fluid and the solid is enhanced, so the average temperature and the maximum temperature both decrease, although at the cost of a higher pressure drop.The average temperature and maximum temperature of the conjugate heat transfer structure with the k-ε-optimized flow channel are lower than those of the Darcy-optimized flow channel and straight flow channel under the same inlet flow rate.The maximum temperature of the Darcy flow channel is as high as 485 K and reduces to 336 K at the highest flow rate.In comparison, the maximum temperature stays below 328 K for the design by the k-ε model, indicating a much better cooling performance in terms of eliminating hot spots.

Heat Transfer Characteristics Analysis
The heat transfer performance of different flow channels is evaluated by comparing the temperatures of the designed structure.Figure 8 shows the average temperature and maximum temperatures with different inlet velocities for the three cooling structures.For the k-ε-optimized flow channel and the straight flow channel, the average temperature and maximum temperature decrease as the inlet velocity increases.In contrast, the temperature in the Darcy-optimized channel initially rises and then decreases with an increasing inlet velocity, possibly due to the occurrence of the reversed flow in the uppermost channel as mentioned above.When the inlet velocity is lower than 0.5 m/s, the coolant mainly flows through the channels that are close to the inlet, so the heat transfer of the flow channels farther away from the inlet is poor, and local hot spots are generated.This disadvantage is exacerbated at the velocity of 0.5 m/s.Generically, when the inlet flow velocity increases, the convection between the fluid and the solid is enhanced, so the average temperature and the maximum temperature both decrease, although at the cost of a higher pressure drop.The average temperature and maximum temperature of the conjugate heat transfer structure with the k-ε-optimized flow channel are lower than those of the Darcy-optimized flow channel and straight flow channel under the same inlet flow rate.The maximum temperature of the Darcy flow channel is as high as 485 K and reduces to 336 K at the highest flow rate.In comparison, the maximum temperature stays below 328 K for the design by the k-ε model, indicating a much better cooling performance in terms of eliminating hot spots.Figure 9 shows the spatial temperature distribution of the structure when the inlet velocity is 0.5 m/s and 1.5 m/s, respectively.Compared with the other two channel configurations, the temperature distribution of the Darcy-derived design is highly non-uniform showing high local temperature gradient.This can be explained by the flow patterns with large vortices and reversed flow as shown in Figure 7.For the k-ε-optimized channel design, the fluid flow has better attachment with the fin in the heat transfer region, and there are no apparent hot spots.

Uin=0.5m/s
Uin=1.5m/s Figure 9 shows the spatial temperature distribution of the structure when the inlet velocity is 0.5 m/s and 1.5 m/s, respectively.Compared with the other two channel configurations, the temperature distribution of the Darcy-derived design is highly nonuniform showing high local temperature gradient.This can be explained by the flow patterns with large vortices and reversed flow as shown in Figure 7.For the k-ε-optimized channel design, the fluid flow has better attachment with the fin in the heat transfer region, and there are no apparent hot spots.velocity is 0.5 m/s and 1.5 m/s, respectively.Compared with the other two channel configurations, the temperature distribution of the Darcy-derived design is highly non-uniform showing high local temperature gradient.This can be explained by the flow patterns with large vortices and reversed flow as shown in Figure 7.For the k-ε-optimized channel design, the fluid flow has better attachment with the fin in the heat transfer region, and there are no apparent hot spots.

Weighting Factor Ratio wi:wj Analysis
Decreasing the pressure drop requires the fluid to reduce friction and collision with the solid region as much as possible during the flow while improving the heat transfer performance means increasing the contact surface between the fin and the fluid, which will lead to the increase in the flow resistance and degradation of hydraulic performance.Hence, it is difficult to improve heat transfer performance and hydraulic performance

Weighting Factor Ratio w i :w j Analysis
Decreasing the pressure drop requires the fluid to reduce friction and collision with the solid region as much as possible during the flow while improving the heat transfer performance means increasing the contact surface between the fin and the fluid, which will lead to the increase in the flow resistance and degradation of hydraulic performance.Hence, it is difficult to improve heat transfer performance and hydraulic performance simultaneously.In this section, the influence of the weighting factor ratio of the heat transfer and flow objective functions on the topology optimization design is discussed.
Figure 10 shows the topology optimization results, the streamline contour, and temperature distribution using different weight factor ratios.It can be seen that with the increase in the weight factor of the heat transfer objective function w i , the number of the branches and fins of the design increases, and the flow channels become narrower.With the same solid volume, the total length of the convective heat transfer boundary is 665.08 mm, 726.79 mm, and 904.72 mm in the cases of 1:4, 1:2, 1:1, respectively.Therefore, increasing the heat transfer weight factor extends the heat transfer boundary between the fluid and the fin, thereby enhancing the heat transfer performance.Figure 11 presents the pressure drop and the average temperature when different weighting factors are used.With the increase in the weight factor ratio, the average temperature of the heat transfer region gradually decreases.However, this better heat transfer performance is obtained at the expense of a higher pressure drop and pumping power consumption.

Three-Dimensional Flow Channel Model Analysis
Three-dimensional models have more degrees of freedom, so the topology optimization for three-dimensional structures needs drastically more computation time.At the same time, the complex 3D topology optimization designs will significantly increase the difficulty of manufacturing.In this section, the 2D topology optimization design is extruded into a 3D model, and the hydraulic and heat transfer performance of the 3D extruded structure in the turbulent regime will be discussed.mm, 726.79 mm, and 904.72 mm in the cases of 1:4, 1:2, 1:1, respectively.Therefore, in-creasing the heat transfer weight factor extends the heat transfer boundary between the fluid and the fin, thereby enhancing the heat transfer performance.Figure 11 presents the pressure drop and the average temperature when different weighting factors are used.With the increase in the weight factor ratio, the average temperature of the heat transfer region gradually decreases.However, this better heat transfer performance is obtained at the expense of a higher pressure drop and pumping power consumption.fluid and the fin, thereby enhancing the heat transfer performance.Figure 11 presents the pressure drop and the average temperature when different weighting factors are used.
With the increase in the weight factor ratio, the average temperature of the heat transfer region gradually decreases.However, this better heat transfer performance is obtained at the expense of a higher pressure drop and pumping power consumption.Figure 12 shows the 3D conjugate heat transfer model of the k-ε-optimized flow channel and straight flow channels.The geometry in the XY plane is the same as the design in 3.2 (Figure 2).The height of the flow channel in the Z direction is 10 mm, and the solid base region is 2 mm thick.The bottom of the heat sink is in contact with a heat source, and the boundary condition is uniform heat flux.The inlet and outlet conditions of the fluid are the same as in Section 3.2.The interface between the fluid and solid domain, and the upper and lower walls of the flow channel are given non-slip boundary conditions.Symmetrical boundary conditions are used to simplify the computation process.
The pumping power represents the external energy required to drive the fluid through the cooling channels.Therefore, it is used to analyze the hydraulic characteristics of the heat exchanger.Pumping power is defined as: where the .
V is the volume flow of the fluid.The pumping power represents the external energy required to drive the fluid through the cooling channels.Therefore, it is used to analyze the hydraulic characteristics of the heat exchanger.Pumping power is defined as: where the  ̇ is the volume flow of the fluid.
Figure 13 shows the pumping power variation for the two channel designs at different inlet velocities.For both the k-ε-optimized and the straight flow channels, as the flow rate increases, the pumping power consumed by the fluid flow increases rapidly.Compared to the straight flow channel, the k-ε-optimized flow channel consumes less pumping power at the same inlet velocity.This advantage of the k-ε-optimized flow channel is more pronounced when the inlet velocity is higher.The total thermal resistance and Nusselt number are used to compare the heat transfer performance of the 3D model [34].The total thermal resistance can be defined as follows: R th =  surf, max −   (22) where the  surf, max is the maximum temperature of the surface of the conjugate heat transfer model.
The Nusselt number is defined as: where the characteristic length Dh depends on the cross-sectional area of cooling channel  ch and the wetted perimeter   and the relation of them can be described as: h avg is the average convective heat transfer coefficient, according to Newton's law of cooling, it is defined as: The total thermal resistance and Nusselt number are used to compare the heat transfer performance of the 3D model [34].The total thermal resistance can be defined as follows: where the T surf, max is the maximum temperature of the surface of the conjugate heat transfer model.The Nusselt number is defined as: where the characteristic length D h depends on the cross-sectional area of cooling channel A ch and the wetted perimeter C w and the relation of them can be described as: Entropy h avg is the average convective heat transfer coefficient, according to Newton's law of cooling, it is defined as: where T w,avg refers to the average wall temperature.T f ,avg is the average fluid temperature and can be calculated as follows: Figure 14 shows the variation of the total thermal resistance and Nusselt number of the heat transfer structure under different flow velocities, respectively.The total thermal resistance decreases and the Nusselt number increases with the increase in inlet velocity.When the fluid velocity in a flow channel becomes higher, the thickness of the thermal boundary layer is reduced and the heat transfer between the fluid and the fin is enhanced.For all the inlet velocities used here, the total thermal resistance of the heat sink with the k-ε-optimized flow channel is lower and the Nusselt number is higher than that of the straight flow channel.Therefore, it can be concluded that the 3D k-ε-optimized flow channel improves the convective heat transfer performance of the heat sink compared with the conventional straight channel.Figure 15 shows the 3D temperature distribution and streamlined distrib ferent heat transfer structures with an inlet velocity of 1 m/s.Similar to the re 2D models (Figure 7), for the heat sink of the straight flow channel, the abrup or contraction of the flow channel at the inlet and outlet of the flow channe vortices and leads to local kinetic energy loss.On the other hand, for the hea the k-ε-optimized flow channel, the fins evenly distribute the fluid in the flo Moreover, the temperature distribution of the entire heat sink is relatively un pared with the straight channels.
Figure 16 shows the temperature distribution at the bottom of the heat sin with the hot surface.For the heat sink with the straight flow channel, the aver ature is 315.49K and the maximum temperature is 328.83K.For the heat sink ε-optimized flow channel, the average temperature is 314.12 and the maxim ature is 317.73K, indicating a cooler and more uniform temperature profile.Figure 15 shows the 3D temperature distribution and streamlined distribution of different heat transfer structures with an inlet velocity of 1 m/s.Similar to the results of the 2D models (Figure 7), for the heat sink of the straight flow channel, the abrupt expansion or contraction of the flow channel at the inlet and outlet of the flow channel generates vortices and leads to local kinetic energy loss.On the other hand, for the heat sink with the k-εoptimized flow channel, the fins evenly distribute the fluid in the flow channels.Moreover, the temperature distribution of the entire heat sink is relatively uniform compared with the straight channels.
Figure 16 shows the temperature distribution at the bottom of the heat sink in contact with the hot surface.For the heat sink with the straight flow channel, the average temperature is 315.49K and the maximum temperature is 328.83K.For the heat sink with the k-ε-optimized flow channel, the average temperature is 314.12K and the maximum temperature is 317.73K, indicating a cooler and more uniform temperature profile.channels.Moreover, the temperature distribution of the entire heat sink is relatively uniform compared with the straight channels.
Figure 16 shows the temperature distribution at the bottom of the heat sink in contact with the hot surface.For the heat sink with the straight flow channel, the average temperature is 315.49K and the maximum temperature is 328.83K.For the heat sink with the k-ε-optimized flow channel, the average temperature is 314.12K and the maximum temperature is 317.73K, indicating a cooler and more uniform temperature profile.

Conclusions
In summary, we have developed a stable topology optimization scheme for conjugate heat transfer involving single-phase flow in the turbulent regime, and its effectiveness is verified by numerical simulation.From the above discussion, we can draw the following conclusions: Although the article proposes a topology optimization method for turbulent conjugate heat transfer, the method still has some limitations.For example, the variable density method is unable to capture the changes in the solid-liquid interface during the optimi-

Conclusions
In summary, we have developed a stable topology optimization scheme for conjugate heat transfer involving single-phase flow in the turbulent regime, and its effectiveness is verified by numerical simulation.From the above discussion, we can draw the following conclusions: The heat transfer structure with the k-ε-optimized flow channel has superior hydraulic and heat transfer performance compared to conventional straight channels and Darcybased TO designs.Elimination of undesirable flow vortices is believed to be an important factor for this performance.

3.
A multi-objective function considering the flow resistance as well as the thermal performance is used, and different weight factors are tested.As the objective function puts more weight on heat transfer, the optimized flow channels become narrower and more branched.4.
The k-ε-optimized flow channel has better hydraulic and heat transfer performance in the practical application of the 3D heat sink, which reduces the pumping power and thermal resistance during the heat transfer process and effectively increases the Nusselt number.

20 Figure 2 .
Figure 2. Schematic of the design model and computational settings.

Figure 2 .
Figure 2. Schematic of the design model and computational settings.

Figure 2 .
Figure 2. Schematic of the design model and computational settings.

Entropy 2023 ,
25, 1299 8 of 19 optimization efficiency and the final optimization object function, the k-ε model yields the best optimization results among the 3 models used.tion process are shown in Figure 5.There are only two flow channels in the design domain in the Navier-Stokes flow model.Such a simple flow channel configuration cannot effectively improve the heat transfer performance of the conjugated heat transfer structure.When the Darcy flow model is used, the shape of the solid region looks like multiple fins.The fluid is evenly distributed by the fins at the inlet and flows into different flow channels which are nearly parallel.In the k-ε turbulent flow model, the solids in the design domain are mostly distributed in blocks, so the flow channel has more branches to enhance the heat transfer effect between the fluid and the solid.

Figure 6 .
Figure 6.Pressure drop and the relative difference in pressure drop between different flow channel designs.

Figure 6 .
Figure 6.Pressure drop and the relative difference in pressure drop between different flow channel designs.

Figure 7 .
Figure 7.The distribution of streamline in 0.5 m/s and 1.5 m/s.

Figure 8 .
Figure 8.The temperature variation of different flow channels under different inlet velocities.(a) Average temperature (b) maximum temperature.

Figure 9 .
Figure 9. Temperature distribution of different flow channels in 0.5 m/s and 1.5 m/s.

Figure 9 .
Figure 9. Temperature distribution of different flow channels in 0.5 m/s and 1.5 m/s.

Figure 10 .
Figure 10.The k-ε-optimized flow channel, the distribution of streamline, and the temperature with different weighting factor ratios wi (heat transfer):wj (pressure drop).

Figure 11 .Figure 10 .
Figure 11.The pressure drop and the average temperature with different weighting factor ratios.

Figure 10 .Figure 11 . 1 Figure 11 .
Figure 10.The k-ε-optimized flow channel, the distribution of streamline, and the temperature with different weighting factor ratios wi (heat transfer):wj (pressure drop).

Entropy 2023, 25 , 1299 15 of 19 3. 2 (
Figure 2).The height of the flow channel in the Z direction is 10 mm, and the solid base region is 2 mm thick.The bottom of the heat sink is in contact with a heat source, and the boundary condition is uniform heat flux.The inlet and outlet conditions of the fluid are the same as in Section 3.2.The interface between the fluid and solid domain, and the upper and lower walls of the flow channel are given non-slip boundary conditions.Symmetrical boundary conditions are used to simplify the computation process.

Figure 13 20 Figure 13 .
Figure13shows the pumping power variation for the two channel designs at different inlet velocities.For both the k-ε-optimized and the straight flow channels, as the flow rate increases, the pumping power consumed by the fluid flow increases rapidly.Compared to the straight flow channel, the k-ε-optimized flow channel consumes less pumping power at the same inlet velocity.This advantage of the k-ε-optimized flow channel is more pronounced when the inlet velocity is higher.Entropy 2023, 25, x FOR PEER REVIEW 16 of 20

Figure 13 .
Figure 13.The pumping power consumed by different channel designs.

Entropy 2023 ,
25,  x FOR PEER REVIEW channel improves the convective heat transfer performance of the heat sink with the conventional straight channel.

Figure 14 .
Figure 14.The variation of the total thermal resistance and Nusselt number: (a) Tota sistance (b) Nusselt number.

Figure 14 .
Figure 14.The variation of the total thermal resistance and Nusselt number: (a) Total thermal resistance (b) Nusselt number.

Figure 15 .
Figure 15.The temperature and velocity distribution of different conjugate heat transfer models.Figure 15.The temperature and velocity distribution of different conjugate heat transfer models.

Figure 15 . 20 Figure 16 .
Figure 15.The temperature and velocity distribution of different conjugate heat transfer models.Figure 15.The temperature and velocity distribution of different conjugate heat transfer models.Entropy 2023, 25, x FOR PEER REVIEW 18 of 20

1 .
It is crucial to add appropriate penalty terms to the turbulent kinetic energy equation and turbulent energy dissipation equation of the k-ε model, respectively.Compared with the Navier-Stokes flow channel and Darcy-optimized flow channel, the k-ε-optimized channel can better balance the relationship between calculation time and design requirements.2. The heat transfer structure with the k-ε-optimized flow channel has superior hydraulic and heat transfer performance compared to conventional straight channels and Darcy-based TO designs.Elimination of undesirable flow vortices is believed to be an important factor for this performance.3. A multi-objective function considering the flow resistance as well as the thermal performance is used, and different weight factors are tested.As the objective function puts more weight on heat transfer, the optimized flow channels become narrower and more branched.4. The k-ε-optimized flow channel has better hydraulic and heat transfer performance in the practical application of the 3D heat sink, which reduces the pumping power and thermal resistance during the heat transfer process and effectively increases the Nusselt number.

Figure 16 .
Figure 16.The temperature distribution of the heating surface.

1 .
It is crucial to add appropriate penalty terms to the turbulent kinetic energy equation and turbulent energy dissipation equation of the k-ε model, respectively.Compared with the Navier-Stokes flow channel and Darcy-optimized flow channel, the k-ε-optimized channel can better balance the relationship between calculation time and design requirements.2.

Table 2 .
The values of computation results.

Table 2 .
The values of computation results.

Table 3 .
Grid independence test results.