Improved Thermal Performance of a Serpentine Cooling Channel by Topology Optimization Infilled with Triply Periodic Minimal Surfaces

: The present study utilizes a density-based topology optimization method to design a serpentine channel under turbulent ﬂow, solving a high pressure loss issue and enhancing heat transfer capability. In the topology optimization, the k – ε turbulence model is modiﬁed by adding penalization terms to reveal turbulence effects. Heat transfer modeling is included by setting the modiﬁed energy equation with additional terms related to topology optimization. The main objective is to minimize pressure loss while restricting heat transfer. The 2D simpliﬁed model is topologically optimized. Then, the optimal solution with intermediate results is extruded in the 3D system and interpreted with triply periodic minimal surfaces (TPMS) to further enhance heat transfer performance. Compared to the baseline serpentine channel, the optimized model inﬁlled with the diamond-TPMS structure lowers pressure loss by 30.8% and signiﬁcantly enhances total heat transfer by up to 45.8%, yielding thermal performance of 64.8% superior to the baseline. The temperature uniformity is also improved. The simulation results show that the curvatures in the optimized model with diamond-TPMS structure eliminate the large recirculation ﬂow and low heat transfer regions. This model diminishes the effect of Dean’s vortices but promotes high turbulent kinetic energy, leading to better uniform ﬂow and heat transfer distributions.


Introduction
Applications of a multipass cooling channel under turbulent flow and heat transfer enhancement themes can be found in an internal cooling gas turbine blade, microelectronic cooling, etc.Generally, the serpentine channel is a standard internal cooling scheme in practical turbine blades, effectively utilizing the coolant in the midchord region [1,2].However, complicated flow fields induced by sharp turns after bend regions cause a nonuniform flow distribution.Rao et al. [3] have reported that the turning regions cause more than 90% of the total pressure drop of the multipass channel.These phenomena in the serpentine channels result in highly nonuniform-flow heat transfer and high pressure loss, decreasing overall thermal performance.Therefore, it is desirable to design multipass channels that generate a lower-pressure loss while maintaining the total heat transfer and achieving better uniformity.
In order to control the flow field and heat transfer distribution in serpentine channels, placing the turning guide vane at the bend region is a common method that has been examined in several previous studies.Schüler et al. [4] and Guo et al. [5] investigated experimentally and numerically different turning guide vane designs in the multipass channel.They found that the pressure loss was reduced by about 13.4%-48.4%when using the turning guide vanes, but the heat transfer deteriorated.Similarly, Lei et al. [6] observed that the pressure loss was largely reduced when installing the turning vane, yet heat transfer decreased because the flow separation was suppressed.This design process was based on a trial-and-error method.To achieve low pressure loss while maintaining total heat transfer in a serpentine cooling channel, various geometrical parameters were varied, causing a cost-and time-consuming design process.
Alternatively, optimization methods can speed up and ease the designs.Layout optimization is usually divided into size, shape and topology optimization (TO).Size and shape optimizations have been used to minimize pressure loss and enhance thermal performance in heat exchangers.The results reveal that the optimal design significantly reduces the pressure drop and achieves higher thermal efficiency than the original model [7].However, it is observed that the design process only adjusts the positions of points or boundaries of the channel to obtain the objective functions without changing the topology of the structure.Moreover, since additive manufacturing is rapidly developed, which can compete for limitations to produce relatively complex geometries, the size and shape optimizations cannot completely exploit this manufacturing scheme's potential.Meanwhile, finding the optimal shapes and optimizing the topology of channels can be obtained simultaneously through topology optimization.Additionally, free-form and efficient geometries can be generated, satisfying solutions for additive manufacturing.Therefore, topology optimization is a more general design method than size and shape optimizations [8].
Topology optimization has been employed in different design systems covered by various disciplines.It has been well developed and affirmed widely due to competitive designs in cost and functionality under various constraints [8][9][10][11][12][13][14].Several methods, such as the homogenization approach, the evolutionary technique, the level-set and the density approach, have been developed for implementation.
Among topology optimization schemes, the density-based method has gained more popularity [15].In the comparisons of Manuel et al. [16], the density-based approach yields better outcomes for most cases in terms of iteration numbers.Additionally, the density method merits rapid convergence and weak dependence on the initial guess [8,16].Furthermore, it could be highlighted that no prior research has considered turbulent flow systems except the density-based method.The key idea of the density-based method is the penalized interpolation function, which continuously forces fictitious porous cells, including fluid and solid, to be fluid or solid subdomains.A comprehensive review and its general formulation of topology optimization can be found in previous work [8].Here, a brief overview of the density-based topology optimization applied in recent research on thermofluid systems is shown in Table 1.
Table 1.An overview of density-based topology optimization in recent studies on thermofluid systems.
As highlighted in Table 1, most previous studies are 2D models with a steady-state laminar flow regime.The FEM is the most implemented discretization method since the topology optimization initiates from solid mechanics.Moreover, it is easier for the discrete adjoint scheme to manage the FEM [15].
Nevertheless, only limited works on turbulent flow topology optimization have been studied, which is primarily needed in most industrial flow applications, e.g., a turboshaft engine's turbine blade serpentine cooling channels.Several works employed a frozen turbulence assumption to simplify sensitivity analysis for the topology optimization of Energies 2022, 15, 8924 3 of 23 turbulent flows [23,24].However, as the frozen turbulence neglects the effect of porosity in the turbulent viscosity calculation, the optimized model negatively affects the accuracy.The shortcomings of this assumption have been presented by Dilgen et al. [25].
Recently, study of the turbulence effects in the topology optimization for turbulent flow has appeared.Dilgen et al. [25] optimized many 2D and 3D complex flow channels using Spalart-Allmaras and k-ω turbulence models, where k denotes turbulent kinetic energy and ω is the specific dissipation rate.The results showed that the generated structure in a U-bend channel significantly reduced the pressure drop across the channel.Later, Yoon [26] proposed the topology optimization with the FEM using the k-ε turbulence model, where ε is the turbulent energy dissipation.However, it is noted that the energy equation for topology optimization has not yet been found, and convective heat transfer optimization has not yet been explored in those studies.
Even though topology optimization is an efficient tool for generating complex structures, it still causes some problems.The results of density-based topology optimization may convert a number of intermediate densities.This type of conversion causes pressure maldistribution in the cooling channel between the design and its realization by additive manufacturing, which could not explore the full performance of the optimized structure.To overcome the issues, the solid domain of the optimized model, including the intermediate results, can be fulfilled with a graded lattice structure.This method has been observed to help manage the intermediate densities in the conventional topology optimization, reduce the structure's weight, and preserve the original design structure [27]; however, this approach is rarely applied to cooling channel applications.
Among complex lattice structures, the repeating network, called triply periodic minimal surface (TPMS), has attracted much research interest due to its compact and significant thermal properties.Additionally, the TPMS is a self-supporting structure, which has been observed to present better mechanical performance than other strut-based lattice structures [28].The TPMS networks have also shown better heat transfer performance in cooling channels than the lattice-frame material due to a smoother surface and larger wetted area [29].The TPMS, with its complex and interconnected structures, provided consistent fluid streamlines, promoting permeation [30].Moreover, many previous works have found that the TPMS structures considerably increase heat transfer with reasonably increased pressure loss, improving the thermal performance of the cooling channel [31][32][33].
This work proposes an optimized design to minimize pressure and increase heat transfer, improving the thermal performance of a serpentine cooling channel.The density-based topology optimization for the thermofluid problem under turbulent flow is firstly utilized to minimize the pressure loss while maintaining heat transfer in the serpentine cooling channel of a turbine blade.The turbulence effects are revealed by adding penalization terms in the k-ε turbulence model programmed by COMSOL.The functions to activate/deactivate the convective and conduction terms in the modified energy equation are also included to account for the heat transfer modeling in the topology optimization.Three-dimensional (3D) serpentine models with the TPMS structures are obtained by infilling the TPMS structures in the solid domain and the intermediate density resulted from the conventional topology optimization.The final models are resimulated in ANSYS FLUENT to show elaborated flow structure, pressure loss and heat transfer advantages over the baseline and the channel generated from the conventional topology optimization.By infilling the TPMS structures in the conventional topology optimization, the intermediate result is eliminated, the original channel structure is preserved, and the heat transfer is further increased due to the complex structure and large wetted area of the TPMS.

Model Description
The design model is inspired by an internal cooling channel of a turboshaft heat exchanger.The channel dimension of the numerical model in this work is consistent with the experimental serpentine cooling channel for a turbine blade with short passes, studied by Guo et al. [5]. Figure 1 illustrates the three-short-pass serpentine channel for (a) a full-scale 3D model and (b) a simplified 2D design.Each channel pass has a cross-section of 120 mm × 60 mm, yielding a hydraulic diameter (D h ) of 80 mm and an aspect ratio (AR) of 2.0.The length of all the straight passages is 245 mm, resulting in a length-to-hydraulic diameter of 3.06.Before entering the heating section, entrance and exit sections extended with a length of 300 mm are created at the channel inlet and outlet, respectively.The first turning region is designed as a rectangular shape with an inner radius of 20 mm since this region near the tip of the turbine blade suffers from a high heat load, which needs a large heat transfer area to cool.Meanwhile, the second bend near the blade platform receives a lower heat load; hence it is designed as an arc shape with outer and inner diameters of 140 and 20 mm, to reduce the total pressure loss in the channel [5].Typically, the topology optimization requires hundreds of iterations to reach the convergence of the final optimal design.Hence, to save the cost of obtaining the optimized design, the 2D design simplified from the original 3D model is firstly used in this work [19][20][21], as demonstrated in Figure 1b.An inlet boundary condition (Γ in ) is set as uniform velocity (U in ) with a prescribed temperature (T in ).The inlet velocity is defined based on the experimental Reynolds number (Re), as follows: where ρ is fluid density and µ is dynamic viscosity.An outlet boundary condition (Γ out ) is set as a zero static pressure (p 0 ).The boundary condition of the heating surfaces on the top and bottom is imposed with uniform heat flux (q 0 ).The remaining walls are thus set as no-slip conditions (Γ wall ) to surround the serpentine channel.Figure 1b also shows the design domains, in which all the heating regions inside design domains (Ω id ) are evolved.Meanwhile, the extended entrance and exit of the channel set as outside design domains (Ω od ) are not considered for optimization.

Thermofluid Modeling
Previous works have considered an incompressible flow for thermofluid topology optimization to design turbine blade cooling channels [24].Moreover, it has been observed that when the Mach number is less than 0.3, the deviation of the compressible and incompressible flow can be neglected.Since the maximum bulk velocity of air in the channel is lower than 15 m/s, the Mach number is less than 0.1, then incompressible and steady-state flow are assumed throughout this work.

Fluid-Flow Modeling
The continuity equation and Navier-Stokes equation assuming the incompressible steady fluid are expressed as follows: where u indicates the fluid velocity vector; ∇ signifies the gradient operator; p denotes the pressure field; I is the identical tensor; µ T is the turbulent dynamic viscosity; * is the transpose; and F represents the Brinkman friction term introduced by Borrvall and Petersson [34], calculated from Equation (4).Here, the body force is not exerted on the fluid outside the design domain (F = 0 in Ω od ).
where α u,max is the maximum inverse permeability of the porous medium, while α u,min is the minimum one.I u (γ) denotes the inverse permeability interpolated function stated in Section 3.1.
From Equation ( 5), α u,max can be defined depending on Darcy number (Da):

Turbulence Modeling
As the topology optimization progressively evolves the design domain to allow material distribution, modifying the turbulence model is also necessary.In the present work, the authors modify the standard k-ε model in COMSOL by adding penalty terms to reveal the effects of turbulence on the serpentine cooling channel for topology optimization, which is modified as follows: The closure coefficients and all parameters in Equations ( 6) and ( 7) are described as follows: where P k signifies the turbulent kinetic energy source term.The last terms in Equations ( 6) and ( 7) are augmented to enforce k and ε to k 0 and ε 0 for the solid regions in the design domain.They are defined as follows [26]: where α k,max and α ε,max are the maximum limits of the penalization for k and ε, respectively, while α k,min and α ε,min are the minimum ones.I k (γ) and I ε (γ) are material interpolated functions in both Equations ( 8) and ( 9).Here, the values of α k,min , α ε,min , k 0 and ε 0 are assigned to zero, while the value of α k,max and α ε,max are set equal to the value of α u,max .The material interpolations are detailed in Section 3.1.

Heat Transfer Modeling
The thermal convection-diffusion equation outside the design domain is given as follows: where c f denotes the fluid heat capacity, k f indicates the fluid thermal conductivity and T is the temperature field in the serpentine channel.
In the design domain, as the distribution of the fluid and solid domains is determined after the optimization, the authors build up the energy equation by adding the functions to control the convective and conduction terms as follows: where I c (γ) is the function used to activate/deactivate the convective term for fluid (I c (γ) = 1) and solid (I c (γ) = 0).I K (γ) is used to interpolate between the thermal conductivity of the fluid and the thermal conductivity of the solid (k s ).Lastly, the term q 0 indicates the constant bottom wall heat flux.
It should be noted that q 0 is the uniform heat flux applied to the heating boundaries in the 2D model.Additionally, the energy equation is weakly coupled to the fluid flow problem to obtain u; then, it substitutes the energy equation to solve T. Furthermore, the Reynolds number and heat flux applied in this work are fixed at Re = 10,000 and q 0 = 1500 W/m 2 , respectively, because varying the velocity or the heat input in the topology optimization could change the topology of the structure in the final result [11,14].

Optimization Procedures
The procedure for obtaining the optimized model with lattice structures in this study is depicted in Figure 2. Firstly, the 2D simplified model, shown in Figure 1b, is topologically optimized using the boundary conditions and governing equations described in Section 2. The optimal solid domain, including the intermediate densities, is then extruded into the 3D model and interpreted with the TPMS structures.Finally, the final designs are recalculated for comparison.
Overview of the procedures in this study.

Topology Optimization
In this section, the details of the topology optimization process to obtain the 2D optimized structure are presented.

Material Distribution
The material distribution in the density-based method described by the design variable (γ) takes a value between 0 and 1 at any point in the design domain.Here, γ = 1 denotes free flow in the channel, while γ = 0 is solid regions.The discrete 0/1 problem is relaxed to use gradient-based methods.As the material distribution takes intermediate values (0 ≤ γ ≤ 1), this relaxation requires the interpolation of transition areas.The interpolations of the fluid flow equation in the design domain are expressed in Equations ( 12)-( 14): Energies 2022, 15, 8924 7 of 23 where q u , q k and q ε are tuning parameters to control the function curvature.
For the thermal convection and conductivity, the interpolation terms I c (γ) and I K (γ) are calculated using solid isotropic material penalization (SIMP) [8] in Equations ( 15) and ( 16), respectively: where n c and n K are penalization power coefficients; δ is fixed at δ = 10 −15 to avoid numerical problems when solving Equation ( 11) [22].

Problem Formulation
In this study, minimizing the pressure loss is set as the objective function (J), while the heat transfer across the serpentine channel is constrained.In minimizing the pressure loss for the cooling channel, the power dissipation within the volume is generally selected to proceed in the optimization process since it shows the flow fluctuation in the considered domain.It provides robust convergence, and better flow uniformity can be expected when minimizing this function [11,12].The power dissipation is expressed as follows: In prior works, the temperature through the boundary has been defined as another objective function to maximize heat transfer [24].This function represents the temperature gain of fluid across the channel.Since the majority of the previous work has observed that the use of the turning vanes deteriorated the heat transfer [4][5][6], the temperature constraint (C) is thus employed to restrict the heat transfer in the present work.It is written as follows: Here, Γ out is the outlet boundary, and n is the normal vector.With this expression, the heat transfer in the channel can be constrained (g) as follows: where C ref is the reference value calculated from the conventional serpentine channel.Hence, the mathematical optimization model is written in Equation ( 20): where S is the vectors with discrete state fields, R indicates the residual vector function evaluated from the discretization of the governing equations, g represents inequality constraints and x denotes the design variable vector.

Density Filter and Projection
During the iterative process, the design variable changes independently in each cell, causing a large difference value of the design variables in the adjacent cells.In the present optimization, a partial differential equation filter [35] is adopted to mitigate such phenomena, which is defined as follows: where γ f signifies a smoothed design variable.R min is a filter radius calculated by = 1.5 × maximum mesh size in the considered domain.
Additionally, a projection technique, which projects the design variable into either 0 or 1, is adopted to mitigate the intermediate results.Here, a smoothed Heaviside projection is applied as follows: where γ p is the projected design variable, β is a parameter used to control the steepness of the projection, and η denotes the projection threshold, fixed at η = 0.5 in this study.

Summary of Topology Optimization Process
The present study uses COMSOL Multiphysics to run the topology optimization process.The flow chart of the topology optimization process is demonstrated in Figure 3.After the design variables and the parameters are initialized, the governing equations are solved to attain physical variable fields.The objective and the constraint functions are then computed.If the convergence criterion is not achieved, the sensitivity is computed by the adjoint method to obtain the gradient information.The design variable field is then updated using GCMMA and regularized.When the constraint is achieved, and the relative change in the objective is less than 10 −4 , or after every 100 iterations, the optimization proceeds to the next continuation steps.The continuation approach is used to avoid the low-quality local minima result and ease the convergence in this optimization process [18].Eventually, the optimization is considered to be completed when all the continuation steps are computed. ---

Infilling Lattice Structures
After obtaining the 2D optimized model, the optimal solid domain within the channel, including intermediate results generated by the conventional topology optimization, is extruded into the 3D model and replaced by TPMS structures.The TPMS networks are used to infill the extruded optimized model with a range of thickness according to the

Results and Discussion
In this section, two main parts are presented.For the topology optimization, the problem setup and 2D simplified model validation are detailed in the first part.The optimal density distribution in the 2D optimized model is extruded into the 3D domain and fulfilled with the TPMS structures.The final designs are then recalculated and validated in the second part.The pressure drop, heat transfer and temperature uniformity are compared with the baseline serpentine channels.

Validation of the 2D Model
Since the 3D channel is simplified into the 2D model to reduce the computational cost for topology optimization, the accuracy of the stated variable field obtained from 2D calculations should be satisfied [19].Air and copper are selected in this work for fluid and solid properties for the simulation.The boundary conditions and governing equations explained in Section 2 are used to compute the flow, turbulence effect and heat transfer.It should be marked that in COMSOL, the out-of-plane heat flux (q 0 ) in the 2D model is equivalent to the uniform heat flux boundary in the 3D system.
The mesh system of the 2D and 3D models is demonstrated in Figure 5.It is suggested that the boundary layer mesh evolves during the optimization for the accurate turbulence calculation, or the mesh should be fine enough to obtain accurate and reliable results.The extra-fine mesh setting at 42,468 elements, containing mainly triangle element type, is thus generated for all 2D models, as shown in Figure 5a.Meanwhile, for the 3D model, the tetrahedral mesh is generated, and the boundary layer mesh with a prism shape is employed in the wall regions.The fine mesh system is selected at 2.38 million elements, as shown in Figure 5b.
To verify the discrepancy between the 2D and 3D models, averaged values of the pressure, velocity and temperature at 12 different positions are plotted in Figure 6.In all 2D cases, the results are averaged from 12-cut lines, whereas in all 3D cases, the outcomes are extracted from the planes, then averaged.Overall, the 2D results have similar trends to the 3D ones for all variables.The pressure values in Figure 6a and velocity magnitudes in Figure 6b for the 2D model are lower than in the 3D model, while the temperature values in Figure 6c are opposite.The 2D model overestimates the temperature by about 0.4%-2.5%,and the average error value is 1.6%.The temperature discrepancy between the 2D and 3D models can attribute to the limitation of a pseudo-3D model, which has also been observed in several investigations [13,19].However, the temperature distributions that slightly decreased along the channel path are well estimated, providing similar trends among them.The variations can be acceptable for the optimization process.Therefore, comparing results with the 3D numerical model, the 2D model can be accurately employed to perform topology optimization in this work.

Optimized Structure
The initial design is a uniform density field of γ 0 = 0.6.The model is evaluated through the continuation parameters, as given in Table 2, to find the optimal structure.After the simulation, the design converges to the optimized structure, as presented in Figure 7.It can be noticed that the model leaves the intermediate fields to satisfy the temperature constraint by distributing the material in the channel.The intermediate results are due to the high conductivity difference between the copper and air used in this work [19].For the 2D optimized design, the total pressure loss between inlet and outlet boundaries decreases by about 33.9% compared to the original 2D model, while the constraint temperature in the channel is maintained.The optimal density distribution within the 2D optimized channel, including intermediate results generated by conventional topology optimization, is extracted with ranges of density filtering from 0.1 to 0.9 and extruded into the 3D model, as shown in Figure 8a.The TPMS structures are then used to infill with a range of thickness between 2-4 mm.The lowest thickness of 2 mm is specified at the density filtering of 0.9, and the thickness is increased accordingly, reaching the highest thickness of 4 mm at the density filtering of 0.1, as demonstrated in Figure 8b.The thickness range is selected based on the printable capability in an actual scale of the cooling channel for a gas turbine blade.Here, the single unit cell size of the diamond-and gyroid-sheet structures is 120 × 120 × 120 mm 3 .

3D Model Comparisons
In this section, the optimized model extracted with a density filtering of 0.5 (TO model) is also included to compare with the baseline channel.All models and the boundary conditions are shown in Figure 9.These models are applicable at the Reynolds number Re = 10,000 and heat flux q 0 = 1500 W/m 2 according to the design settings in the conventional topology optimization, as mentioned in Section 2. For all cases, the channel dimensions are consistent with the baseline serpentine channel, as described in Figure 1a.The inlet and outlet boundaries are the same as the baseline channel.Here, the red zones are treated by the uniform heat flux, including the top, bottom and extruded areas from the optimization and TPMS structures.The remaining surfaces are set as no-slip walls.The flow and heat transfer for all models in this section are simulated using the commercial CFD solver ANSYS FLUENT 18.0, which provides lower computational costs and consumes lesser computer memory for 3D complex cooling channels than COMSOL.The overview of mesh systems for the optimized model with the diamond-sheet structure is demonstrated as an example in Figure 10.The hybrid mesh involving tetrahedron and prism elements is created to discretize the entire domain.The tetrahedron meshes are finely dense in the small-and high-curvature regions, while the boundary layer meshes with a prism shape are used near the no-slip walls.The meshes near the wall are refined to ensure that y+ values are less than one.

3D Model Validation
For the validation, the conventional three-short-pass smooth channel investigated by Guo et al. [5] is used as the referee design.The SIMPLEC algorithm is utilized for the pressure-velocity coupling, and the second-order upwind method is adopted.The convergence criterion is less than 10 −4 for the continuity, velocity and turbulence quantities and less than 10 −6 for the energy equation [5].
Three sets of the meshes at 4.95, 7.31 and 11.6 million elements are simulated to check mesh independence for the baseline channel.The numerical results are steady and converge to a certain value as the mesh elements increase.The grid convergence index (GCI) [36] is evaluated using the spanwise-averaged surface temperature to estimate the numerical accuracy of the mesh at 7.31 million, as demonstrated in Figure 11.The line denotes the highest mesh solution, while the red bars represent the error.It is found that the maximum error of the samples is less than 0.5%, and the overall GCI value on the spanwise-averaged temperature is about 0.1%, which indicates that the numerical results are independent of the meshes.Therefore, the mesh element is selected at 7.31 million for the baseline serpentine channel, while this mesh setting is also used for the other models.In order to select the turbulence model with the highest accuracy, several Reynoldsaveraged Navier-Stokes equations (RANS), i.e., SST k-ω, realizable k-ε and RNG k-ε, are executed and compared with the baseline serpentine channel.SST means the shear stress transport turbulence model, and RNG denotes renormalization group methods.Here, the globally averaged Nusselt number of each pass (Nu) is used to compare the heat transfer with the experimental data from Guo et al. [5].It can be evaluated as follows: where Q tot denotes the total heat transfer, A w is the wetted area, and ∆T signifies the log mean temperature difference of the heating wall, calculated as follows: where T w is the average wall temperature, while T in and T out are the inlet and outlet fluid temperatures.
The globally averaged Nusselt number simulated from different turbulence models for the baseline serpentine channel at the Reynolds number of 10,000 is shown in Figure 12.It should be noted that Guo et al. [5] have found that the SST k-ω turbulence model could accurately simulate the heat transfer of the three-pass channel with ribs.This is because the SST k-ω turbulence model captures well the reattachment flow from the ribbed walls.However, in this work, SST k-ω and RNG k-ε considerably overestimates the results in the second pass by 25.9% and 21.4%, respectively.Meanwhile, the average Nusselt number simulated using realizable k-ε agrees well with the experimental results for all smooth passages due to correctly predicting the large recirculation flow.The mean error deviates by about 10.7% compared to the experimental data.Hence, the realizable k-ε turbulence model with enhanced wall treatment is selected to simulate all models in this study.

Pressure Loss
The relative pressure coefficient (C p ) is employed to characterize the pressure loss from different models.The relative pressure coefficient along the channel flow path is calculated as follows: where p loc is the local static pressure at measuring points and p ref indicates the pressure value at the channel inlet.
The relative pressure coefficient at three different locations at the Reynolds number of 10,000 is demonstrated in Figure 13.In the first pass, as the structure for all models is the same without changing from the baseline channel, the values are comparable because the coolant flows smoothly with minor flow resistance.The losses are higher when the coolant enters the second and third passes.Here, the baseline serpentine channel causes the highest losses, followed by the TO model with the gyroid-TPMS structure.In contrast, the TO model provides the lowest values, which are lower than the baseline model by 14.5% and 44.3% for the second and third passages.The friction factor calculated to compare the total pressure loss in the channel can be expressed as follows: where ∆p is calculated from the pressure at the inlet boundary minus the pressure at the outlet boundary.L is the total length of the serpentine cooling channel.
The friction factor for all models at the Reynolds number of 10,000 is compared in Figure 14.All cases cause lower friction loss than the baseline channel, and the TO model provides the lowest value.It can be observed that the TO model with diamond-and TO model with gyroid-TPMS structures cause slightly higher friction loss than the TO model attributed to their complex curvatures.Compared to the baseline channel, the TO model,

Heat Transfer
The globally average Nusselt number, including two turning regions, at the Reynolds number of 10,000 is presented in Figure 15.It is evident that all topology optimization models provide higher heat transfer than the baseline channel, particularly in the TO model with the diamond-TPMS structure.This result is because the heat transfer is constrained during the optimization, and the heat transfer is further increased when interpreting the intermediate results with the TPMS structures.The globally averaged Nusselt number in the TO model, TO model with diamond-and TO model with gyroid-TPMS structures are higher than in the baseline channel by about 9.6%, 21.1% and 17.6%, respectively.
- The total Nusselt number (Nu tot ) presents the total heat transfer performance based on the flat base heating area.The total Nusselt number includes the heat transfer contribution from all wetted surfaces, as the total heat transfer is experienced by the turbine blade's suction and pressure walls, and it is calculated as follows: where A h is the flat heated area.It is noted that the total wetted area of the TO model, TO model with diamond-and TO model with gyroid-TPMS structures are increased by 8.5%, 15.6% and 13.3% compared to the baseline channel.The total Nusselt number at the Reynolds number of 10,000 is shown in Figure 16.The total Nusselt number is higher than the globally averaged Nusselt number, particularly in the TO model with the diamond-TPMS structure.This substantial increase is because the diamond-sheet structure has a larger contact area than the TO model curvature and the gyroid-sheet structure.The total Nusselt number for the TO model, TO model with diamond-and TO model with gyroid-TPMS structures increases by 25.3%, 45.8% and 36.6%,respectively, compared to the baseline channel.To show the heat transfer enhancements in the serpentine channel, the Dittus-Boelter equation [37] is employed to normalize the globally averaged and total Nusselt numbers, defined as follows: Nu 0 = 0.023Re 0.8 Pr 0.4 (28 where Pr is the Prandtl number, and the value of Pr for the air in this study is 0.7.
The globally averaged and total heat transfer enhancements at the Reynolds number of 10,000 are presented in Figure 17a,b, respectively.The enhanced heat transfer is observed for all models since the serpentine channels generate high heat transfer, particularly from the first turning region.Among the comparative models, the heat transfer enhancement is more evident for the TO model with diamond-TPMS structure since this model considerably increases total heat transfer in the channel.

Thermal Performance
To compare the thermal performance of the serpentine cooling channel, the heat transfer performance requires an assessment after involving friction loss.In this study, the thermal performance of the serpentine channels can be calculated as follows: where f 0 is the Blasius correlation, defined in Equation (30): The thermal performance for all cases at the Reynolds number 10,000 is presented in Figure 18.It is evident that the baseline serpentine channel shows the lowest thermal performance due to causing the highest pressure loss and lowest heat transfer in the channel.Meanwhile, the TO model with the diamond-TPMS structure exhibits the highest value because the optimized curvatures help minimize the pressure loss in the serpentine channel, while the diamond-sheet structure enhances high heat transfer.The thermal performance for the TO model, TO model with diamond-and TO model with gyroid-TPMS structures improves by 47.0%, 64.8% and 46.5%, respectively, compared to the baseline channel. - Thermal performances for all serpentine models at Re = 10,000.

Temperature Uniformity
The high temperature difference between leading and trailing surfaces generated in the gas turbine channel is observed when the channel rotates at high speed [1,2].It is essential to be concerned about the temperature uniformity that affects the temperature difference on the surface.The uniformity index of temperature (θ) is introduced in this section to determine surface temperature uniformity.It can be calculated as follows: where T i denotes the local temperature, and A i is the local area.T avg is the average temperature of the considered region, calculated as follows: The uniformity index of surface temperature for all serpentine channels at the Reynolds number of 10,000 is presented in Figure 19.In the first pass, all cases exhibit comparable values, while the uniformities drastically decrease in the second pass for the baseline model.These negative results are due to the strong turning effect at the first rectangular bend, causing nonuniform flow distribution throughout the second pass [5].However, the uniformities are improved in the TO model, TO model with diamond-and TO model with gyroid-TPMS structures owing to the curvatures that manipulate smoother flow distributions.For all optimized models, the temperature uniformity on the surface improves by about 4.5%-5.6% in the second passage and 0.4%-2.0% in the third passage, superior to the baseline channel.
- In short, the optimized model infilled with the diamond-sheet structure shows the best thermal performance because it provides low pressure loss and enhances the highest heat transfer.This structure also improves the temperature uniformity throughout the channel compared to the baseline model.Therefore, this study selects the TO model with diamond-TPMS structure as the optimal serpentine channel.

Detailed Flow and Heat Transfer Characteristics of the Optimal Serpentine Channel
According to the analysis of the pressure loss, heat transfer and temperature uniformity in Section 4.2, the optimized serpentine channel with diamond-sheet structure is selected to be the optimal model, and the detailed flow and heat transfer characteristics are demonstrated and compared to the baseline serpentine channel.The pressure distribution on the middle section of the channel at the Reynolds number of 10,000 is demonstrated in Figure 21 for (a) the baseline model and (b) the TO model with diamond-TPMS structure.For the baseline channel, the high-pressure loss emerges in the two bend regions, which could be attributed to the generation of the Dean's vortices [5].However, after inserting the diamond-sheet structure into the optimized serpentine channel, the pressure loss is significantly reduced.It can be observed that the pressure loss is reduced, particularly in the second passage inside the diamond-sheet structure, and the area of low pressure continues to the third passage.The TO model with the diamond-TPMS structure manipulates the fluid well, reducing the recirculation flow, which lowers the total pressure loss in the serpentine channel.
The flow mechanisms and turning effects can be discussed with the dimensionless turbulent kinetic energy (TKE/U in 2 ) and streamlines on different cross-sections of the serpentine channel, as presented in Figure 22.For both models, the turbulent kinetic energy and streamlines exhibit the same characteristics in the first passage.The higher turbulent kinetic energy values can be observed in the second and third passages.
For the baseline model, as shown in Figure 22a, the Dean's vortices can be observed on the outer wall when the fluid enters the second passage on plane B3.These vortices are prominent in plane B2, while the fluid mixes at the end of the second passage, resulting in high turbulent kinetic energy in the middle region.The coolant turning into the third passage further increases the turbulent kinetic energy and continues until it leaves the channel, as seen on planes C1-C3.Even though the strong mixing flow could enhance high heat transfer, this effect causes substantial differences in the flow distribution, eventually leading to non-uniform heat transfer in the channel.For the TO model with the diamond-TPMS structure, as shown in Figure 22b, the Dean's vortices at the inner wall, as found in the baseline model, are diminished because the large solid volume of the diamond-sheet structure distributes the flow in the second passage, as demonstrated on planes B3-B1.However, the high turbulent kinetic energy areas can be observed.The main coolant is forced to flow along the curvature of the diamond-sheet structure, intensifying high turbulent kinetic energy.Moreover, the complex flow generated inside the diamond-sheet structure expands large regions of high turbulence.The turbulent kinetic energy in the third passage is also more uniform than in the baseline model.These phenomena in the TO model with the diamond-TPMS structure lead to higher heat transfer and better uniformity than in the baseline channel.

Conclusions
In this work, the concept of density-based topology optimization is utilized to design the serpentine cooling channel under turbulent flow.In the topology optimization, the penalization terms are added to the k-ε turbulence model to reveal the turbulence effects.The functions to activate/deactivate the convective and conduction terms in the energy equation are also carried out.The main objective is to minimize the pressure drop while maintaining heat transfer.The 2D simplified models are topologically optimized.The optimized solution, including the intermediate results, is then extruded into the 3D system and infilled with TPMS structures.The final models are resimulated and validated to compare with the baseline serpentine channel.The main conclusions are drawn as follows: 1.
The conventional topology optimization model achieves a lower pressure loss and provides higher heat transfer than the baseline channel.Infilling the diamond-and gyroid-sheet structures in the optimal solution, including the intermediate results of the conventional topology optimization model, further enhance high heat transfer and maintains lower pressure loss than the baseline channel.

2.
The 3D optimized models with the TPMS structures provide lower pressure loss by about 19.0%-30.8%and higher total heat transfer by 36.6%-45.8%,compared to the 3D baseline channel.The optimized mode infilled with the diamond-TPMS structure is selected as the optimal serpentine channel in this study since it provides the best thermal performance, up to 64.8%, superior to the baseline channel.The temperature uniformity on the surface is also improved, particularly in the second passage.

3.
The 3D optimized model with the diamond-TPMS structure significantly eliminates the recirculation flow and the low heat transfer regions at the second and third passages.This model also reduces the influence of Dean's vortices but maintains high turbulence kinetic energy, leading to better uniform flow and heat transfer distributions.
Overall, the method in this work provides an advanced tool to manage the intermediate results in the conventional topology optimization of turbulent flow and heat transfer systems.Complex structures generated by this approach can be fabricated with additive manufacturing.Even though the present optimization models could be employed at a specific Reynolds number and heat input, this method can be applied to design cooling channels at higher Reynolds numbers and different cooling channel shapes, achieving higher thermal performance for gas turbine blades, heat sinks for high-power electronics, etc.  Uniform heat flux (W/m 2 ) q u , q k , q ε Tuning parameters in material interpolated functions

Figure 1 .
Figure 1.Numerical model of a three-short-pass serpentine channel: (a) full 3D model; (b) 2D simplified design from the 3D model.

Figure 3 .
Figure 3. Flow chart of the topology optimization process.

Figure 4 .
Figure 4.A single unit cell of (a) diamond-sheet structure; (b) gyroid-sheet structure for fulfilling the optimal results in this study.

Figure 5 .Figure 6 .
Figure 5. (a) The mesh system for 2D model validation; (b) mesh system in the 3D model.

Figure 9 .
Figure 9. Models for comparison: (a) optimized model extracted with density filtering of 0.5 (TO model); (b) optimized model infilled with diamond-sheet structure (TO model with diamond-TPMS); (c) optimized model infilled with diamond-sheet structure (TO model with gyroid-TPMS).

Figure 10 .
Figure 10.An example of mesh systems for the 3D models.

Figure 11 .
Figure 11.The spanwise-averaged temperature and grid convergence index on the surface of the conventional serpentine channel at Re = 10,000.

Figure 12 .
Figure 12.Comparisons of globally averaged Nusselt number in each pass for different turbulence models at Re = 10,000.

Figure 13 .
Figure 13.Comparisons of relative pressure coefficient for all serpentine channels at Re = 10,000.

Figure 14 .
Figure 14.Friction factors for all serpentine channels at Re = 10000.

Figure 15 .
Figure 15.Globally averaged Nusselt numbers for all serpentine channels at Re = 10,000.

Figure 16 .
Figure 16.Total Nusselt numbers for all serpentine channels at Re = 10,000.

Figure 17 .
Figure 17.Heat transfer enhancements for all serpentine models at Re = 10,000: (a) globally averaged Nusselt number ratio; (b) total Nusselt number ratio.

Figure 19 .
Figure 19.Temperature uniformity on the surface for all serpentine models at Re = 10,000.

4. 3
.1.Flow CharacteristicsThe velocity distributions and streamlines on the middle section (Z = 30 mm) of the serpentine channel at the Reynolds number of 10,000 are shown in Figure20for (a) the baseline channel and (b) the TO model with diamond-TPMS structure.It is clearly seen that for the baseline serpentine channel, the large recirculation flow is generated at the second bend region and third passage.The high velocity can be observed at the outer wall of the third passage.These phenomena cause nonuniform distribution of the flow.Meanwhile, the recirculation flow in the TO model with the diamond-TPMS structure is eliminated, as seen in Figure20b.The curvature in this model induces the coolant to flow near the inner wall of the second passage.The flow in the diamond-sheet structure is also distributed evenly, leading to a more uniform velocity distribution throughout the channel.

Figure 20 .
Figure 20.Velocity distributions and streamlines in the middle section of the serpentine channel: (a) baseline model; (b) TO model with diamond-TPMS structure.

Figure 21 .Figure 22 .
Figure 21.Pressure distributions in the middle section of the serpentine channel: (a) baseline model; (b) TO model with diamond-TPMS structure.
Fluid heat capacity (J/(kg K)) C p Relative pressure coefficient C ref Reference value in constraint function Da Darcy number D h Hydraulic diameter of the channel (m) F Fictitious body-force (N/m 3 ) Inequality constraints I u , I k , I ε , I c , I K Material interpolated functions J Objective function k, k 0 Turbulent kinetic energy of the fluid and solid domains (m 2 /s 2 ) k f , k s Thermal conductivity of the fluid and solid (W/(m K)) L Total length of the serpentine channel (m) n Normal vector n c , n K Penalization power coefficients in material interpolated functions Nu, Nu tot Globally and total Nusselt numbers Nu 0 Dittus-Boelter equation P k Turbulent kinetic energy source term (W/m 3 ) p 0 Zero static pressure (Pa) p loc Local static pressure at measuring points (Pa) p ref Pressure value at the channel inlet (Pa) Pr Prandtl number Q tot Total heat transfer (W) q 0 Nondimensional distance from the endwall to the first element node Greek Letters α u,max , α u,min Maximum and minimum inverse permeabilities of the porous medium (Pa s/m) β Projection slope γ Design variable γ f , γ p Smoothed and projected design variables ε, ε 0Turbulent energy dissipation for fluid and solid regions (m 2 /s 3 )