Optimization of a Wavy Microchannel Heat Sink with Grooves

: In this study, a wavy microchannel heat sink with grooves using water as the working ﬂuid is proposed for application to cooling microprocessors. The geometry of the heat sink was optimized to improve heat transfer and pressure loss simultaneously. To achieve optimization goals, the average friction factor and thermal resistance were used as the objective functions. Three dimensionless parameters were selected as design variables: the distance between staggered grooves, groove width, and groove depth. A modiﬁed Latin hypercube sampling (LHS) method that combines the advantages of conventional LHS and a three-level full factorial method is also proposed. Response surface approximation was used to construct surrogate models, and Pareto-optimal solutions were obtained with a multi-objective genetic algorithm. The modiﬁed LHS was proven to have better performance than the conventional LHS and full factorial methods in the present optimization problem. A representative optimal design showed that both the thermal resistance and friction factor improved by 1.55% and 3.00%, compared to a reference design, respectively.


Introduction
Microprocessors generate high heat flux and thermally interact with their surroundings. Because they are composed of many integrated components, their efficiency and performance are significantly influenced by temperature. The temperature must be kept between 363 K and 383 K to maintain the best performance [1]. Therefore, it is essential to develop an effective cooling system to maintain a proper temperature even with high heat generation.
Cooling systems are being made less noisy and smaller, and it is estimated that heat sinks capable of cooling at more than 1000 W/cm 2 will be required in the near future [2]. Both air and water can be used for cooling systems, but as the heat generation increases with the development of microprocessors, air cooling systems have a limitation in maintaining effective cooling performance [1]. In addition, to increase air-cooling performance, the fan speed must be increased, which increases noise.
Water-cooling microchannel heat sinks have been widely used to reduce noise and meet increased requirements for cooling. Numerical and experimental studies on these heat sinks have been actively conducted. Tuckerman and Pease [3] experimented with a microchannel heat sink consisting of straight channels with heat flux of 790 W/cm 2 using water as a coolant. They confirmed that the water had great heat transfer characteristics. Wang et al. [4] carried out experiments and numerical analyses on a microchannel heat sink consisting of straight channels with ribs and grooves. They found that secondary flows occurring behind the ribs and grooves prevented the formation of the thermal boundary layer and promoted fluid mixing and heat transfer.
Ansari et al. [5] and Farhanieh et al. [6] investigated the cooling performance of straight microchannels with grooves. They found that the interfacial area of heat transfer is increased by the grooves, thereby increasing the cooling performance. Ansari et al. [5] suggested that high local Nusselt numbers are obtained near the upstream and downstream regions of the groove structures. Farhanieh et al. [6] confirmed that although the heat transfer performance was low due to recirculating flow in the grooved area, the groove structure prevented the formation of the thermal boundary layer, enhancing the overall heat transfer performance.
Greiner et al. [7] evaluated the friction coefficient and cooling performance of flow paths with triangular grooves through experimental and numerical analyses. In the case of laminar flow, the friction coefficient decreased as the hydraulic diameter was increased by the grooves. As the Reynolds number increased, the flow over the triangular grooves became more complex. Consequently the heat transfer performance improved. Xie et al. [8] confirmed that in the case of grooved channels, the heat transfer performance improved in the area, where the fluid velocity increased due to the narrowed channel.
Recently, some studies were performed on curved microchannel heat sinks [9,10]. Gong et al. [9] proposed wavy microchannels for a heat sink and studied its cooling performance. They found that the wavy microchannel heat sink caused a vortex at the trough and crest sections in each cycle, resulting in a significant improvement in cooling performance compared to a straight-microchannel heat sink. The pressure losses did not increase significantly. Sui et al. [10] investigated a wavy microchannel heat sink for various flow conditions and amplitudes. They compared three-dimensional numerical analyses and experiments in all cases. The numerical results of the cooling performance and friction coefficient were reliable, and vortices were developed at the trough and crest sections in each cycle.
With the rapid development of computing power, optimization designs that can handle a huge amount of data have become practical [11][12][13]. In particular, optimization based on a surrogate model has been widely used to reduce the computing time [14][15][16]. A surrogate model is constructed using sample data at several selected points in the design space. Design of experiments (DOE) is used to extract the sample data. The prediction accuracy of surrogate models is affected by the sample data, so DOE should be carried out carefully [17]. DOE can typically be classified into two categories according to the extraction method: factorial design [18] based on orthogonal extraction and Latin hypercube sampling (LHS) [19], which uses random extraction.
Factorial designs are experiments that combine all levels of each factor with all levels of all other factors in an experiment [18]. This method is very intuitive, and the number of sample data is determined by the number of design variables. Therefore, it is easy to use because it can be applied without considering the distribution and number of sample data. However, the features within the design space are considered less because the sample data are focused on the boundaries of the design space.
LHS uses a random extraction method for sample points within the design space. This method is widely used because the space-filling quality of the sampling points is good. In addition, it can create any allotted number of sampling points [19]. However, since LHS extracts sample data inside the design space, the boundary values of each design variable are not considered. Therefore, a surrogate model based on LHS often makes predictions that are too high at the bounds of design variables.
In the present work, a modified LHS that uses the advantages of orthogonal and random extraction methods is proposed. The modified LHS extracts sample data by applying LHS inside the design space and prevents excessive prediction of the surrogate model by applying the three-level full factorial method to the boundaries of the design space. A wavy microchannel heat sink improves the heat transfer efficiency compared with straight microchannels. In this study, a wavy microchannel heat sink with grooves attached to the channel walls is proposed. A numerical analysis of the laminar flow and heat transfer in the microchannels was performed using three-dimensional Navier-Stokes equations.
Multi-objective optimization of the wavy microchannels was also performed to simultaneously enhance the heat transfer efficiency and reduce the pressure loss. The proposed modified LHS was compared with conventional DOE methods to determine the effec- tiveness of the proposed sampling method in constructing surrogate models. For the optimization, response surface approximation (RSA) [20] and a genetic algorithm [21] were used as a surrogate model and searching algorithm, respectively.

Numerical Analysis
The computational domain and design variables are shown in Figure 1. The heat sink consists of 62 wavy microchannels, and the computational domain includes one of them, which is composed of 10 cycles. Two grooves are attached to each channel wall in one cycle of the channel, as shown in Figure 1. The amplitude of the wavy channel is 138 µm, and the total length (20P) is 25,000 µm. The thickness of the side wall (2t) is 193 µm, the channel width (W) is 207 µm, and the height of the flow path (h) is 406 µm.
ultaneously enhance the heat transfer efficiency and reduce the pressure loss. The proposed modified LHS was compared with conventional DOE methods to determine the effectiveness of the proposed sampling method in constructing surrogate models. For the optimization, response surface approximation (RSA) [20] and a genetic algorithm [21] were used as a surrogate model and searching algorithm, respectively.

Numerical Analysis
The computational domain and design variables are shown in Figure 1. The heat sink consists of 62 wavy microchannels, and the computational domain includes one of them, which is composed of 10 cycles. Two grooves are attached to each channel wall in one cycle of the channel, as shown in Figure 1. The amplitude of the wavy channel is 138 μm, and the total length (20P) is 25,000 μm. The thickness of the side wall (2t) is 193 μm, the channel width (W) is 207 μm, and the height of the flow path (h) is 406 μm.
The information for the reference channel is shown in Table 1. In the reference channel, the horizontal locations of grooves on both the walls are the same (D = 0). Grooves are attached to crests and troughs when D = 0. As D increases, grooves on the left wall (located at z = 400 μm in Figure 1) move in the -x direction, and grooves on the right wall (located at z = 0 in Figure 1) move in the +x direction from a crest (or a trough).   The information for the reference channel is shown in Table 1. In the reference channel, the horizontal locations of grooves on both the walls are the same (D = 0). Grooves are attached to crests and troughs when D = 0. As D increases, grooves on the left wall (located at z = 400 µm in Figure 1) move in the −x direction, and grooves on the right wall (located at z = 0 in Figure 1) move in the +x direction from a crest (or a trough). Conjugate heat transfer analysis was carried out on the flow channel and solid domain for laminar flow in steady state using three-dimensional Navier-Stokes equations. The commercial computational fluid dynamics code ANSYS CFX 15.0 ® (Version 15.0, ANSYS Inc., Canonsburg, PA, USA, 2013) [22] was used for the analysis. The boundary conditions are shown in Figure 2. The boundary conditions were used in the same way as in a previous study [10]. The working fluid was water at 300 K, and the material of the solid wall was copper.
Processes 2021, 9, x FOR PEER REVIEW 4 of 20 Conjugate heat transfer analysis was carried out on the flow channel and solid domain for laminar flow in steady state using three-dimensional Navier-Stokes equations. The commercial computational fluid dynamics code ANSYS CFX 15.0 ® (Version 15.0, AN-SYS Inc., Canonsburg, PA, USA, 2013) [22] was used for the analysis. The boundary conditions are shown in Figure 2. The boundary conditions were used in the same way as in a previous study [10]. The working fluid was water at 300 K, and the material of the solid wall was copper. The inlet Reynolds number ( = /μ) was 700, and the corresponding flow rate was assigned at the inlet. The average velocity of water at the inlet is 1.992 m/s. The static pressure was used as an exit boundary condition. Periodic conditions for temperature were applied to both side boundaries (i.e., wavy surfaces) considering the neighboring microchannels. Uniform heat flux conditions (50 W/cm 2 ) were applied to the bottom boundary (substrate), and the upper boundary was assumed to be adiabatic.
The fluid and solid computational domains consist of hexahedral and unstructured tetrahedral meshes, respectively. Since the flow velocity changes near the groove, dense meshes are placed there. Near the solid wall, dense meshes were also placed in anticipation of large temperature and velocity gradients due to the boundary layer. Figure 3 shows an example of the grid system. Convergence conditions were set so that the root-meansquared residual values of all parameters fell to 1.00 × 10 -6 . Each computation took about 5-6 h on a computer with an Intel Core i7-4790K 4GHz CPU. The inlet Reynolds number (Re = ρUD h /µ) was 700, and the corresponding flow rate was assigned at the inlet. The average velocity of water at the inlet is 1.992 m/s. The static pressure was used as an exit boundary condition. Periodic conditions for temperature were applied to both side boundaries (i.e., wavy surfaces) considering the neighboring microchannels. Uniform heat flux conditions (50 W/cm 2 ) were applied to the bottom boundary (substrate), and the upper boundary was assumed to be adiabatic.
The fluid and solid computational domains consist of hexahedral and unstructured tetrahedral meshes, respectively. Since the flow velocity changes near the groove, dense meshes are placed there. Near the solid wall, dense meshes were also placed in anticipation of large temperature and velocity gradients due to the boundary layer. Figure 3 shows an example of the grid system. Convergence conditions were set so that the root-meansquared residual values of all parameters fell to 1.00 × 10 −6 . Each computation took about 5-6 h on a computer with an Intel Core i7-4790K 4GHz CPU.

Optimization Procedure
The multi-objective optimization problem was formulated as follows:

Optimization Procedure
The multi-objective optimization problem was formulated as follows: is the vector of real-valued objective functions, x is a vector of the design variables, and LB and UB indicate the vectors of the lower and upper bounds, respectively [23]. Figure 4 shows a flowchart for the multi-objective genetic algorithm (MOGA) optimization process using a surrogate model. Firstly, the objective functions and constraints are defined according to the design goals. Secondly, the design variables and their ranges are selected.

Optimization Procedure
The multi-objective optimization problem was formulated as follows: is the vector of real-valued objective functions, x is a vector of the design variables, and LB and UB indicate the vectors of the lower and upper bounds, respectively [23]. Figure 4 shows a flowchart for the multi-objective genetic algorithm (MOGA) optimization process using a surrogate model. Firstly, the objective functions and constraints are defined according to the design goals. Secondly, the design variables and their ranges are selected. The full factorial method, LHS, and modified LHS were used in DOE to select design points (i.e., sample designs). Values of the objective functions were evaluated by numerical analysis at these design points. Next, to approximate the objective functions, surrogate The full factorial method, LHS, and modified LHS were used in DOE to select design points (i.e., sample designs). Values of the objective functions were evaluated by numerical analysis at these design points. Next, to approximate the objective functions, surrogate models were constructed using these objective function values. A genetic algorithm (GA) was used to find the global optima. Finally, Pareto-optimal solutions (a collection of nondominated solutions) were derived using MATLAB (Release 14, the Math Work Inc., Natick, MA, USA, 2004) [24].

Design Variables and Objective Functions
For optimization, three geometric parameters were selected as the design variables through a preliminary parametric study: the ratio of the groove depth to half of the side wall thickness (d/t), the ratio of the groove width to half of the cycle length (w/P), and the distance between staggered grooves on the opposite walls to half of the cycle length (D/P). The depth and width of the grooves were expected to affect the vortices occurring around the grooves and thus have sensitive effects on the heat transfer. The distance between staggered grooves affects the disturbance of the main flow.
The average Nusselt number Nu is defined as follows [10]: where D h is the hydraulic diameter of the microchannel, and k w is the thermal conductivity of water. h is the average convective heat transfer coefficient, which is defined as follows: where q is the heat flux, A b and A c are the bottom area and the side area of the flow channel, and T w and T m are the average temperature at the solid wall and the average of the inlet and outlet temperatures, respectively.
The friction factor f is defined as follows [10]: where dp, ρ, and U are the pressure difference between the inlet and outlet, the density of water, and the average velocity at the inlet, respectively. The thermal resistance R th is defined as follows [5]: where T s,max and T f,inlet are the highest temperature at the bottom substrate and the average temperature of the cooling fluid at the inlet, respectively, and A is the area of the microchannel substrate. The local Nusselt number (Nu x ) is defined as follows: where q l and T w,l are the local heat flux and local temperature at the surface of the solid wall, respectively. R th and f were selected as objective functions for the multi-objective optimization: F Rth = R th and F f = f. The thermal resistance R th is related to the highest local temperature, which affects the performance of micro devices. The friction factor f was used to reduce the pressure drop through the microchannel. A parametric study was carried out for the performance functions using three design variables: D/P, w/P, and d/t. Based on the parametric study, the ranges of the three design variables were selected, as shown in Table 2.

Modified LHS
Factorial design [18,25] is a classical DOE method that explores the design space. 2-level and 3-level full factorial designs are widely used to estimate interactions between design variables. Figure 5 shows examples of 2-and 3-level full factorial designs for two design variables. In the 2-level full factorial design, the sample points are located at the ends of each boundary, as shown in Figure 5a. In the 3-level full factorial design, the sample points are located at the ends and middle of each boundary, as shown in Figure 5b. In these full factorial designs, the distribution and number of the sample points are determined according to the number of design variables when the level is determined. level and 3-level full factorial designs are widely used to estimate interactions between design variables. Figure 5 shows examples of 2-and 3-level full factorial designs for two design variables. In the 2-level full factorial design, the sample points are located at the ends of each boundary, as shown in Figure 5a. In the 3-level full factorial design, the sample points are located at the ends and middle of each boundary, as shown in Figure 5b. In these full factorial designs, the distribution and number of the sample points are determined according to the number of design variables when the level is determined. LHS [19] is one of the most popular DOE methods for random sample distributions. To allocate p samples using LHS, the range of each parameter is separated into p bins, which yields a total number of p n bins for n design variables in the design space. The samples are randomly chosen in the design space, each sample is randomly arranged inside a bin, and for all one-dimensional projections of the p samples and bins, there is exactly one sample in each bin, as shown in Figure 6. Therefore, LHS is relatively incapable of handling samples at the boundaries of the design space compared to the full factorial designs.
Since the surrogate model is built using the data at the sample points, the distribution of the sample points has a very significant influence on the prediction accuracy of the surrogate model. Therefore, when using LHS with sample points concentrated inside the design space, it is possible to predict the interaction well inside the design space, but predictions that are too high may occur at the boundaries where there are no data. On the other hand, in the full factorial design, the sample points are focused on the boundaries of the design space, so it is possible to make a relatively accurate prediction at the boundaries, but there is a problem in the prediction inside the design space.
To solve this problem, a modified LHS is proposed. In the modified LHS, sample points are extracted using the 2-level full factorial method at the boundaries of the design space, and the LHS method is used to select sample points inside the design space. An example of the modified LHS for two design variables is shown in Figure 7. In this method, the surrogate model does not show high predictions at the boundaries of the design space. Furthermore, by selecting the sample points inside the design space, the shortcomings of the full factorial design can be overcome. MATLAB [24] was used to extract the sample points. Three-level full factorial design, LHS, and the modified LHS were tested for the same optimization problem. Twenty-seven sample points were extracted for all these DOE methods. LHS [19] is one of the most popular DOE methods for random sample distributions. To allocate p samples using LHS, the range of each parameter is separated into p bins, which yields a total number of p n bins for n design variables in the design space. The samples are randomly chosen in the design space, each sample is randomly arranged inside a bin, and for all one-dimensional projections of the p samples and bins, there is exactly one sample in each bin, as shown in Figure 6. Therefore, LHS is relatively incapable of handling samples at the boundaries of the design space compared to the full factorial designs.   Since the surrogate model is built using the data at the sample points, the distribution of the sample points has a very significant influence on the prediction accuracy of the surrogate model. Therefore, when using LHS with sample points concentrated inside the design space, it is possible to predict the interaction well inside the design space, but predictions that are too high may occur at the boundaries where there are no data. On the other hand, in the full factorial design, the sample points are focused on the boundaries of the design space, so it is possible to make a relatively accurate prediction at the boundaries, but there is a problem in the prediction inside the design space.
To solve this problem, a modified LHS is proposed. In the modified LHS, sample points are extracted using the 2-level full factorial method at the boundaries of the design space, and the LHS method is used to select sample points inside the design space. An example of the modified LHS for two design variables is shown in Figure 7. In this method, the surrogate model does not show high predictions at the boundaries of the design space. Furthermore, by selecting the sample points inside the design space, the shortcomings of the full factorial design can be overcome. MATLAB [24] was used to extract the sample points. Three-level full factorial design, LHS, and the modified LHS were tested for the same optimization problem. Twenty-seven sample points were extracted for all these DOE methods.

Surrogate Model and Searching Algorithm
The surrogate model was configured based on the sample points obtained using DOE methods. Response surface approximation (RSA) [20] was used as the surrogate model. MOGA coupled with the RSA model was used to obtain Pareto-optimal solutions [24].
The RSA model is multivariate polynomial model, and a continuous response y(x) is usually modeled as follows [20]: where x is a vector of design variables, fj(x) (j = 1, …, N) are the terms of the model, βj (j = 1, …, N) are the coefficients, and the error ε is assumed to be uncorrelated and distributed with a mean of 0 and constant variance [20]. A second-order polynomial is used for the RSA model, and the model can be expressed as follows:

Surrogate Model and Searching Algorithm
The surrogate model was configured based on the sample points obtained using DOE methods. Response surface approximation (RSA) [20] was used as the surrogate model. MOGA coupled with the RSA model was used to obtain Pareto-optimal solutions [24].
The RSA model is multivariate polynomial model, and a continuous response y(x) is usually modeled as follows [20]: where x is a vector of design variables, f j (x) (j = 1, . . . , N) are the terms of the model, β j (j = 1, . . . , N) are the coefficients, and the error ε is assumed to be uncorrelated and distributed with a mean of 0 and constant variance [20]. A second-order polynomial is used for the RSA model, and the model can be expressed as follows: The model involves an intercept, linear terms, quadratic interaction terms, and squared terms (from left to right). R 2 and R adj 2 are used to decide the goodness of the fit and should be close to 1 for a good fit [20].
GA is a random global search technique that solves problems based on natural evolution. An initial population of individuals is defined to represent a part of the solution to a problem [21]. Before starting the search, a set of chromosomes is randomly selected from the design space to obtain the initial population. Through subsequent computations, the individuals adapt in a competitive way. The initially selected set of chromosomes is called the parental generation, and the subsequent selected set of chromosomes is called the child generation. In this process, genetic search operators (selection, mutation, and crossover) are used to obtain chromosomes that are superior to the previous generation [21]. MATLAB [24] was used to invoke the GA for multi-objective optimization.

Grid Dependency Test and Validation of Numerical Results
A grid dependency test was carried out for the reference shape based on Richardson's extrapolation method and grid convergence index (GCI), which represents numerical uncertainty by estimating the discretization error according to the procedure presented by Roache [26] and Celik and Karatekin [27]. Table 3 shows the results of calculating the discretization error for Nu. The number of grid nodes was adjusted by setting the grid segmentation index to 1.3, and three different grid systems were analyzed. When N 2 was used, the extrapolation error (e 21 ext ) was about 0.3%, and GCI 21 f ine was about 0.4%, which indicate small numerical uncertainty. Therefore, N 2 was selected as the optimal grid system. To verify the numerical results, they were compared with experimental data obtained by Sui et al. [10] for the Nusselt number and friction factor in a wavy microchannel under the same boundary conditions, as shown in Figure 8. As shown in Figure 8, the numerical results for the friction factor show good agreement with the experimental data, except at the lowest Reynolds number. The numerical results for the Nusselt number deviate slightly from the experimental data over the whole Re range but show the same qualitative tendency. At Re = 300, the errors are relatively large because the pressure drop and the flow rate are relatively small, as discussed by Sui et al. [10].
Processes 2021, 9, x FOR PEER REVIEW 10 of 20 about 8.34%, and Rth decreased by about 2%, but the friction factor f also increased by about 1.25% in comparison with the smooth wavy channel. This means that the grooves largely enhance the heat transfer but with less increase in the friction. This improvement in the heat sink performance with grooves is expected to be further increased by optimization.    Table 4 shows the comparison of the performance parameters between the smooth and reference wavy microchannels. In the wavy channel with grooves, Nu increased by about 8.34%, and R th decreased by about 2%, but the friction factor f also increased by about 1.25% in comparison with the smooth wavy channel. This means that the grooves largely enhance the heat transfer but with less increase in the friction. This improvement in the heat sink performance with grooves is expected to be further increased by optimization. The temperature distributions on the wavy wall on the right side of the reference and smooth microchannels are shown in Figure 9. In the case of the reference design, it can be seen that the temperature increase in the flow direction is smaller than that of the smooth microchannel, resulting in lower maximum temperature. This shows improved heat transfer performance on the sidewalls and confirms the results shown in Table 4.   The local Nusselt number (Nu x ) distributions on the wavy side walls are shown in Figure 10. High Nu x regions are found between the crest and the trough (e.g., x/2P = 5.75-6.25) on the left wall, but they are found between a trough and crest (e.g., x/2P = 6.25-6.75) on the right wall. In addition, most of the high Nu x regions are distributed near the top and bottom of the flow path. In the case of the smooth microchannel, a low Nu x region is found near the middle height of the flow path immediately after each crest (e.g., x/2P = 5.75). In the reference design, the high Nu x regions are found just downstream of the grooves, and the total area of the high Nu x regions is larger than that of the smooth channel. This results in high Nu in the grooved microchannel, as shown in Table 4.    Figure 11 shows the flow fields of the reference design and the smooth wavy microchannel. Figure 11a shows that the velocity gradients near the left wall are larger than those near the right wall between a crest and trough (x/2P = 5.75-6.25), but vice versa between the trough and crest (x/2P = 6.25-6.75). This phenomenon occurs due to the fluid inertia. The regions with high velocity gradients and those with high Nu x shown in Figure 10 are nearly the same. Thus, it can be inferred that the high velocity gradient near the wall promotes the heat transfer and enhances Nu x .   Figure 11b shows the velocity vectors in the y-z plane. Vortices are found near the top and bottom of the flow channel in the crest. In the case of the reference design with grooves, a complicated flow structure is found near the left wavy wall at the edge of a groove (x/2P = 5.875) due to the upward flow escaping from the groove, which promotes mixing of the fluid (and thereby heat transfer) in these regions. This is due to sudden contraction of the flow area just downstream of a groove and provides a reason for the high Nux regions downstream of the grooves shown in Figure 10b. Even though the grooves are at the same locations on both the wavy walls, the flow fields shown in Figure  11b are not symmetric in the z direction because the main flow upstream of the groove proceeds in the +z direction. Figure 12 shows the temperature distributions in the y-z plane at the inflection point of the wavy channel (x/2P = 6). The temperature gradient is relatively small near the upper Figure 11. Velocity vectors on x-z plane and y-z plane (a) velocity vectors on the x-z plane (y/H = 0.6) (b) velocity vectors on the y-z plane. Figure 11b shows the velocity vectors in the y-z plane. Vortices are found near the top and bottom of the flow channel in the crest. In the case of the reference design with grooves, a complicated flow structure is found near the left wavy wall at the edge of a groove (x/2P = 5.875) due to the upward flow escaping from the groove, which promotes mixing of the fluid (and thereby heat transfer) in these regions. This is due to sudden contraction of the flow area just downstream of a groove and provides a reason for the high Nu x regions downstream of the grooves shown in Figure 10b. Even though the grooves are at the same locations on both the wavy walls, the flow fields shown in Figure 11b are not symmetric in the z direction because the main flow upstream of the groove proceeds in the +z direction. Figure 12 shows the temperature distributions in the y-z plane at the inflection point of the wavy channel (x/2P = 6). The temperature gradient is relatively small near the upper and lower sides of the flow path in common. This is thought to be due to the strong vortices shown in Figure 11b. These low temperature gradients also contribute to the distribution of high Nu x in these regions ( Figure 10). Figure 12 shows that the temperature on the left side in the reference design is still low, even at the medium height, unlike in the smooth wavy microchannel. This is due to the fluid mixing caused by the strong secondary flow downstream of the grooves. and lower sides of the flow path in common. This is thought to be due to the strong vortices shown in Figure 11b. These low temperature gradients also contribute to the distribution of high Nux in these regions ( Figure 10). Figure 12 shows that the temperature on the left side in the reference design is still low, even at the medium height, unlike in the smooth wavy microchannel. This is due to the fluid mixing caused by the strong secondary flow downstream of the grooves.  Figure 13 shows the results of the parametric study for , Rth, and f. When one parameter was changed, the other parameters were fixed at the reference values shown in Table 2. With the change of parameters, the friction factor f shows small variations of less  Figure 13 shows the results of the parametric study for Nu, R th , and f. When one parameter was changed, the other parameters were fixed at the reference values shown in Table 2. With the change of parameters, the friction factor f shows small variations of less than 2.5%. At d/t = 0.5, the maximum Nu and f and minimum R th are found, which indicates that there is an optimum groove depth for heat transfer enhancement.  Figure 13 shows the results of the parametric study for , Rth, and f. When one parameter was changed, the other parameters were fixed at the reference values shown in Table 2. With the change of parameters, the friction factor f shows small variations of less than 2.5%. At d/t = 0.5, the maximum and f and minimum Rth are found, which indicates that there is an optimum groove depth for heat transfer enhancement. R th and Nu are inversely correlated with the variation of d/t. Nu and f have maximum values at w/P = 0.15, but R th has a minimum value at w/P = 0.25. R th decreases as w/P increases in the tested range. This means that wide grooves are effective in reducing thermal resistance. In the case of D/P, Nu and f increase as D/P increases, but R th shows a maximum at D/P = 0 (non-staggered grooves). This indicates that if the absolute value of D is fixed, the relative locations of grooves between the two wavy walls do not affect R th unlike Nu and f.

Parametric Study
At D/P = 0.5, Nu shows the maximum value of 20.2. This is the value improved by 23.2% from the value (Nu = 16.4) of the wavy microchannel heat sink without grooves predicted by Sui et al. [10] under the same geometric and Reynolds number conditions. Figure 14 shows the distribution of sample points of three different DOEs. In Figure 14a, the 3-level full factorial design shows that the sample points are located on only the boundaries of the design space. In the case of the conventional LHS, the sample points are distributed in only the design space, as shown in Figure 14b. However, in the modified LHS, the sample points are located on the boundaries and inside of the design space, as shown in Figure 14c. The RSA models for the objective functions, F f and F Rth , were formulated in terms of the design variables normalized between 0 and 1 for three different DOEs as follows:  Figure 15 shows the Pareto-optimal solutions using three different DOEs. The Paretooptimal solutions are the best solutions that can be achieved for one objective without disadvantaging another objective and are sensitive to the constructed surrogate model [28]. Pareto-optimal fronts obtained using the conventional and modified LHS methods have similar smooth curves, as shown in Figure 15. However, the full factorial design shows a curve that has two inflection points, unlike the other curves. The Pareto-optimal front from modified LHS predicts the lowest optimum values of the two objective functions among the tested DOEs in most of the range. The Pareto-optimal fronts of the two LHS methods cover wider ranges than that of the full factorial design.

Optimization Results with Different DOEs
To compare the optimization results, three Pareto-optimal designs (PODs) were extracted from each Pareto-optimal front using K-means clustering [29], as presented in Figure 15. The predicted objective function values at the PODs and numerical calculations at the same PODs are compared in Table 5. In case of the full factorial design, the PODs are close to the boundaries of one or two design variables. However, the PODs of conventional LHS are found inside the design space. In this case, a POD close to a boundary (D/P = −0.8802 in POD A) yields a larger relative error between the predicted and calculated objective function values than the other PODs. As mentioned earlier, the surrogate model obtained using conventional LHS over-predicts the values at the boundary of the design space, and the error increases near the boundary. 14a, the 3-level full factorial design shows that the sample points are located on only the boundaries of the design space. In the case of the conventional LHS, the sample points are distributed in only the design space, as shown in Figure 14b. However, in the modified LHS, the sample points are located on the boundaries and inside of the design space, as shown in Figure 14c. The RSA models for the objective functions, Ff and FRth, were formulated in terms of the design variables normalized between 0 and 1 for three different DOEs as follows:  To compare the optimization results, three Pareto-optimal designs (PODs) were extracted from each Pareto-optimal front using K-means clustering [29], as presented in Figure 15. The predicted objective function values at the PODs and numerical calculations at the same PODs are compared in Table 5. In case of the full factorial design, the PODs are close to the boundaries of one or two design variables. However, the PODs of conventional LHS are found inside the design space. In this case, a POD close to a boundary (D/P = −0.8802 in POD A) yields a larger relative error between the predicted and calculated objective function values than the other PODs. As mentioned earlier, the surrogate model obtained using conventional LHS over-predicts the values at the boundary of the design space, and the error increases near the boundary.  In the case of the modified LHS, even the POD located close to the boundary shows good prediction with maximum relative error less than 1.5%. Thus, the modified LHS shows the best prediction accuracy among the tested DOEs. The full factorial design and conventional LHS generally over-predict the objective function values with positive relative errors at the three PODs, but the modified LHS generally under-predicts the values. Therefore, there is not much difference in the calculated objective function values at the PODs among the tested DOEs. In the case of the modified LHS, even the POD located close to the boundary shows good prediction with maximum relative error less than 1.5%. Thus, the modified LHS shows the best prediction accuracy among the tested DOEs. The full factorial design and conventional LHS generally over-predict the objective function values with positive relative errors at the three PODs, but the modified LHS generally under-predicts the values. Therefore, there is not much difference in the calculated objective function values at the PODs among the tested DOEs. The R 2 and adjusted R 2 values of the RSA models constructed using three different DOEs are listed in Table 6. As mentioned earlier, values closer to 1 indicate a better surrogate model. In this respect, the modified LHS shows the best results in Table 6. This is consistent with the results shown in Table 5. In addition, the RSA model with the full factorial design has the worst performance. Applying a DOE method that properly locates the sample points on the boundaries and inside of the design space makes it possible to construct a more accurate surrogate model and obtain superior optimization results.

Analysis of the Optimized Design
POD 2 obtained with the modified LHS was selected for further analysis because the values of both the objective functions were improved compared to the reference design. In POD 2, R th and f are reduced by 1.55% and 3.00%, respectively, compared with those of the reference design. The Nu x distributions on the wavy walls are shown in Figure 16, which compares the heat transfer performance between POD 2 and the reference design. In the case of the reference design, there are high Nu x regions on the left wall between the grooves in the crest and the trough (x/2P = 5.75-6.25). In the case of POD 2, high Nu x regions are shown on the left side wall between the inflection point where the groove is located and the trough (x/2P = 5.50-6.25). The right wall between x/2P = 7.00 and 7.50 shows the same Nu x distribution as the left wall between x/2P = 5.50 and 6.25. Unlike the reference design, the high Nu x distribution for POD 2 from the inflection point to the crest (x/2P = 5.50-5.75) seems to be one of the important factors in reducing R th compared to the reference design.  Figure 17 shows the streamlines and velocity vectors. Figure 17a shows that the recirculating flow develops in the grooves in the reference design. These recirculating flow regions correspond to the low Nux regions in the grooves in Figure 16. This occurs because the flow recirculation hinders the heat transfer on the wall. Figure 17b shows the velocity vectors at the cross sections perpendicular to the flow direction. As mentioned earlier, in the reference design, a strong interaction between the vortical flow in the channel and upward flow from the groove is found near the left wall at the edge of the groove (x/2P = 5.875). In the case of POD 2, this phenomenon is also found near the left wavy wall at the edge of the groove (x/2P = 5.593) but is weaker than in the reference design. This seems to be affected by the location of the groove (which is the crest in the reference design but an inflection point in POD 2) and the fact that the two grooves on both the walls are attached at the same location in the reference design.  Figure 17 shows the streamlines and velocity vectors. Figure 17a shows that the recirculating flow develops in the grooves in the reference design. These recirculating flow regions correspond to the low Nu x regions in the grooves in Figure 16. This occurs because the flow recirculation hinders the heat transfer on the wall. Figure 17b shows the velocity vectors at the cross sections perpendicular to the flow direction. As mentioned earlier, in the reference design, a strong interaction between the vortical flow in the channel and upward flow from the groove is found near the left wall at the edge of the groove (x/2P = 5.875). In the case of POD 2, this phenomenon is also found near the left wavy wall at the edge of the groove (x/2P = 5.593) but is weaker than in the reference design. This seems to be affected by the location of the groove (which is the crest in the reference design but an inflection point in POD 2) and the fact that the two grooves on both the walls are attached at the same location in the reference design.
However, in POD 2, the vortices in the channel become stronger at a location downstream (x/2P = 5.875). This is consistent with the Nu x distributions shown in Figure 16. In the reference design, the high Nu x region persists from the groove edge (x/2P = 5.875) to a location far downstream but disappears before the next groove. However, in POD 2, the high Nu x region is relatively narrow at the edge of the groove (x/2P = 5.593) but grows downstream and becomes widest at the upstream edge of the next groove. Thus, the total area of the high Nu x regions is larger in POD 2 than that in the reference design.
the reference design, a strong interaction between the vortical flow in the channel and upward flow from the groove is found near the left wall at the edge of the groove (x/2P = 5.875). In the case of POD 2, this phenomenon is also found near the left wavy wall at the edge of the groove (x/2P = 5.593) but is weaker than in the reference design. This seems to be affected by the location of the groove (which is the crest in the reference design but an inflection point in POD 2) and the fact that the two grooves on both the walls are attached at the same location in the reference design.

Conclusions
A wavy microchannel heat sink with grooves was optimized using RANS analysis of the flow and conjugate heat transfer. In the reference wavy channel with grooves, Nu increased by about 8.34%, and R th decreased by about 2%, but the friction factor f also increased by about 1.25% compared to the smooth wavy channel. Thus, using grooves, the enhancement of the heat transfer surpasses the increase in the friction.
For optimization, the distance between staggered grooves on opposite wavy walls, the groove depth, and the groove width were selected as design variables. The thermal resistance (R th ) and friction factor (f ) were used as objective functions. A modified LHS that uses the advantages of conventional LHS and the three-level full factorial method was also proposed. The optimization performance of three DOE methods was estimated. Surrogate models of the objective functions were constructed by RSA with each DOE method. The corresponding Pareto optimal solutions were derived, and three representative optimal solutions were selected to compare the predictions of the DOE methods.
The results showed that the optimal solutions using modified LHS methods have the best predictions with less than 1.5% error compared to the numerical calculations. They also had the largest R 2 and adjusted R 2 values, which indicate the best statistical accuracy of the RSA models. For one of the representative optimal solutions, POD2, R th and f decreased by 1.55% and 3.00%, respectively, compared to the reference design, indicating that both the objective functions were improved. Therefore, the multi-objective optimization with modified LHS could effectively improve the performance of the wavy microchannel heat sink with grooves.