Optimization of Microjet Location Using Surrogate Model Coupled with Particle Swarm Optimization Algorithm

Abstract

This study aimed to present the design methodology of microjet heat sinks with unequal jet spacing, using a machine learning technique which alleviates hot spots in heat sinks with non-uniform heat flux conditions. Latin hypercube sampling was used to obtain 30 design sample points on which three-dimensional Computational Fluid Dynamics (CFD) solutions were calculated, which were used to train the machine learning model. Radial Basis Neural Network (RBNN) was used as a surrogate model coupled with Particle Swarm Optimization (PSO) to obtain the optimized location of jets. The RBNN provides continuous space for searching the optimum values. At the predicted optimum values from the coupled model, the CFD solution was calculated for comparison. The percentage error for the target function was 0.56%, whereas for the accompanied function it was 1.3%. The coupled algorithm has variable inputs at user discretion, including gaussian spread, number of search particles, and number of iterations. The sensitivity of each variable was obtained. Analysis of Variance (ANOVA) was performed to investigate the effect of the input variable on thermal resistance. ANOVA results revealed that gaussian spread is the dominant variable affecting the thermal resistance.

1. Introduction

A high-heat flux removal from microdevices is still a major challenge even after more than four decades of research and development. The micro-cooling devices are being investigated both experimentally and in combination with the computational study. The experimental results are not synchronous, but the computational results appear to minimize the anomaly in the experimental analysis of micro-devices termed as “scaling effects,” which are difficult to capture in experimental investigation [1]. Various computation-intensive studies have been performed in the past exploring the effect of fluid flow characteristics and the heat-transfer effects. The micro-cooling devices encompass the use of microchannels [2,3,4,5], porous materials [6], and microjets [7,8,9,10,11,12,13,14] with a promising future in all types of designs. However, each of them has its own set of challenges, from the design stage to the manufacturing stage. The substrate temperature non-uniformity along the length of the microchannel heat sink is one of its inherent features. In order to counter it, Samal and Moharana [15] provided a novel concept of recharging design of microchannel, which improves temperature uniformity and reduces back axial conduction. The porous materials have to tackle the challenge of the high-pressure drop [6]. The microjet impingement needs to deal with the problem of material erosion [16], but relaxes on the constraint of temperature non-uniformity and pressure drop. Most of the investigation of micro-cooling assumes uniform heat flux generation; however, in practical cases, microprocessor chips generate non-uniform heat flux [17]. The flexibility in microjet arrangement allows it to be directly embedded within the device to alleviate the hotspots and thermal stresses.

Many investigations were conducted over different designs to improve the heat transfer effects in the heat sink within the constraints such as maximum pumping power, maximum pressure drop, a constant flow rate of cooling fluid, constant heat flux, constant interface area, or constant cross-sectional area. Hajmohammadi et al. [7] investigated to obtain minimum peak temperature in the substrate, for optimal location and size of heat sink attachments, with the total area of heat sink remaining constant. It was suggested that dividing the large heat sink into a group of smaller ones enhances thermal performance by reducing peak temperatures. Husain et al. [8] numerically investigated a hybrid design of microchannel comprising of a microchannel pillar and jet impingement technique, which enhances heat transfer rate as the stagnation effect under jet diminishes due to channel flow. Zhang [9] combined a slot jet with microchannel design and analyzed the hybrid design numerically. For the same cross-sectional area, three channel shapes were identified, such as circular, trapezoidal, and rectangular. Results showed that the module with a circular channel had the highest pressure drop while the trapezoidal channel shape has the least substrate temperature. The highest pressure drop in the circular shape is attributed to the strong vorticities formed perpendicular to the direction of flow. Peng et al. [10] performed a comparative analysis between traditional microchannel heat sink (TMC) and multi-jet microchannel (MJMC) heat sink, using numerical techniques. They found that temperature uniformity improves significantly with the help of an increasing number of jets. The cooling performance of MJMC surpasses TMC. Husain et al. [11] compared the effect of different flow spent schemes with and without extraction ports. They presented that fluid removed from edges in unconfined flows have higher temperature uniformities and lower pressure drop. Qidwai et al. [4] found that temperature uniformity is obtained at the trade-off with Nusselt number in their study of a diverging channel. Han et al. [12] combined the microchannel with the microjet to capture the merits of both designs. The full-length trenches take away the spent fluid, which helps to reduce the crossflow effect. The maximum temperature distribution remains within the range of 25 °C. Deng et al. [14] investigated slot jet array and hybrid microchannel slot jet array heat sinks in the range of Reynolds number 230 to 2760. They affirmed that a hybrid heat sink has better performance for the same flow rate. Xie et al. [18] presented an argument that a uniform heat flux condition is acceptable for chips to surface only, rather than for heat sink surface, as assumed in many studies. The study included a combination of flow arrangements and different positions of heat-generating chips. The investigation highlighted that chip position significantly affects the temperature distribution. Lelea [19] found that the position and number of jets strongly influence hydrothermal performance. For fixed pumping power, the minimum temperature is higher for heat sinks with multiple inlets, and the temperature difference across the substrate is reduced. Wiriyasart and Naphon [20] considered parametric study of different fin structures liquid jet impingement to provide guidelines for heat sink design for minimum thermal resistance.

In most of the investigations, hybrid designs are analyzed for uniform heat flux conditions. For non-uniform heat flux, Yoon et al. [21], considered thermal resistance as an indicator to obtain the optimal position of partial heating for uniform and non-uniform heat flux conditions. They concluded that partial heating located at the centre behind the heat sink is the optimal position irrespective of the non-uniformity of heat flux. Sharma et al. [17] proposed a one-dimensional semi-empirical modeling approach for quick determination of initial design of microchannel heat sink targeting hot-spots. Design parameters in a microchannel heat sink were optimized by Hadad et al. [22] for non-uniform heat flux, using the RSM and JAYA algorithms together. The other work by Hadad et al. [23] investigated the shape optimization of a water-cooled impingement micro-channel heat sink, including manifolds.

1.1. Regression Based Surrogate Model Optimization Techniques

Parno et al. [24] implemented a Design and Analysis of Computer Experiments (DACE) surrogate model with PSO to optimize the solution of the groundwater problem. The sample of DACE was obtained using Latin Hypercube sampling to provide wide distributions of initial particles to yield robust convergence. Further, different variations of surrogate-assisted PSO excelled in different applications such that among Local PSO, adaptive PSO, and Global PSO; Local PSO outperformed the other two for problems with more than 10 dimensions. To improve heat transfer characteristics in a microchannel heat sink with the least trade-off of friction factor, Seo et al. [25] put forward another method wherein the objective functions were formulated by two surrogate models: Response Surface Approximation and the Kriging model. The Pareto Optimal Front obtained using a multi-objective Genetic Algorithm for two surrogate models were compared based on convergence and divergence criteria. The results of the two models were combined using K-means clustering to obtain the best Pareto Optimal front.

Many researchers proposed a design methodology intending to increase thermal efficiency with low-pressure drops. Shi et al. [26] considered Response Surface Approximation (RSA) combined with the NSGA-II algorithm to obtain Pareto solutions to optimize the selected design parameters. They found that the ratio of secondary channel width to microchannel width is the most influential parameter. Park et al. [27] evaluated design parameters in the wavy channel with grooves with the help of a multi-objective Genetic algorithm coupled with RSA surrogate model. Kulkarni et al. [28] proposed a design methodology for the heat sink to obtain lower thermal resistance. They used RSA with a multi-objective genetic algorithm to obtain Pareto-optimal design parameters.

1.2. Machine Learning-Based Surrogate Model Optimization Techniques

Machine learning algorithms evolve constantly with the aim for better prediction and the least amount of time and effort. They have been implemented in several fields such as intelligent communications [29], 3D printing technology [30], hydrological science [31], and many more. Many studies are conducted, in order to take advantage of machine learning, to improve heat transfer in heat sinks. Some investigations are focused on shape optimization, such Krzywanski [32], who introduced the design methodology of a bio-inspired falling-film heat exchanger using the combination of Genetic Algorithm (GA) and Artificial Neural Network (ANN). They reported that the AGENN model can successfully determine the required design parameters and operating conditions to generate the desired total heat transfer rate of the heat exchanger. Han et al. [33] used surrogate-based shape optimization of the wing-body of aircraft. They used Latin Hypercube sampling; however, to choose new sample points they used infill-criterion, so as to generate new designs based on the known designs. To explore the sensitivity of the change in design parameters, the investigation comprised of Maximum-Likelihood Estimation (MLE), with its only drawback that large sample points are needed to estimate the correlation. They suggested using surrogate-based models several times to reach the global optimum.

Xi et al. [34] used the MATLAB function nonlinear fitting model ‘nlinfit’ to obtain a correlation of 90 experimental data points to determine the Nusselt number. Then, the Back Propagation Neural Network (BPNN) was used to train the model with a Genetic Algorithm (GA) to improve heat transfer performance of ribbed channels. In the combined GA-BPNN algorithm, GA was used to optimize the weights and threshold of BPNN so that, together, they could reach the global optimum solution. Singh [35] addressed the problem of turbine-blade cooling using a slot jet cooling mechanism. The location of slot jet impingement on a concave surface was optimized using hybrid feed-forward Artificial Neural Network (ANN) and Genetic Algorithm (GA). ANN was used to train nonlinear data on 175 cases in the design space. The predicted optimal solution of ANN-GA was simulated in CFD simulation, for which the difference in output was negligible. Shi et al. [36] proposed a methodology for modifying heat exchanger design using CFD coupled with the Radial Basis Neural Network surrogate model and genetic algorithm to achieve maximum flow uniformity.

1.3. Optimization Techniques

GA algorithms can be sub-divided into three operations: Mutation, Crossover, and Selection. Tan et al. [37] suggested that parameters in GA are chosen based on the rule of thumb. There is no standard to select a proper strategy or parameter or their combination, which may lead to poor convergence and trapping in local optima. Additionally, these algorithms need to be run a thousand times before reaching global optima. Samii [38] compared the two different approaches and concluded that PSO is simpler to use based on the fact that it has higher convergence compared with GA, since, when the individual possesses the same genetic material, the crossover effect is almost eliminated. For structural optimization Eagle strategy, the PSO algorithm combines the advantages of global search and intense local search [39]. The Eagle strategy uses the Lévy walk method to identify the search domain, then PSO is applied in the second stage. This efficient combination reduces the computation time and increases the probability of achieving global optima.

Bello-Ochende et al. [40] used the Constructal theory in combination with CFD to reduce maximum wall temperature; however, the originator of the Constructal theory considers it to be a universally valid theory and not an optimality statement [41].

Entropy Generation minimization (EGM) coupled with CFD is another methodology, based on the first principle, applied in many investigations [42,43].

Some studies preferred to use the simplified conjugate gradient method [44,45], which uses sensitivity analysis to establish the conjugate direction. However, in this method, the probability of trapping in the local optima is quite high [46].

1.4. Hybrid Design Schemes

Some studies with non-uniform heat flux focussed on minimizing the maximum substrate temperature. Cho et al. [47] performed an experimental and numerical analysis of the different combinations of channel geometries and header designs. Three types of thermal loads were chosen to obtain peak temperature and temperature distribution in the substrate. Their results showed that a reasonable trade-off between peak temperature and pressure drop is a must, as pressure drop is higher in the diverging channel with a trapezoidal header design with the lowest peak temperature and better temperature distribution. Abo-Zahhad et al. [48] implemented hybrid schemes with microjet and microchannel for thermal management in high concentrator photovoltaic solar cells. They used exergy efficiency (electrical and thermal) for the performance comparison among different layouts of cell structure. The lowest temperature non-uniformity was found to be 5.8 °C. The main objective of their study was to avoid exergy loss, i.e., to have uniform cooling at the lowest possible flow rates. Muszynski and Mikielewicz [49] investigated the arrangement of an array of jets on heat transfer in the heat sink. The fluid interaction of fully submerged jets is very complex; therefore, they preferred an equal-nozzle spacing design for the analysis.

1.5. Research Gap

The non-uniformity in temperature distribution arises due to the non-uniform heat flux generated in the processors. However, due to the miniature design and high thermal conductivity of the substrate material, it is convenient to assume uniform heat flux boundary in most studies, but it is not an absolute true case. According to the literature review, most of the design strategies assumed uniform heat flux boundary conditions in the numerical analysis. However, in experimental conditions, the miniaturized heat sinks have to deal with non-uniform heat flux conditions. The microjet heat sink geometry dealing with hot-spots assumes equal jet spacing to eliminate the complexity of multi-jet interaction in the fluid domain. To the best of the authors’ knowledge, no heat sink designed with non-uniform jet spacing is investigated to date. The objective of this study was to minimize thermal resistance while keeping a check on temperature non-uniformity of the heat sink by providing a strategy for practical design of a water-cooled microjet heat sink, which utilizes the significance of machine learning techniques to empower the design process to yield a better design to tackle ever-increasing challenges.

2. Computational Model

2.1. Flow Configuration and Numerical Methods

Figure 1 shows the configuration of the complete computational domain. The dimensions of the heat sink are 4000 × 3000 × 250 μm (L × W × H), with a copper substrate thickness of 50 μm. The inlet of fluid takes place from the holes provided at the top plate with a fixed diameter of 400 μm, and the spent fluid carries away from one side only. Non-uniform heat flux was applied at the bottom of the heat sink as per the power map adopted from the study of Sharma et al. [17]. The design has inherent limitations of fixed geometrical location hot spots and fixed region of inlet holes. Operational limitations include fixed coolant inlet temperature, constant mass flow rate (0.00025 kg/s), and the power dissipated at the substrate base. The thermo-physical properties of copper and water are given in

Table 1. Thermophysical properties of liquid and solid state at 300 K.

	${\mathit{c}}_{\mathit{p}}(\mathbf{J}/\mathbf{kg}·\mathbf{K})$	$\mathit{k}$ (W/m·K)	$\mathit{\rho }$ (kg/m3)	$\mathit{\mu }$ (kg/m·s)
Water	4181	0.614	997	$2.414\times {10}^{-5}\times {10}^{247.8/(T-140)}$ [50]
Copper	381	388	8978

.

The following assumptions were made for the computational domain.

• Continuum regime applies, in which Navier-Stokes equation remains valid.
• The flow is three-dimensional, incompressible, and steady-state for energy and fluid flow.
• Thermo-physical solid and fluid properties are constant, except fluid viscosity.
• Buoyancy and viscous effects are negligible.

The above assumptions guide the governing Equations (1)–(4) for both cooling fluid and solid substrate domain as follows:

For coolant (liquid phase).

Conservation of mass for coolant (liquid phase).

$$\nabla ·\overrightarrow{\mathrm{V}}=0$$

(1)

Momentum for coolant (liquid phase).

$${\rho }_f(\overrightarrow{\mathrm{V}}·\nabla )\overrightarrow{\mathrm{V}}=-\nabla \mathrm{p}+{\mu }_f{\nabla }^2\overrightarrow{\mathrm{V}}$$

(2)

Energy for coolant (liquid phase).

$${\rho }_f{c}_{p,f}(\overrightarrow{V}·\nabla )T_f=k_f{\nabla }^2{T}_f$$

(3)

Energy for substrate (solid phase).

$${\rho }_f{c}_{p,f}(\overrightarrow{V}·\nabla )T_f=k_f{\nabla }^2{T}_f$$

(4)

The Ansys Mesh modeler was used for discretization of the computational domain. The conservation equations of mass, momentum, and energy equations were solved using Ansys CFX. The transport equations were solved with the laminar model. The Rhie-Chow algorithm was used to couple pressure and velocity fields. The convergence criterion was set to be ${10}^{-6}$ for all the equations. Since there is no symmetry in the geometry, complete domain analysis is necessary. Non-uniform heat flux was applied at the bottom wall of the heat sink, as shown in Figure 1b. All the other sides of the heat sink domain were considered insulated, and heat is convected away from the heat sink by fluid only. The mass flow inlet boundary condition was applied with a uniform temperature field across the inlet cross-sections, and the pressure outlet condition is applied at the outlet cross-section. At the interface between the heat sink and the coolant, a no-slip condition was chosen with heat transfer coupling between the two domains. The boundary conditions are summarized in

Table 2. Boundary condition of numerical method.

Boundary Condition	Momentum	Thermal
Fluid domain	Mass flow inlet/Pressure outlet	$T_{in}=290 K$
Heat sink domain at y = 0	$\overrightarrow{V}=0$	${-k_s\frac{\partial T}{\partial \eta }|}_b=q_k$
Heat sink domain at x = 0, x = W, y = H, z = 0, z = L	$\overrightarrow{V}=0$	$\frac{\partial T}{\partial \eta }=0$ (adiabatic)
Interface between fluid and substrate	$\overrightarrow{V}=0$ (No slip condition)	$T_{f,\mathsf{\Gamma }}=T_{s,\mathsf{\Gamma }}$ $-k_f{\frac{\partial T_f}{\partial \eta }|}_{\mathsf{\Gamma }}={-k_s\frac{\partial T}{\partial \eta }|}_{\mathsf{\Gamma }}$

.

The Reynolds number was found to be within the laminar regime at both inlet and outlet cross-section for the given mass flow rate. Based on the respective hydraulic diameter given in Equation (5), the Reynolds number was calculated as given in Equations (6) and (7).

$$d_h=\frac{4A_c}{\varrho }$$

(5)

$$R{e}_{in}=\frac{{\dot{m}}_{in}{d}_{h,in}}{A_{c,in}{\mu }_f}$$

(6)

$$R{e}_o=\frac{{\dot{m}}_o{d}_{h,o}}{A_{c,o}{\mu }_f}$$

(7)

2.2. Grid Independence Test

In order to save computing time and memory, the size of spatial gird size was determined by the trial-and-error method. The gird independence test was considered for thermal resistance and mean substrate temperature, as shown in Figure 2. The difference in results with grid systems 6 million, 9.1 million, and 10.44 million elements is much less than the maximum allowed limit. Thus, a grid with 9.1 million elements was selected for further calculations.

2.3. Validation

The current study’s numerical values were compared with the experimental results for a single-phase microjet silicon heat sink (k = 148 W/m K) as reported by Wang et al. [51]. The test cases with the number of nozzles 13 and 4 with a nozzle diameter of 40 μm and 50 μm were subjected to flow rates of 2 mL/min and 8 mL/min, for a power range of 3.8 to 23.8 W, respectively. A direct comparison between experimental and numerical results for maximum temperature rise in the substrate as a function of power is shown in Figure 3, with a deviation within 15%. This deviation is attributed to several factors such as simplification of the physical model and experimental uncertainties.

2.4. Computational Complexity and Implementation Cost

The initial cost of the CFD setup is very high, which includes a high-end computing machine and annual license fees along with the operating cost associated with manpower. However, this setup can be used to solve many problems associated with thermal and fluid interactions, which reduces its overall expenditure per problem.

In the CFD simulation, the Navier-Stokes equation was discretized over a mesh, formulated separately for each case. This mesh may vary in the number of elements, based on which mesh can be categorized between “coarse” to “very fine”. The “coarse” mesh takes less time compared with “very fine” mesh for simulating the same case. However, the results are far from reality for coarse mesh, whereas results of “very fine” mesh approach values of exact solutions, at the cost of consumption of a large amount of computation power and time. Therefore, the number of elements in a mesh needs to be selected appropriately, and the same was carried out by performing a Grid Independence test.

The system used for CFD simulation comprised an Intel Xeon Processor E5-1620 v4 (4 Core, 3.5 GHz, 3.8 GHz Turbo, 2400 MHz, 10 MB, 140 W) and RAM 32 GB (8 × 4 GB) 2400 MHz DDR4. On this machine, each case took more than 10 h of simulation, as the time taken for each case depends on the complexity of each case, the number of elements, and the combination of physical models used. Further, the result at discrete values of design variables can only be obtained one at a time. In this scenario, it is not feasible to search optimum values in continuous space with the help of the CFD technique alone. To tackle this challenge, different techniques were used, as discussed in the literature.

To save a large amount of computation cost and time, previous researchers have used CFD in combination with AI models, which reduces time to reach optimum values within a few minutes. This is the huge advantage of using combinatorial algorithms.

3. Optimization Method

3.1. Design Space, Sample Points, and Objective Function

The flow configuration remains fixed in earlier studies and assumes uniform base heat flux with modifications in geometry for parametric investigation of design variables without considering the local heat-transfer effects. For each new parameter of the design variable, a complicated numerical analysis must be performed to obtain its flow and heat transfer characteristics. The top plate was provided with 9 inlet holes in this numerical model, as shown in Figure 4. Each hole remains completely in the region marked by section lines A, B, A’, and B’. From the lower-left corner, the A and B locations are 0.96 mm and 1.93 mm, while for A’ and B’, they were 1.3 mm and 2.6 mm, respectively. The center of each hole was placed such that at least one jet was available for cooling in the marked section. The design space was set as:

$0.2<x_{11},x_{21},x_{31}<0.725$;	$0.2<y_{11}, y_{12}, y_{13}<0.975$;
$1.125<x_{12}, x_{22}, x_{32}<1.1775$;	$1.375<y_{21}, y_{22},y_{23}<2.525$;
$2.175<x_{13}, x_{23},x_{33}<2.7$;	$2.925<y_{31},y_{32},y_{33}<3.7$;

The overall thermal resistance and temperature uniformity, which is defined as the standard deviation from the mean substrate temperature [11], are given as in Equations (8) and (9).

$$f_1=\frac{T_{b,max}-T_{f,in}}{q{A}_b}$$

(8)

$$f_2=\sqrt{{{\displaystyle \sum }}_{k=1}^{\tilde{N}}\frac{{(T_k-T_{mean})}^2}{\tilde{N}}}$$

(9)

where, ‘k’ represents the different heat flux regions ($q_{1,}{q}_{2,}{q}_{3,}{q}_{4,}{q}_5)$.

When considering two objective functions to be minimized, the weighted objective function may be defined. However, there is no basis for defining any weightage in the current case. Therefore, in this work, $f_1$ will be the prime objective function and $f_2$ will be checked simultaneously after optimization.

With the help of the Latin hypercube sampling technique, different positions of the hole were obtained, which are given in

Table 3. Sample of centre of inlet holes on the top plate of the substrate.

Sample	11		12		13		21		22		23		31		32		33
	x	y	x	y	x	y	x	y	x	y	x	y	x	y	x	y	x	y
1	0.46	0.73	1.74	0.53	2.53	0.49	0.60	1.55	1.21	2.52	2.50	2.16	0.55	3.50	1.58	3.62	2.36	3.46
2	0.64	0.24	1.76	0.93	2.35	0.62	0.54	2.21	1.19	1.68	2.46	1.42	0.72	3.05	1.71	3.03	2.42	3.23
3	0.66	0.39	1.25	0.28	2.62	0.87	0.41	1.62	1.44	1.78	2.66	2.05	0.30	3.24	1.42	3.19	2.18	3.63
4	0.26	0.92	1.32	0.73	2.58	0.76	0.70	1.83	1.48	2.10	2.39	1.66	0.34	3.34	1.39	3.69	2.27	2.93
5	0.51	0.89	1.63	0.86	2.38	0.84	0.48	2.04	1.76	1.46	2.18	2.52	0.67	3.37	1.35	3.26	2.29	3.16
6	0.71	0.28	1.14	0.89	2.35	0.32	0.50	2.51	1.72	2.26	2.20	1.52	0.24	3.57	1.15	3.65	2.21	3.56
7	0.37	0.30	1.36	0.62	2.57	0.74	0.25	1.93	1.54	1.64	2.62	1.98	0.26	3.56	1.33	3.29	2.52	3.01
8	0.24	0.83	1.65	0.60	2.40	0.54	0.65	1.69	1.64	1.99	2.70	1.91	0.54	3.62	1.70	3.65	2.59	3.58
9	0.21	0.49	1.31	0.77	2.70	0.48	0.42	1.88	1.30	1.73	2.56	1.87	0.22	3.31	1.25	3.33	2.34	3.21
10	0.35	0.85	1.68	0.40	2.48	0.21	0.72	1.84	1.18	1.98	2.38	2.35	0.28	3.21	1.51	3.00	2.58	3.37
11	0.55	0.56	1.21	0.22	2.26	0.66	0.62	2.34	1.62	2.05	2.55	1.60	0.45	3.67	1.17	3.10	2.69	3.26
12	0.32	0.68	1.56	0.37	2.60	0.78	0.30	1.43	1.50	1.55	2.29	1.93	0.51	3.01	1.67	3.44	2.37	3.07
13	0.59	0.42	1.28	0.50	2.51	0.27	0.28	2.47	1.73	1.83	2.65	2.46	0.37	3.43	1.77	3.24	2.67	3.30
14	0.49	0.48	1.60	0.31	2.28	0.41	0.20	2.39	1.60	1.49	2.61	1.99	0.44	3.14	1.38	3.56	2.39	3.50
15	0.45	0.61	1.36	0.92	2.33	0.92	0.35	2.44	1.71	2.28	2.32	1.38	0.67	3.69	1.32	2.96	2.23	2.97
16	0.72	0.78	1.29	0.26	2.66	0.37	0.32	2.01	1.16	1.88	2.44	1.76	0.42	3.48	1.75	3.13	2.22	3.47
17	0.27	0.80	1.54	0.35	2.31	0.94	0.33	2.10	1.53	2.15	2.68	2.32	0.52	3.32	1.66	3.38	2.56	3.11
18	0.38	0.33	1.62	0.44	2.22	0.53	0.43	1.47	1.23	2.42	2.34	1.75	0.32	3.60	1.61	3.49	2.63	3.35
19	0.52	0.62	1.42	0.68	2.43	0.85	0.39	1.59	1.41	2.12	2.36	2.11	0.40	3.53	1.46	3.46	2.61	3.40
20	0.58	0.46	1.69	0.83	2.18	0.45	0.63	2.25	1.35	1.84	2.23	2.07	0.63	3.46	1.54	3.18	2.32	3.60
21	0.61	0.21	1.47	0.66	2.45	0.68	0.48	2.32	1.40	2.21	2.31	1.63	0.69	3.04	1.58	3.54	2.65	3.14
22	0.22	0.97	1.15	0.82	2.41	0.25	0.68	1.52	1.68	1.60	2.27	2.44	0.49	2.93	1.45	3.02	2.41	3.33
23	0.63	0.71	1.51	0.56	2.64	0.34	0.57	2.16	1.32	2.40	2.53	1.56	0.64	2.99	1.64	3.52	2.53	3.65
24	0.43	0.56	1.44	0.46	2.47	0.81	0.56	1.75	1.13	2.48	2.23	2.26	0.59	3.28	1.26	3.08	2.44	3.04
25	0.55	0.95	1.49	0.77	2.23	0.57	0.66	1.66	1.57	2.32	2.48	1.80	0.31	3.08	1.20	3.23	2.30	2.98
26	0.29	0.76	1.18	0.41	2.50	0.61	0.26	2.28	1.31	1.38	2.58	2.19	0.21	2.95	1.48	2.93	2.48	3.52
27	0.43	0.65	1.71	0.24	2.30	0.72	0.23	1.77	1.36	1.93	2.25	1.68	0.60	3.17	1.52	3.34	2.51	3.27
28	0.32	0.37	1.21	0.97	2.67	0.96	0.46	1.96	1.26	1.70	2.52	2.40	0.47	3.13	1.28	3.11	2.25	3.69
29	0.40	0.53	1.39	0.71	2.20	0.28	0.37	1.41	1.46	1.43	2.44	2.28	0.39	3.19	1.23	3.40	2.47	3.09
30	0.69	0.26	1.53	0.57	2.55	0.39	0.52	2.13	1.66	2.35	2.42	1.47	0.58	3.41	1.16	3.59	2.65	3.43

. Latin hypercube sampling (LHS) is a modified Monte Carlo sampling method that divides sample design variables with nonoverlapping intervals [52]. It provides a random sampling combination of each design variable in the design space. In the present case, 30 sample points were obtained from sampling to satisfactorily train the surrogate model using CFD results [36].

3.2. Surrogate Model

RBFN is successfully implemented to train models [53]. The Radial Basis Neural Network (RBNN) function is established using MATLAB 2015, with newrb function, whose maximum number of neurons needs to be provided as an input. It is a three-layer network that includes the input layer, hidden layer, and output layer, as shown in Figure 5. It has a strong ability to predict the objective function with a rapid convergence rate. The distance between selected points in the design space and the sample data points forms the basis to obtain the influence of the data points. This distance is used for training of AI model on the sample data points to calculate the inertia weight (${\mathrm{w}}_{\mathrm{m}}^{\prime })$. The $\mathsf{\gamma }$ constant determines the radial influence of the Gaussian function, which is of less significance in the case of a large sample in the design space. However, for small data sets, the value of $\mathsf{\gamma }$ needs careful consideration before any value is assigned to it. The RBNN function can be summarized in the following Equation (10) [54].

$$h(x,y)={{\displaystyle \sum }}_{m=1}^{Ɲ}\underset{Basis function}{\underbrace{w_m^{\prime } \exp (\underset{Radial}{\underbrace{-\gamma {\Vert X-\upsilon \Vert }^2}})}}$$

(10)

3.3. Particle Swarm Optimization Algorithm (PSO)

The Particle Swarm Optimization (PSO), first proposed by Kennedy and Eberhart [55], is a metaheuristic algorithm inspired by the combined intelligence of a swarm of birds in search of food and shelter. This evolutionary algorithm initializes the optimization process with a population of random solutions, which is then passed on to each particle, based on which each particle updates its velocity and position. As soon as the new velocities and positions are updated, the new optimum solution appears in the design space. Again, based on the new optimum solution, the velocities and position of each particle are updated. This process continues until the convergence criterion is reached by using Equations (11) and (12):

$$V_{ij}^{t+1}=w^{\prime }{V}_{ij}^t+C_1{r}_1^t(Pbes{t}_{ij}-X_{ij}^t)+c_2{r}_2^t{(gbes{t}_j-X)}_{ij}^t$$

(11)

$$X_{ij}^{t+1}=X_{ij}^t+V_{ij}^{t+1}$$

(12)

The procedure of PSO is summarized in Figure 6. The initial solution is represented by the RBNN function’s parameters, including input weights, output weights, and values of design parameters. For each iteration, the position and velocities of the particles were updated and for which the optimal solution was obtained, which is given by Equation (8).

The overall methodology can be summarized in the following steps:

• Identify the design and response parameters.
• Determine the constraint values of design parameters.
• Identify the fixed parameters.
• Obtain values of design parameters using a suitable sampling method.
• Calculate the values of response parameters at all design points using CFD techniques.
• Use design and response parameters to train the Radial Basis Neural Network model.
• Search for optimal design parameters using the Particle Swarm Optimization algorithm.
• The optimal design points are used in CFD analysis for verifying the result of PSO algorithm.

3.4. Parameter and Design for RBNN and PSO

In most of the studies, a Bayesian approach is preferred for comparing the performance of algorithms [56] and to check that the algorithm converges to a global minimum solution. To ensure that optimization does not stop far away from the global minimum, another similar statistical technique, Analysis of variance (ANOVA), was used. The coupled algorithm has variable inputs that include Gaussian Spread ($\gamma$), Number of particles (P), and Number of iterations (ɨ) at user discretion. Each variable at three levels was used in this study to understand the influence on the optimized results.

Table 4. Input parameters and their levels (L).

Factor/Parameter	Symbol	L-1	L-2	L-3
Gaussian Spread	γ	0.3	0.6	0.9
No. of Particles	P	50	125	200
No. of Iterations	ɨ	100	250	500

presents the process parameters with their levels. After performing simulations, thermal resistance and temperature uniformity were measured, which served the purpose of response variables. Nine experimental runs with three replicates were performed as per the L9 (33) orthogonal array (OA) and presented in

Table 5. Simulation plan according to L₉ (3³) Orthogonal array.

Run	γ	P	ɨ
1	1	1	1
2	1	2	2
3	1	3	3
4	2	1	3
5	2	2	1
6	2	3	2
7	3	1	2
8	3	2	3
9	3	3	1

.

4. Results and Discussion

The CFD approach simulates the heat transfer and flow characteristics in the heat sink with a cooling jet with laminar flow. These results were validated with the experimental results [51]. The CFD results were used to train the surrogate model RBNN to form a continuous design space, on which the particles of PSO algorithm were allowed to search for optimum design parameters.

Table 6. Predicted optimized position of inlet holes by the PSO algorithm.

11		12		13		21		22		23		31		32		33
x	y	x	y	x	y	x	y	x	y	x	y	x	y	x	y	x	y
0.20	0.84	1.31	0.59	2.70	0.28	0.62	1.65	1.13	1.88	2.45	2.17	0.20	3.12	1.28	2.93	2.39	3.37

shows the predicted values of an optimal solution obtained by the algorithm. Using these optimal solutions, the CFD case was prepared and simulated, for which the temperature contour was plotted (Figure 7) and in which the maximum temperature was found to be 343.38 K. As shown in

Table 7. Comparison of values predicted by the AI model and CFD results.

	f1	f2
Predicted Value	6.205	8.019
Actual Value	6.170	7.914
%e	0.565	1.323

, for the predicted values and CFD-calculated values of thermal resistance and temperature uniformity, the percentage error was less than 2%. The metaheuristic algorithms are a powerful tool that provide the optimal solution rapidly. However, careful validation with the actual results is necessary for these algorithms since other variables influence the optimal solution, which is obtained by the trial-and-error method. Therefore, multiple solutions must be obtained prior to validation.

4.1. Effect of PSO Parameters

In Figure 8, the value of the objective function thermal resistance is found to be reducing with each iteration of the PSO algorithm. The searching for optimal parameters in the design space with infinite points would be computationally expensive if a metaheuristic algorithm would not have been applied. However, even with quick search and fast convergence features, these algorithms have some issues such as unknown variables that significantly influence the solution. Additionally, there is a chance that solutions can be trapped in local optima. In order to avoid trapping in local optima, it is necessary to run multiple simulations. ANOVA was performed on the measured values of response variables, as presented in

Table 8. Measurement of the response variable after performing simulations as per Taguchi L9 OA.

Run No.				Exp1	Exp2	Exp3	Average
	γ	P	ɨ	f1	f1	f1	f1
1	0.3	50	100	6.357	6.352	6.340	6.350
2	0.3	125	250	6.312	6.522	6.522	6.452
3	0.3	200	500	6.464	6.501	6.522	6.496
4	0.6	50	500	6.226	6.226	6.226	6.226
5	0.6	125	100	6.232	6.229	6.227	6.229
6	0.6	200	250	6.226	6.226	6.226	6.226
7	0.9	50	250	6.225	6.222	6.215	6.221
8	0.9	125	500	6.211	6.206	6.211	6.209
9	0.9	200	100	6.222	6.207	6.212	6.214

. This was carried out to understand the effect of unknown variable inputs in the search algorithm and to obtain the most dominant parameter that significantly affects the thermal resistance (f₁).

4.2. Analysis of Variance (ANOVA)

ANOVA (

Table 9. ANOVA Results.

Source of Variation	Degrees of Freedom	Sum of Squares	Mean Square	F-Ratio	p-Value	Percentage Contribution
Gaussian spread	2	0.169181	0.084591	20.53	0.046	89.10
No. of population	2	0.006097	0.003049	0.74	0.575	3.21
No. of iterations	2	0.006358	0.003179	0.77	0.564	3.35
Error	2	0.008241	0.004121			4.34
Total	8	0.189877				100

) was used to determine which parameter significantly influences the thermal resistance.

ANOVA results revealed that gaussian spread is the most significant parameter to influence the thermal resistance, followed by the number of iterations and number of populations. Gaussian spread affects the values of the output of nearby values of sample data points, which happened to be the center of the normal distribution curve from the center, as determined by the RBNN function, and this distance is multiplied by the gaussian spread. Since the gaussian spread is at user discretion, careful consideration is a necessity for obtaining adequate heat transfer to obtain good thermal performance.

5. Conclusions

In this work, a 3-D numerical study of microjets’ position was conducted, based on the samples provided by Latin Hypercube sampling. The results obtained from physical simulation were used to train the surrogate model RBNN coupled with the PSO algorithm to obtain the optimized location of microjets which minimize substrate temperature for non-uniform heat flux distribution. The predicted values by surrogate models were compared with the results of CFD simulation, for which the percentage error was 0.53% for the target function thermal resistance (f₁), and for temperature uniformity (f₂) it was 1.3%. The advantage of this method is that earlier studies considered equal jet spacing, since fluid interaction in the multi-jet domain is complicated and difficult to quantify. This method accommodates complications of fluid interaction within the design space. The Particle Swarm Optimization is a metaheuristic algorithm which is used to search the design space for the optimal solutions. The influence of three parameters in the coupled surrogate model was analyzed. Results from the coupled algorithm vary slightly for each complete run. Therefore, three experiments were performed for each combination of selected parameters. ANOVA was performed on the measured values of thermal resistance, and it was found that gaussian spread is the most dominant parameter affecting the thermal resistance, followed by the number of iterations and number of populations for the given range of process parameters.