Numerical thermal–hydraulic analysis and multiobjective design optimization of a printed circuit heat exchanger with airfoil overlap fin channels

In this study, we numerically investigated the performance of a printed circuit heat exchanger (PCHE) featuring National Advisory Committee for Aeronautics (NACA) 0020 airfoil overlap fins integrated into its straight flow channels. Numerical simulations were performed by solving Navier–Stokes and energy equations to analyze liquid water at temperatures ranging from 30 to 90°C, which is possible for use with a geothermal heating system, for example. The optimization of the PCHE airfoil fin geometry is performed using nondominated sorting genetic algorithm II (NSGA‐II) with three design features: front, back, and overlap lengths. The effectiveness and pressure drop of the PCHE are evaluated as two objective functions. The present study shows that, when a Reynolds number falls between 1088 and 5000, the overlapping airfoils exhibit lower pressure drops and higher convective heat flux relative to the separate airfoils, and an increase in overlap length always decreases the pressure drop positively. It is also found that using overlapping airfoils of front length of 0.1 mm, back length of 0.119 mm, and overlap length of 0.203 mm leads to almost the same effectiveness and a 22% decrease in pressure drop (which is generally believed to preferably maintain the high heat transfer performance and minimize pressure drop, as a favorable design), compared to a base design of the airfoils with front and back lengths of 0.2 mm and an overlap length of 0.1 mm.

this lost heat as possible. 1,2Heat exchangers effectively transfer heat between the geothermal fluid and a working fluid, and they are used for this purpose.Among all heat exchangers, printed circuit heat exchangers (PCHEs) [3][4][5] have received more attention than conventional plate heat exchangers (PHEs) due to their high compactness and effectiveness. 6,7raditional heat exchangers are often made by layering metal, leading to fluid leaks due to perforations caused by corrosion (Fan et al. 8 ).This issue can be resolved by diffusion bonding, which is a manufacturing process where plates are joined under high-mechanical-pressure and high-temperature conditions to create seamless metal layers. 4,9This process allows PCHEs to operate at high temperatures (reaching 1173 K) and pressures (reaching 60 MPa). 10 Diffusion bonding eliminates the need to heat the metal to a liquid phase or to use welding materials, reducing plastic deformation and increasing mechanical strength. 7Through diffusion bonding, PCHEs can be manufactured in a practical and efficient manner.
Integrating airfoil profiles into PHEs or PCHEs can improve thermal performance by increasing the heat transfer area and maintaining low pressure loss due to the streamlined airfoil design.Many studies have been conducted on the use of airfoils within PHEs.First, by introducing an airfoil fin, Kim et al. 11 optimized the thermal and hydraulic performance of a PCHE and showed that the pressure loss per unit heat transfer rate of the airfoil fins is approximately one-twentieth of that of zigzag fins.Tang, Cui, and Sundén 12 investigated the effects of layouts with four different airfoil configurations-National Advisory Committee for Aeronautics (NACA) 0015, NACA 0018, NACA 0021, and NACA 0024-on the thermal and hydraulic performance levels of PCHEs, showing that NACA 0024 can be exploited to achieve the highest overall heat transfer performance.By considering a microgas turbine recuperator for extended-range electric vehicles, Wang et al. 13 proposed a multiobjective optimal design of NACA airfoil with three targets-heat transfer rate, pressure drop, and compactness-based on nondominated sorting genetic algorithm II (NSGA-II).Pidaparti et al. 14 experimentally explored the thermal-hydraulic performance of a set of chemically etched discontinuous offset rectangular and airfoil fin surface patterns for supercritical CO 2 power cycles, proposing empirical correlations for the friction factor and the Nusselt number on a systematic basis of the experimental data.Zhang et al. 15 presented an airfoil fin PCHE, indicating that the airfoil fin geometry has a high hydraulic performance; the pressure drop is approximately 1/6 of that of a zigzag PCHE with a comparable heat transfer rate.Nguyen et al. 16 performed a numerical study on an airfoil corrugated PHE with commercial and geometrical parameter analysis over Reynolds number ranges of 4600-8400, showing that NACA 0030 airfoil profile can be used to achieve the highest thermal performance, with an overall heat transfer coefficient that is 12% higher than those of the NACA 0020 and NACA 0025 profiles; the researchers reported that a pressure drop increases by approximately 33.7% for NACA 0030.In an experimental study, Chang et al. 17 investigated the heat transfer performance of asymmetric airfoil-fin PCHEs based on the average thermal-resistance ratio, reporting that the reverse flow barely improves the heat transfer and significantly deteriorates the pressure drop.Since four various microchannels with airfoil fins, circular fins, trapezoidal fins, and triangular fins are employed for a counter flow heat exchanger, Cao et al. 18 found it favorable to incorporate beryllium oxide (BeO) ceramic as the microchannel heat exchanger material, showing that the optimum structure is obtained by a microchannel with trapezoidal fins; the compelling thermal features of BeO allow a more uniform temperature distribution, with an exciting 219% improvement in the effectiveness.Nevertheless, room remains for improvement regarding effectiveness and pressure drop reduction with novel airfoil fin structural modifications for PCHEs, providing additional benefits for industrial applications including waste heat recovery systems in geothermal engineering.
With the inclusion of zigzag, sinusoidal, airfoil, S-shaped, rhombic and sinusoidal fin channel structures, as reported in Tsuzuki et al., 19 Xie et al., 7 Yang et al., 20 and Liu et al., 21 heat transfer enhancement is often accompanied by sufficient mixing, resulting from changes in fluid flow direction.However, flow diversion around bending points can result in a pressure drop.The careful design of channel structures while considering geometric factors, such as bending angle, entrance length, and fillet radius, and/or thermoentropic factors is crucial for optimization, and it can be facilitated with robust optimization methods.Due to NSGA-II, which is an effective optimization method that was first proposed by Deb et al. 22 to demonstrate the strong ability to gauge the Pareto frontier with precise and relatively fast convergence, 23 special attention has recently been paid to maximizing the efficiencies of heat exchangers.For example, Hadibafekr et al. 23 employed NSGA-II to perform multiobjective optimization of a wavy lobed heat exchanger tube, showing the thermoentropic analysis and thermal performance of a heat exchanger built with a wavy twisted tri-lobed tube.By exploring the NSGA-II optimization of geometrical parameters in PCHEs, Jin et al. 24 indicated that the appropriate core length of the zigzag channel changes with the changing design of the outlet temperature.Soleimani and Eckels 25 developed a multiobjective optimization of 3D microfins using NSGA-II, presenting a unique intersection of the drag reduction approach and heat transfer enhancement; the researchers found unique operating geometries and a particular optimal design geometry at Re = 49,000, where drag decreases by 8% and heat transfer increases by 21% relative to the smooth surface.Sai and Rao 26 employed a hybrid optimization method, incorporating NSGA-II, to optimize the design of a shell-and-tube heat exchanger, reporting that the hybrid method has a 4.85% and a 1.51% lower total cost in their two cases.By utilizing the Slipcevie method and experimental verification with one hybrid method of NSGA-II and multiobjective particle swarm optimization, Xu et al. 27 built a performance prediction model.Therefore, NSGA-II is still an efficient alternative for practical heat exchanger applications.
The aims of this study are to investigate the effects of overlapping airfoil fin channels on PCHEs and optimize their design through NSGA-II with three parameters: front, back, and overlap lengths.The motivation for this study is the lack of research on the effect and optimization of the overlapping airfoils (as illustrated in Figure 1) and their potential for creating a more uniform temperature distribution and possibly improving heat transfer relative to regular separate airfoils.The organization of the present paper is shown in the following sections.In Section 2, the mathematical model for the heat transfer exchanger, the governing equations, the boundary conditions, and the NSGA-II descriptions are shown.In Section 3, the results and discussions are described.In Section 4, the summary and future outlook are presented.

PHYSICAL AND MATHEMATICAL MODEL DESCRIPTIONS AND OPTIMIZATION
The entire computational domain considered in Figure 2 is separated into two main regions: the solid domain is made of stainless steel-UNS S30415 (upper panel) because of its thermal conductivity, and the fluid domain (lower panel), including two cold stream channels (shown in blue) and one hot stream channel (shown in red), which is between the two cold plates, represent the empty space that allows for the flow of liquid water for convective heat transfer.Within each of the three layers of the flow channels, there are three lines of 98 airfoil fin structures, as shown in Figure 2. The total length of our model is based on the straight channel from Xie et al., 7 while the pitch ratio of the fins is determined by referencing the article by Yang et al. 28 In the present study, the airfoils are arranged for each of the three flow plates in two manners.In the first scheme, airfoils are placed with equal spacing distributions (Figure 2B).In the second scheme, the trailing portion of every airfoil overlaps with the leading edge of its adjacent portion into alignment, as shown in Figure 2C.The former is treated as a benchmark design to show the thermal-hydraulic performance comparison with the latter, as discussed in Section 3.1.It is important to note that the Reynolds number of both schemes are set at 1088, if not otherwise specified, and the inlet temperatures of the hot and cold streams are set at 90 and 30 • C, respectively.When choosing the inlet temperatures for both hot and cold streams, we carefully considered the temperatures observed in real hot springs within the temperature range employed in a relevant study conducted by Maruoka et al.Reducing the spacing between airfoils is possible to enhance heat transfer, albeit at the expense of increased pressure drop compared to overlapping airfoils.To investigate this further, we modified the horizontal spacing between the fins from 1.84 mm (see Figure 2B) to 0.42 mm while keeping other conditions constant.The results reveal that reducing the horizontal spacing improves heat transfer; however, it also results unfavorably in a 29% increase in pressure drop compared to the overlap base design.As a result, we decided to use the spacing of the fins of 1.84 mm in the present study.
Because the size of the PCHE with fin structures is sufficiently small, the incompressibility of the working fluid is physically ensured for the fluid domain in the present study.Accordingly, the steady state mass, momentum transport, and energy conservation without an external heat source or gravitational effect are described mathematically as follows: where  is the fluid density, v is the fluid velocity, K is the viscous stress, p is the static pressure, k f is the fluid thermal conductivity, T f is the fluid temperature, and c p is the specific heat capacity of the fluid.For the solid domain, the energy equation is referred to as thermal conduction.
where k s is the solid thermal conductivity.The conjugate conditions on the solid/fluid interface are adopted to provide temperature continuity for the thermal fields, heat fluxes of a body and flow near the interface.Given a high Reynolds number, reaching approximately several thousands, we consider the fluid within the turbulence regime; the fluid requires the use of a turbulence model to calculate the fluctuations of flow motion.In the present study, the shear stress transport (SST) turbulence model with a two-equation eddy-viscosity model 30,31 is exploited to solve typical engineering problems with two parameters, such as turbulence kinetic energy k and dissipation of the eddies ω; the SST model is suitable for properly capturing the effects of the adverse pressure gradient that often occurs in the flow fields of PCHEs, as reported by Xie et al. 7 with experimental validation.In other words, the SST model is commonly acknowledged for effectively incorporating the benefits of handling free streams derived from the k-ε model, 32 along with the crucial characteristics of accurate numerical treatment for near-wall regions from the k-ω model. 33To illustrate the effects of the use of the SST turbulence model, the results obtained for two simulations under otherwise identical conditions are shown in Figure 3.In the first simulation, a relatively small flow eddy near the object is observed in laminar flow.In the second simulation, clear eddy structures appear in the SST turbulent model near the concave area or neck (represented geometrically by the trailing part of an airfoil that overlaps the head of an identical airfoil) on either side.Meanwhile, regardless of whether we turned on the turbulent modeling, we found that variations in the outlet temperatures, collected from the present simulations and used to calculate the effectiveness in Equation ( 7) the pressure drop defined in Equation ( 8), were always small.Notably, in Section 3.1, the SST turbulence model is considered in simulations referring to the thermal-hydraulic analysis, and in Section 3.2, throughout NSGA-II optimization, we consider no turbulence model to avoid the failure of convergence due to the poor mesh resolution near the airfoil overlap surface when a change in the overlap length parameter occurs.
Wang et al. 34 noted that Tsuzuki et al. 19 reported the generation of local eddy currents and a significant increase in stress due to fluid flow in zigzag flow channels, which could potentially jeopardize equipment safety.On the other hand, incorporating rounded corners in the steering parts of the zigzag flow channels leads to a decrease in both eddy currents and stress phenomena. 34In the present study, the absence of sharp convex corners near the low-lying areas of the overlapping airfoils results in the reduced stress, and low-speed eddy currents are shown in Figure 3.The two essential features are generally believed to make equipment operation less susceptible to damage, thereby ensuring safety.
The resulting turbulence governing equations are shown below.
where P k is the effective rate of production of k, σ k and β * 0 are constants, μ t is turbulent viscosity ) σ ω , β, and σ ω2 are constants, P ω is the rate of production of ω, and f v1 is the blending function.These parameters associated with the turbulent model are those same as the ones used in our previous study in Xie et al. 7

Computational methods and boundary conditions
To solve Equations ( 1)-( 6), the use of finite element analysis to calculate the velocity and temperature fields is implemented through the commercial code COMSOL Multiphysics.Both the fluid and solid domains are decomposed into small elements, including tetrahedrons, pyramids, prisms, triangles, and quads, which are used to accurately represent the curved airfoil interfaces between the fluid and solid domains; this representation allows us to ensure an adequate number of elements (more than 1,657,925 elements) for meshing the three main geometrical sections (two fluid domains and one solid domain) in a manner that always guarantees grid independence results by collecting the data on the exit temperature of the hot plate and the pressure drop.During optimization, the computational domains are re-meshed automatically within COMSOL Multiphysics whenever the NSGA-II optimization changes into one of the geometric design factors considered herein; this factor is changed through SolidWorks, which is used to build the present 3D PCHE model for numerical simulation.The fully developed flow boundary condition is ensured in the inlets of both hot and cold streams in the present study to represent the more or less actual condition, and the periodic boundary condition is prescribed at the lateral boundaries of the computational domain, due to actual heat exchangers that have more than one set of flow channels.Additionally, we use the pseudo time-stepping approach 35 to accelerate convergence rather than the Anderson's algorithm, 36 which allows us to have slightly faster convergence; the former approach makes us always ensure convergence in the simulations, which differs from the latter, even when considering the SST turbulence model.

Model validation
The present model validation is carried out with the experimental data from Liu et al., 37 where there are 100 total straight microchannels made of stainless steel-316L and arranged as a counter flow on the hot and cold sides.Due to the fluid flowing evenly into each channel, only one pair of flow channels (one for hot fluid and the other for cold fluid) is established in this validation, as shown in Figure 4.The hot working fluid is supercritical carbon dioxide and the cold one is liquid water.The properties of the supercritical carbon dioxide refer to the data from NIST fluid data.A fully developed flow is ensured at the inlets of both hot and cold streams, and the lateral sides of the computational domain are set with periodic boundary conditions, as both used in Liu et al. 37 The SST turbulence model is used in the validation with the size and stretching factor of the mesh thickness of the first boundary layer of 0.01 mm and 1.2, respectively. 37And the other Input parameters used for numerical validation are listed in Table 1.Given that Reynolds numbers ranging between 1088 and 5000, considered in this study, result in a relatively thin viscous sublayer of the boundary layer, making it challenging to use a sufficient number of grid points to resolve the viscous sublayer (Ferziger and Peri ć38 ), we selected a suitable mesh, as depicted in Figure 4.The thickness of the first boundary layer in the mesh corresponds to a y + value of approximately 1.This choice enables us to adequately represent the viscous sublayer in turbulent flow 39 and ensure numerical convergence.The results show that the simulation results are in good agreement with the experimental results, 37 with the relative errors of around 1.5% and 2.4% for the outlet temperatures of the hot and cold steams, respectively.
F I G U R E 4 Geometrical configuration of a single layer of straight channels with counterflow arrangement for validation.

Multiobjective optimization using NSGA-II
Due to a desire to solve the multiobjective optimization problem, multiobjective evolutionary algorithms are created, [40][41][42] which involve the Pareto front to obtain the best solution through mutual constraint.This elitist strategy that ensures a better solution is included throughout the iterative process.The first generation of NSGA has the following negative points 22 : (1) the objective function of each individual's solution must be compared with every other function at each dominated level, thus making the computation expensive; (2) a sharing parameter has to be specified; and (3) the rudimentary NSGA typically lacks an effective elite strategy to retain good solutions.NSGA-II was made properly in response to these three disadvantages.Since the improvement, there is no longer a need to share parameters during computation; in nondominated sorting, each solution is compared once, thus increasing the sorting efficiency.
Providing the maximum heat transfer and minimum pressure loss is central to heat exchanger design.An increase in heat transfer with the use of wavy flow channels typically results in large pressure losses (i.e., large frictional losses); small pressure loss represented by the flowing fluid through straight channels decreases heat transfer due to the insufficient contact surface area.As a result, two opposing objective functions, that is, effectiveness and pressure drop, are considered simultaneously in the present study as follows: where Q actual is the actual heat transfer rate, Q max is the maximum heat transfer rate, ṁ is the mass flow rate, c ph is the specific heat capacity of the hot fluid, c p,min is the lowest specific heat capacity, and T h i , T h o , and T c i represent the temperatures of the hot inlet, hot outlet, and cold inlet, respectively, p in is the inlet pressure, and p out is the outlet pressure.The Pareto optimal front 43 must weigh two objective functions so that the corresponding optimal PCHE geometry can be obtained.As discussed in Section 3.1, the feasible design variables that are responsible for determining the performance of PCHEs are the front, back, and overlap lengths, as illustrated in Figure 2.For the current optimization, the front length and back length range between 0.1 and 0.3 mm, respectively, while the overlap length varies between 0.05 and 0.55 mm.It is important to note that neither the front length nor the back length can exceed 0.4 mm.This limitation is primarily due to the presence of narrow passages between the overlapping airfoils that may impede the flow of working fluid, which may cause heat transfer performance deterioration.Therefore, we restrict the range of both the front length and the back length to be between 0.1 and 0.3 mm.
The NSGA-II algorithm 22 starts with the creation of the initial population based on the defined constraints.Then, based on the nondomination condition, sorting the population that is subsequently ranked with front numbers is performed.By following the obtained distances and ranks and by making use of a binary tournament selection where the selection criterion is based on a less crowding distance and low domination rank, parents are chosen from the population.The new parents undergo crossover to generate a new offspring population to mutate afterwards.Finally, the parent population is replaced with the best individuals of the population.This process is repeated until the maximum number of iterations is reached.The abovementioned overall procedure is shown in the flowchart, as shown in Figure 5, and the present algorithm is implemented through MATLAB code.
Generally, nondominated and crowding distance sorting are central to NSGA-II, 22 as shown in Figure 6.The essential points are summarized below.After mutation, the offspring population O n is formed, creating a combination of the original parent population P n and the offspring population.This combination is twice the parent population size.After nondominated sorting and crowding distance sorting are applied, the nondominated set (Rank 1) is the best solution among the population.If the size of Rank 1 is smaller than the population size, all of the members of Rank 1 are selected into the new population P n+1 ; the rest of the empty members in P n+1 come from Rank 2, Rank 3, and so forth, until the whole new population is fully accommodated. 22If P n+1 does not provide sufficient space to accommodate all members of Rank L, the better solutions within a less crowded region (according to crowding distance sorting) are selected for P n+1 . 22

RESULTS AND DISCUSSION
In the present study, the NACA 0020 airfoil is employed in the simulations.The first two digits of NACA 0020 (00) represent the symmetry of the airfoil shape, whereas the last two digits (20) represent the percentage of the chord that is composed of the airfoil thickness.And the present PCHEs with and without airfoil overlap of the chord length discussed in this section refer to the corresponding configurations shown in Figure 2B 28 The mass flow rates are determined based on the Reynolds numbers considered in Table 2.

Thermal-hydraulic analysis
The flow behavior around each of the non-overlapping and overlapping airfoils of the PCHE is a complex phenomenon that involves the interactions of various physical parameters.In particular, the velocity distribution plays a crucial role in determining the heat transfer efficiency of an airfoil used in PCHEs.Figure 7 shows the advantages of using overlapping airfoils within PCHEs by comparing the total heat flux, velocity contour, and streamline patterns.Both PCHEs equipped with separate and overlapping airfoils have two small passages between the upper wall and the first row of the airfoils and between the bottom wall and the last row of airfoils.These long narrow passages typically have high-speed areas, resulting in high heat flux and extensive heat transfer enhancement.The maximum fluid velocity is near the area of maximum thickness for each airfoil.As the thickness of the airfoil decreases, the flow velocity gradually decreases around either side of its rear part due to the reduction in airfoil cross-sectional area that the flow encounters, thus creating a gradual decrease in the velocity of the flow.High flow speed areas are observed in both the separate and overlapping airfoils, which supposedly create additional heat transfer.In the present study, the high-speed domains of the overlapping airfoils are wider than those of the separate airfoils, particularly in the spaces between the upper walls and first rows of the airfoils and between the bottom walls and last rows of the airfoils, increasing the effectiveness by 3.7% and maintaining the pressure drop at a similar level.Since the high velocity areas of the separate airfoils are shorter than those of the overlapping airfoils, adding overlapping airfoils increases heat transfer.Thus, introducing overlapping airfoils is likely to enhance heat transfer.Interestingly, the concept of total heat flux, as shown in Figure 7, can be broken down into two components: convective and conductive heat flux.A noteworthy observation can be made from Figure 8, where the downstream areas of each high-conduction airfoil are clearly larger than the rest of the areas.In general, the convective heat flux magnitude is  significantly greater than the conductive heat flux magnitude.This phenomenon is attributed to the high flow velocity and strong heat convection that the flow channels with overlapping airfoils generate.The magnitude of the convective heat flux, which is three orders of magnitude greater than that of the conductive heat flux, is indicative of its dominance over the conductive heat flux.Hence, the enhancement in heat transfer occurs primarily due to the convective heat flux rather than the conductive heat flux.The observed changes in the fluid flow patterns resulting from variations in the Reynolds number are noteworthy.As shown in Figure 9, there is a clear reduction in the low fluid velocity area when the Reynolds number increases from 1088 to 4000.This effect is most pronounced near the surfaces of the upper and bottom rows of the overlapping airfoils, where the low-speed boundary layer is barely seen.Interestingly, the low-speed areas of Re = 4000 thickens near the rear surfaces of the middle rows of the overlapping airfoils.Even so, the high speed areas dominantly bring in more heat transfer, thus resulting in higher heat flux.As a result, the overall effectiveness of the heat exchange process increases from 0.867 to 0.97, representing an improvement reaching 11.9%.This result underscores the effectiveness of introducing overlapping airfoils in enhancing heat exchange when the Reynolds number increases.Overall, these observations demonstrate the importance of the Reynolds number in determining the fluid flow patterns and heat exchange efficiencies in overlapping airfoil configurations.

F I G U R E 10
The correlation between convective heat flux and Reynolds number for both the overlapping and separate airfoils considered in the present study, using a base design of f l and b l both equal to 0.2 mm, and O l of 0.1 mm for the overlapping airfoils.
Figure 10 further shows the relationship between convective heat flux and Reynolds number for both the overlapping and separate airfoils considered in the present study.The graph visually demonstrates that as the Reynolds number rises, the convective heat flux also increases.Notably, this trend is more prominent in the case of overlapping airfoils than in the case of separate airfoils.As the Reynolds number continues to climb, the distinction between separate airfoils and overlapping airfoils becomes increasingly apparent, highlighting the clear advantage of utilizing overlapping airfoils, particularly at the high Reynolds numbers considered herein.
As shown in the previous paragraphs, the velocity distribution in the flow around the airfoil is crucial in determining the effectiveness of the present PCHE.In the present study, the overlapping airfoil shape is essentially made through the three geometric factors-the front length (f l ), the back length (b l ), and the overlap length (O l )-as depicted in Figure 2. Let us take a close look at changes in these lengths which can impact the effectiveness of the present PCHE (as given by Equation ( 7) and shown in Figure 11).
Figure 11 indicates that greater front and back lengths can enhance effectiveness, as demonstrated by the red area in the top-right corner of each panel.In other words, without being affected by the overlap length, an increase in either the front length f l or back length b l generally results in a clear increase in effectiveness.Meanwhile, as depicted in Figure 2, the length of overlap between two airfoils, referred to as O l , can also impact the effectiveness of heat transfer, as shown in Figure 11.Interestingly, the area of high effectiveness, denoted by the deep red region, is the largest when O l equals 0.2 mm.This result indicates that a well-designed O l of 0.2 mm can possibly provide superior heat transfer performance.Nonetheless, as the overlap length surpasses 0.2 mm, the area of high effectiveness gradually diminishes due to a decrease in the heat transfer area of the overlapping airfoils.When O l reaches 0.5 mm, the region of high effectiveness values exceeding 0.88 completely disappears.Since the relationship between the increase in overlap length and its impact on effectiveness is not always straightforward, as reported in Figure 11, it is essential to find a balance with the three geometric factors to achieve optimal heat transfer performance.
On the other hand, the peak pressure is typically observed at the head of the airfoil, followed by a gradual decrease in the flow pressure along the mainstream direction.This decrease generates a detrimental increase in pressure loss, 44,45 which then creates a high pressure drop.The high pressure drop generated in a system can create an unfavorable situation in which the working fluid is likely to be pressurized through a pump before entering the heat exchanger to maintain the high flow speed at the outlet, leading to an undesirable consequence: an increase in electricity used to power the pump.
The correlation between pressure drop and Reynolds number for both the overlapping and separate airfoils is demonstrated in Figure 12.As Reynolds number increases, pressure drop always increases.When Reynolds number is 1088, there is no discernible difference in pressure drop between the overlapping and separate airfoils.However, at a given Reynolds number falling between 1088 and 5000, the overlapping airfoils exhibit lower pressure drops compared to the separate airfoils.The smaller pressure drops due to the overlapping airfoil's low-lying area (see Figure 1) appear as if ones look at the similar equivalent of a dimpled golf ball that reduces more drag at a certain Reynolds number range than a smooth ball 46 or a lobed hailstone falling in air with various levels of surface roughness that shows a reduction in drag within a specific range of Reynolds number. 47This suggests that using overlapping airfoils can help to reduce pressure drop more effectively than using separate airfoils in the current study.
Changes in the front length (f l ) and the back length (b l ), and the overlap length (O l ) can also impact the pressure drop (as given by Equation 8) are shown in Figure 13.It is revealed that an increase in overlap length favorably impacts decreases in the high pressure drop area, as represented by the deep red region of each panel of Figure 13.As the overlap length increases, the area of high pressure drop gradually decreases; this phenomenon is evident by the gradual change in the blue color.Specifically, when the overlap length reaches 0.5 mm, the peak pressure drop decreases to approximately 115,000 Pa, due to small hollow places in the surface of the overlapping airfoils, indicating that increasing the overlap length can be advantageous for reducing the pressure drop, but the decrease in front and back lengths can offset the  benefits of increasing the overlap length.Therefore, it is necessary in the present study to carefully balance the tradeoffs between the overlap length and the front and back lengths to maximize heat transfer enhancement and minimize the resultant pressure drop to achieve the desired heat transfer performance.

Multiobjective optimization via NSGA-II
The present numerical simulation results, as discussed in Section 3.1, indicate that increasing the maximum thickness of the airfoil in general improves the heat transfer effectiveness; while the front length f l and back length b l of the overlapping airfoils increase, the relationship between the increase in overlap length and the increase in effectiveness is not always possible.Moreover, the improvement in effectiveness often comes at the expense of pressure loss, 48 which is also reported in the results of the separate airfoils reported in previous literature. 12,49,50Therefore, a tradeoff exists between the heat transfer effectiveness and pressure loss in the overlapping airfoils considered herein with the three geometric parameters (i.e., f l , b l , and O l ).It is crucial to strike a balance between these competing factors to optimize the performance levels of the overlapping airfoils in practical applications.The effectiveness and pressure drop of the present PCHE are employed as two objective functions.Additionally, according to Deb et al., 22 the NSGA-II's initial population of the present study is generated randomly, ensuring that it is sufficiently capable of finding optimal results within a reasonable timeframe.Figure 14 shows the Pareto front obtained through the present study.In the context of Pareto optimization, the solutions identified as Pareto-optimal are considered to be non-dominant. 51This implies that there is no single solution that can be considered superior to the Pareto-optimal solution in terms of both objective parameters (i.e., ε and Δp) at the same time.Essentially, the Pareto-optimal solution represents a trade-off between the objective parameters, where improving one objective parameter would lead to a reduction in the other objective parameter.Therefore, the Pareto-optimal solution represents the best possible outcome when considering both objective parameters simultaneously.Observing Figure 14, it is noticeable that the Pareto front curve appears to be slightly increasing until the effectiveness reaches 0.877.However, after the effectiveness surpasses 0.885, the curve becomes steeper, causing the Pareto-optimal data points to be dispersed.This confirms that heat transfer enhancement can increase with pressure drop.Table 3 lists some optimal solutions selected from the Pareto front based on both effectiveness and pressure drop criteria.
As presented in Table 3   effectiveness and an 82% increase in pressure drop relative to the use of an overlap base design ( = 0.862 and Δp=78,152 Pa) of the airfoils with f l = b l = 0.2 and O l = 0.1 mm, which represents that a system with the best effectiveness (or heat transfer) performance ( = 0.897) can be developed with a penalty of worse pressure loss (Δp=142,185 Pa), which is not preferred because relying on a pump that consumes a lot of electricity to pressurize the working fluid and maintain a high level of pressure at the inlet is uneconomical.When minimizing pressure loss, the use of overlapping airfoils of f l = 0.1, b l = 0.119 and O l = 0.203 mm leads to almost the same effectiveness and a 22% decrease in pressure drop compared to the same base design.The latter optimal design that minimizes pressure loss is deemed to be more practical, mainly because it can effectively maintain high heat transfer performance (ε = 0.851) while minimizing pressure drop (Δp = 61,000 Pa).Moreover, Table 3 shows that changes in pressure drop are more sensitive than changes in effectiveness.This occurs because modifying the overlap length value substantially impacts the pressure drop.

CONCLUSIONS AND OUTLOOK
In this study, we presented a thermal-hydraulic analysis and the multiobjective design optimization of a PCHE with airfoil overlap fin channels for improving heat transfer.In the flow channels of the PCHE, the overlapping airfoils consisted of the trailing part of one airfoil that overlapped with the head of another airfoil.To optimize the airfoil fin geometry of the PCHE, three design geometric features were considered: front, back, and overlap lengths (i.e., f l , b l , and O l ).These features were analyzed through a parametric study and then through the NSGA-II.The current optimization involves front and back lengths within the respective ranges of 0.1 to 0.3 mm, and an overlap length that varies from 0.05 to 0.55 mm.The effectiveness and pressure drop of the PCHE were measured as the two objective functions to evaluate the performance of the PCHE with overlapping airfoils.
The key points of this study were as follows: 1. Heat transfer enhancement in general occurred near either side of the maximum thickness of the front and rear airfoils, when f l and b l become greater.2. It is not always possible that an increase in the overlap length (O l ) creates a consistent increase in heat transfer in the present study, as the effectiveness gradually decreases owning to a decrease in the heat transfer area of the overlapping airfoils, when the overlap length surpasses 0.2 mm. 3.An increase in overlap length (O l ) always decreases the pressure drop in a positive way that can reduce a heavy reliance on a pump required to pressurize the working fluids before the working fluids enter the heat exchangers.The present study of using overlapping airfoils equipped within the flow channels of the PCHEs showed sensible promise in keeping heat transfer efficiency at a high level and reducing pressure drop for systems.However, more future research is required to optimize the design and performance of overlapping more than two airfoils and carry out experimental studies that can provide valuable insights for practical implementation.This research will determine if further heat transfer enhancement can be achieved for a wide range of practical problems.

F I G U R E 1
Upper panel of a schematic diagram illustrating the flow past a curved mass, which is represented by the top part of a single airfoil.Lower panel showing the flow past a bulge made of the trailing portion of an airfoil overlapping the head of another airfoil.The overlap forms an eddy pattern near a low-lying area, increasing the flow mixing, as represented by red eddy structures, and boosting heat transfer.

F I G U R E 2
Heat exchanger configurations with two cold plates (blue) and one hot plate (red) between the two cold plates and the operating conditions and dimensions considered herein: (A) solid exterior of the computational domain; (B) PCHE with no airfoil overlap of the chord length (L fin ); (C) PCHE with airfoil overlap of the entire length L fino through optimization with three geometric design parameters, that is., front length (f l ), overlap length (O l ), and back length (b l ), for an overlap base design of f l and b l both equal to 0.2 mm, and O l of 0.1 mm.

F I G U R E 3
Flow streamtrace patterns and velocity contours around a fin with two overlapping airfoils for the laminar (left panel) and shear stress transport (SST) turbulent models (right panel).

F I G U R E 5
NSGA-II flow chart used in the present study, in which the two shaded areas represent the COMSOL Multiphysics simulation framework with the three input design parameters and two output parameters.F I G U R E 6 NSGA-II associated with nondominated and crowding distance sorting.

F I G U R E 7
Total heat flux and velocity contours and streamline patterns of separate and overlapping airfoils.F I G U R E 8 Total heat flux decomposed into convective and conductive heat fluxes for separate and overlapping airfoils.

F I G U R E 9
Velocity contours of overlapping airfoils, with Re = 1088 (left) and Re = 4000 (right).

F I G U R E 11
Effectiveness of the present heat exchanger configuration with overlapping airfoils, as shown in Figure 2C, as a function of the geometric parameters of front length (f l ) and back length (b l ) for different values of overlap length (O l ).The asterisks indicate the locations of the maximum effectiveness.

F I G U R E 12
The correlation between pressure drop and Reynolds number for both the overlapping airfoils and separate airfoils considered in the present study, using a base design of f l and b l both equal to 0.2 mm, and O l of 0.1 mm for the overlapping airfoils.

F I G U R E 13
Pressure drop of the present heat exchanger configuration with overlapping airfoils, as shown in Figure 2C, as a function of the geometric parameters of front length (f l ) and back length (b l ) for different values of overlap length (O l ).The asterisks indicate the locations of the minimum pressure drop.
through the present NSGA-II optimization, when maximizing heat transfer enhancement, the use of overlapping airfoils of f l = 0.295, b l = 0.292, and O l = 0.222 mm in PCHEs results in a 4.1% improvement in F I G U R E 14 The Pareto front obtained through the present NSGA-II optimization.TA B L E 3 Results taken from a set of nondominated solutions (Pareto front) being chosen as optimal through the present NSGA-II optimization.Front length f l (mm) Overlap length O l (mm) Back length b l (mm) Effectiveness  Pressure drop P (Pa)

4 .
When the main goal was to increase heat transfer as much as possible via NSGA-II optimization, using overlapping airfoils in PCHEs of f l = 0.295, b l = 0.292, and O l = 0.222 mm led to a 4.1% improvement in effectiveness and an 82% increase in pressure drop relative to the utilization of an overlap base design (ε = 0.862 and Δp = 78,152 Pa) of the airfoils with f l = b l = 0.2 and O l = 0.1 mm. 5.When the main goal was to reduce pressure loss through NSGA-II optimization, utilizing overlapping airfoils of f l = 0.1, b l = 0.119 and O l = 0.203 mm resulted in more or less the same effectiveness and a 22% decrease in pressure drop relative to the base design mentioned in the previous point.6.The preferred design is achieved with f l = 0.1, b l = 0.119, and O l = 0.203 mm, meaning that such a judicious choice can maintain the effectiveness and have a minimum of the pressure drop (Δp = 61,000 Pa).
fluid density [kg/m 3 ] v fluid velocity [m/s] p fluid static pressure [Pa] K viscous stress tensor [N m −2 ]T f fluid temperature [K] T s solid temperature [K] k f fluid thermal conductivity [W m −1 k −1 ] k s solid thermal conductivity [W m −1 k −1 ] c p fluid specific heat at constant pressure [J kg −1 K −1 ] μ fluid viscosity [Pa s] μ t turbulent viscosity [Pa s] P k effective rate of production of k [kg m −1 s −3 ] P ω effective rate of production of ω [kg m −1 s −3 ] f v1 blending function [−] k turbulence kinetic energy [m 2 s −2 ] ω dissipation rate [m 2 s −3 ]  kinematic viscosity [m 2 s −1 ] ṁ mass flow rate [kg s -1 ] heat transfer rate [Watt] Q max maximum heat transfer rate [Watt] f l front length [mm] b l back length [mm] O l overlap length [mm]  effectiveness [−] P n original parentpopulation [−] P n+1 new parent population [−] O n offspring population [−] f 1 first objective used for Pareto front [−] f 2 second objective used for Pareto front [−] GREEK LETTERS  fluid density [kg/m 3 ] μ fluid viscosity [Pa s] μ t turbulent viscosity [Pa s] P k effective rate of production of k [kg m −1 s −3 ] P ω effective rate of production of ω [kg m −1 s −4 ] ω dissipation rate [m 2 s −3 ]  kinematic viscosity [m 2 s −1 ] σ ω constant [−] σ k constant [−] σ ω2 constant [−] minimum fin fin without airfoil overlap fino fin with airfoil overlap 37put parameters used for numerical validation through the experimental data taken from Liu et al.37 TA B L E 1 ,C, respectively, if not otherwise specified.The thermal conductivity and specific heat capacity of the solid domain made of UNS S30415 are 21 W/mK and 480 J/kgK, respectively.It is worth mentioning that the Reynolds numbers depicted in Figures10 and 12, spanning from 1088 to 5000, correspond to the mass flow rates listed in Table2.For Reynolds number calculation, the fluid dynamic viscosity is 7.98 × 10 −4 (Pa s) at 30 • C and 3.15 × 10 −4 (Pa s) at 90 • C and the fluid density is 995.71 kg∕m 3 at 30 • C and 965.06 kg∕m 3 at 90 • C, respectively, and the characteristic length of the airfoils is 4.85 × 10 −4 m, whose definition can be found in Yang et al.
Reynolds Numbers and corresponding mass flow rates.