# NSGA-II-Based Microchannel Structure Optimization Problem Study

Authors: Linlin Wang
Journal: Scientific Programming
 (2022)

---

## Abstract

A multiobjective genetic algorithm is used to optimize the structure of circular concave cavities and microchannels. The objective functions of thermal resistance and pumping power are constructed by the response plane approximation method according to the simulation results, and then a mathematical model of multiobjective genetic optimization with the structural parameters of microchannels as variables is established. The Pareto optimized solution sets of thermal resistance and pumping power are calculated by the nondominated ranking genetic algorithm NSGA-II, and the comprehensive heat transfer performance is evaluated by the enhanced heat transfer factor. The results show that the multivariate statistical coefficientsR2 of the thermal resistance and pumping power objective functions are 0.932 9 and 0.996 6, respectively, indicating the high accuracy of the fitted functions. The optimized channel structure (e1=0.036.8mm, e2=0.019.3mm) was used to achieve a more uniform temperature field distribution and better integrated heat transfer performance (enhanced heat transfer factor η=1.23). When the thermal resistance is larger or the pump work is larger, the comprehensive heat transfer effect is not as good as the working condition when the thermal resistance and pump work are more uniform.


---

## Body

## 1. Introduction

With the development of microfabrication technology, one of the effective methods of using microchannel high power intensive microelectronic devices with high surface-to-body ratio is widely used in advanced engineering fields such as microelectronics, energy, military nuclear energy, and biochemical industry [1–4]. According to microchannel convective heat transfer theory, the heat transfer performance of microchannel heat sink can be improved by increasing the heat transfer surface area of the channel or improving the heat transfer performance of the fluid, or use boiling heat transfer, such as nanofluids [5–9].It has been shown [10–12] that the heat transfer performance of complex structured microchannel heat sinks is excellent because the continuous variation of the channel cross section can interrupt the thermal boundary layer development and keep the flow and heat transfer in an underdeveloped state, which serves to enhance the heat transfer; moreover, the setting of rough elements inside the channel enhances the internal disturbance, which in turn has an important impact on the flow and heat transfer characteristics [13, 14].In [15], the structures of microchannel heat sink with fan-shaped cavities and microchannel heat sink with triangular-shaped cavities were simulated and optimized, and it was concluded that the cavity action can destroy the microchannel flow and heat transfer boundary layer, thus enhancing the microchannel heat transfer. In [16], the effects of rib spacing ratio to pipe diameter (p/d), rib height ratio to pipe diameter (e/d), and Reynolds number on the flow and heat transfer in circular tubes with microchannels were experimentally studied, and the Nusselt (Nu) number and friction coefficient correlation equations with the three parameters of p/d, e/d, and Reynolds number as variables were proposed to guide engineering applications. Yang and Cao [17] analyzed in detail the effect of sinusoidal structured rough elements on the flow and heat transfer in microchannels and analyzed that their enhanced heat transfer is due to the change in flow direction and the increase in area. In [18], a complex structured microchannel heat sink with a combination of fan-shaped cavities and inner ribs was proposed, and its single-phase flow and heat transfer were investigated by experimental and simulation methods, and the results showed that its integrated heat transfer performance was significantly better than that of an equal-section microchannel at Reynolds numbers greater than 300. Therefore, complex structured microchannels have obvious advantages in heat transfer.A reasonable channel structure design can also improve the heat transfer effect. Therefore, heat transfer efficiency and pressure drop are two important parameters to measure the comprehensive heat transfer performance of microchannels. The engineering design generally involves multiparameter optimization, called multiobjective parameter optimization [19]. Hemmat Esfe et al. [13] used a nondominated ranking genetic algorithm (NSGA-II) with Nusselt number and friction coefficient as optimization objectives for structural optimization of ribbed and cavity channels with rib height, cavity diameter, height and location as variables.In [15], the flow and heat transfer in variable cross section and microchannel heat sink with rough elements were investigated, and two objective functions of thermal resistance and pumping power were optimized using a multiobjective genetic algorithm to obtain the Pareto optimization solution set. In [6], the multiobjective genetic algorithm was used to optimize the heat transfer coefficient and pressure drop for microchannels with concave structure using channel height, cavity spacing, and height as variables. Theoretically, the genetic algorithm is able to find the optimal solution of the objective function probabilistically and randomly, which is more adaptable than the conventional optimization algorithm.To overcome the limitations of straight tube microchannels, researchers have conducted extensive studies for split-recombination (SAR) microchannels. The upstream flow in a SAR microchannel is split into two equal streams, which then hit each other within the mixing element. A series of microchannel units with the same structure interconnect and further enhance the material process. Benefiting from the combination of multiple microchannel units, the synthesis rate in SAR microchannels is significantly greater than that of straight microchannels with only a single microchannel unit [11, 12].The most common problem with current microchannel structures is the inhomogeneous mixing of reactants caused by insufficient turbulence. The key to promote the mixing effect within microchannels is to further enhance the split-recombination effect of microchannels. At the same time, the SAR microchannel structure still has great potential to further optimize and enhance the passive mixing process. However, insufficient grid accuracy and inaccurate phase interface capture techniques result in the current numerical simulation methods in microchannels to remain imperfect. Therefore, based on the CLSVOF multiphase flow model and the SST k-ω turbulence model, this study investigates a series of microchannels with complex structures.Therefore, based on the CLSVOF multiphase flow model and SST k-ω turbulence model, this study simulates a series of microreactors with complex structures by using the CLSVOF multiphase flow model for two-dimensional numerical simulation of gas-liquid two-phase flow and carries out multivariate optimization of the throat size and baffle length by the proxy modeling method to predict and analyze the optimal structure of the microreactors.


## 2. Numerical Simulation Methods

### 2.1. CLSVOF Multiphase Flow Model

For two incompatible fluids, it is very difficult to capture the free phase interface during dynamic motion [13]. Therefore, Sussman and Puckett proposed the CLSVOF method [14]. The principle of this method is to use the convection equation of the VOF function to achieve mass conservation and the level-set function to calculate the curvature and normal vectors to capture the phase interface smoothly. In the CLSVOF method, the phase interface is reconstructed by a segmented linear scheme [15, 16], and the initialization of the distance function of the level-set equation is achieved based on the reconstructed phase interface. For the level-set model, the mass of the transport equation is not conserved and the normal vector of the curvature calculation in the VOF model is not accurate enough [15, 17, 18]. The CLSVOF model remedies the computational defects of the mass nonconservation and the inaccuracy of the curvature calculation in the level-set and VOF models and combines its computational advantages to capture the phase interface more accurately. Hadad et al. [19] simulated the underwater bomb shock wave by CLSVOF multiphase flow model. The CLSVOF model is able to capture the phase interface more accurately under strong interaction conditions.In the level-set model, the phase interface is defined as an equivalence surface belonging to an implicit functional equationϕx,t, and the value of ϕx,t at any moment can be determined by solving the controlling equation to obtain the zero equivalence surface, i.e., the location of the phase interface. ϕx,t function is defined as follows:(1)ϕ=dx,Γt<0,0,dx,Γt>0.The position of the phase interface at timet is represented by t, and dx,Γt is the distance function from the center of the grid to the phase interface Γt. In the VOF model, different fluids are incompatible with each other. The volume function of fluid within each cell in the VOF model is defined as follows:(2)c=fluid volume in cellunit volume.IfC = 1, the fluid of the grid is in continuous phase, C = 0 is in discrete phase, and 0 < C < 1 means that the grid is at the phase interface. Therefore, the position and shape of the phase interface in space can be determined by solving for the volume function C.In an incompressible fluid, there are(3)DCDt=∂C∂t+V⟶⋅∇⋅C=0,ρ∂V⟶∂t+ρV⟶⋅∇V⟶=−∇P+∇⋅V⟶∇V⟶+∇V⟶T+ρg+Fsv.The above equation is the convective transport equation for the volume functionc as well as the momentum conservation equation. Fsv is the surface tension term of the momentum equation.The CLSVOF model compensates the deficiencies of mass conservation and phase interface parameters of the level-set model and VOF model and combines the advantages of the two models. The CLSVOF model can combine the advantages of both models to capture the phase interface more accurately.


## 2.1. CLSVOF Multiphase Flow Model

For two incompatible fluids, it is very difficult to capture the free phase interface during dynamic motion [13]. Therefore, Sussman and Puckett proposed the CLSVOF method [14]. The principle of this method is to use the convection equation of the VOF function to achieve mass conservation and the level-set function to calculate the curvature and normal vectors to capture the phase interface smoothly. In the CLSVOF method, the phase interface is reconstructed by a segmented linear scheme [15, 16], and the initialization of the distance function of the level-set equation is achieved based on the reconstructed phase interface. For the level-set model, the mass of the transport equation is not conserved and the normal vector of the curvature calculation in the VOF model is not accurate enough [15, 17, 18]. The CLSVOF model remedies the computational defects of the mass nonconservation and the inaccuracy of the curvature calculation in the level-set and VOF models and combines its computational advantages to capture the phase interface more accurately. Hadad et al. [19] simulated the underwater bomb shock wave by CLSVOF multiphase flow model. The CLSVOF model is able to capture the phase interface more accurately under strong interaction conditions.In the level-set model, the phase interface is defined as an equivalence surface belonging to an implicit functional equationϕx,t, and the value of ϕx,t at any moment can be determined by solving the controlling equation to obtain the zero equivalence surface, i.e., the location of the phase interface. ϕx,t function is defined as follows:(1)ϕ=dx,Γt<0,0,dx,Γt>0.The position of the phase interface at timet is represented by t, and dx,Γt is the distance function from the center of the grid to the phase interface Γt. In the VOF model, different fluids are incompatible with each other. The volume function of fluid within each cell in the VOF model is defined as follows:(2)c=fluid volume in cellunit volume.IfC = 1, the fluid of the grid is in continuous phase, C = 0 is in discrete phase, and 0 < C < 1 means that the grid is at the phase interface. Therefore, the position and shape of the phase interface in space can be determined by solving for the volume function C.In an incompressible fluid, there are(3)DCDt=∂C∂t+V⟶⋅∇⋅C=0,ρ∂V⟶∂t+ρV⟶⋅∇V⟶=−∇P+∇⋅V⟶∇V⟶+∇V⟶T+ρg+Fsv.The above equation is the convective transport equation for the volume functionc as well as the momentum conservation equation. Fsv is the surface tension term of the momentum equation.The CLSVOF model compensates the deficiencies of mass conservation and phase interface parameters of the level-set model and VOF model and combines the advantages of the two models. The CLSVOF model can combine the advantages of both models to capture the phase interface more accurately.


## 3. NSGA-II-Based Scheme Solution Steps

### 3.1. NSGA-II Optimization Process

NSGA-II is a multiobjective genetic algorithm using fast sorting and elite mechanism, which has the characteristics of fast operation speed and good convergence [16, 17] and can effectively solve the multiobjective optimization problem in this paper. The Pareto solution set can be obtained by NSGA-II, and the solution steps are as follows.(1)
Input the parameters of the distribution system and generate the derivative matrix.(2)
Initialize the algorithm parameters, population numberp0, crossover parameter fga, variation parameter Cga, and maximum number of iterations N.(3)
Randomly create the initial populationp0 and calculate the objective function value for each individual.(4)
Sort the populations according to nondominance and calculate the crowding factor, and each solution is assigned a fitness value equal to its own nondominance rank.(5)
Use a binary tournament mechanism to selectp0/2 individuals as the parent population and perform genetic inheritance, crossover, and mutation to produce p0/2 offspring populations.(6)
Invoke the trend calculation program to calculate the objective function values of the offspring populations.(7)
Merge the resulting offspring populations with the original parent populations, rank the nondominated cases, and calculate the crowding coefficients.(8)
Selectp0 dominant individuals as the next generation population based on hierarchical analysis, rank nondominance, and calculate crowding coefficients.(9)
Determine whether the maximum number of iterations N is reached; if yes, output the Pareto solution set; otherwise, proceed to (4) loop.


### 3.2. Selection of the Multiobjective Optimal Solution

The final optimization result of NSGA-II is a set of Pareto solution sets from which the decision maker needs to select the optimal solution based on the preference information, which is essentially a multiattribute decision problem. In this paper, we adopt the approximate ideal solution ranking method (TOPSIS) to perform the optimal solution selection.The approximate ideal solution ranking method (TOPSIS) [18], in essence, is a decision-making method by calculating the distance between the alternative solution and the ideal solution and the negative ideal solution, so that the selected solution has the minimum distance from the ideal solution and the maximum distance from the negative ideal solution. Firstly, all the noninferior solutions in the Pareto solution set are obtained to form N alternative solutions x1,x2,x3,…,xN, and the number of attributes of the solution is n, i.e., the number of objective functions; then, the mth attribute value of solution xi is fmxi. First, the attribute values of all solutions are converted to dimensionless attributes, and the mth attribute value of the processed solution xi is(4)fm′xi=fmxi∑mNfm2xi\∗MERGEFORMAT.The relative distancedxi for scenario xi is calculated as follows:(5)dxi=d+xid+xi+d.xi\∗MERGEFORMAT,d+xi=∑m=1nλmfm′xi−λmfm+′2\MERGEFORMAT,d−xi=∑m=1nλmfm′xi−λmfm−′2\∗MERGEFORMAT,where d+xi、 d−xi represent the distance from solution xi to the ideal solution and the negative ideal solution; λm is the weight of attribute fm′, which is taken as equal weight in this paper; and fm−′、 fm+′ represent the optimal and inferior values of all solutions after normalization of attributes. The solution with the smallest relative distance among all solutions is the optimal solution in the multiobjective optimization.


## 3.1. NSGA-II Optimization Process

NSGA-II is a multiobjective genetic algorithm using fast sorting and elite mechanism, which has the characteristics of fast operation speed and good convergence [16, 17] and can effectively solve the multiobjective optimization problem in this paper. The Pareto solution set can be obtained by NSGA-II, and the solution steps are as follows.(1)
Input the parameters of the distribution system and generate the derivative matrix.(2)
Initialize the algorithm parameters, population numberp0, crossover parameter fga, variation parameter Cga, and maximum number of iterations N.(3)
Randomly create the initial populationp0 and calculate the objective function value for each individual.(4)
Sort the populations according to nondominance and calculate the crowding factor, and each solution is assigned a fitness value equal to its own nondominance rank.(5)
Use a binary tournament mechanism to selectp0/2 individuals as the parent population and perform genetic inheritance, crossover, and mutation to produce p0/2 offspring populations.(6)
Invoke the trend calculation program to calculate the objective function values of the offspring populations.(7)
Merge the resulting offspring populations with the original parent populations, rank the nondominated cases, and calculate the crowding coefficients.(8)
Selectp0 dominant individuals as the next generation population based on hierarchical analysis, rank nondominance, and calculate crowding coefficients.(9)
Determine whether the maximum number of iterations N is reached; if yes, output the Pareto solution set; otherwise, proceed to (4) loop.


## 3.2. Selection of the Multiobjective Optimal Solution

The final optimization result of NSGA-II is a set of Pareto solution sets from which the decision maker needs to select the optimal solution based on the preference information, which is essentially a multiattribute decision problem. In this paper, we adopt the approximate ideal solution ranking method (TOPSIS) to perform the optimal solution selection.The approximate ideal solution ranking method (TOPSIS) [18], in essence, is a decision-making method by calculating the distance between the alternative solution and the ideal solution and the negative ideal solution, so that the selected solution has the minimum distance from the ideal solution and the maximum distance from the negative ideal solution. Firstly, all the noninferior solutions in the Pareto solution set are obtained to form N alternative solutions x1,x2,x3,…,xN, and the number of attributes of the solution is n, i.e., the number of objective functions; then, the mth attribute value of solution xi is fmxi. First, the attribute values of all solutions are converted to dimensionless attributes, and the mth attribute value of the processed solution xi is(4)fm′xi=fmxi∑mNfm2xi\∗MERGEFORMAT.The relative distancedxi for scenario xi is calculated as follows:(5)dxi=d+xid+xi+d.xi\∗MERGEFORMAT,d+xi=∑m=1nλmfm′xi−λmfm+′2\MERGEFORMAT,d−xi=∑m=1nλmfm′xi−λmfm−′2\∗MERGEFORMAT,where d+xi、 d−xi represent the distance from solution xi to the ideal solution and the negative ideal solution; λm is the weight of attribute fm′, which is taken as equal weight in this paper; and fm−′、 fm+′ represent the optimal and inferior values of all solutions after normalization of attributes. The solution with the smallest relative distance among all solutions is the optimal solution in the multiobjective optimization.


## 4. Microstructure Model

For nonlinear objective functions, the ability of the proxy model to fit is particularly impressive and has great potential for application in various research areas [7].After obtaining the sample datasets,ys, the kriging model will perform a linear weighted interpolation for the sample function response values, i.e.,(6)y^x=∑i=1nλiyi.λ is the weighting factor, which is the most critical influence parameter in the agent model building process. The optimal weighting factor in the kriging model should satisfy the mean squared error:(7)MSEy^x=EλTYS−Yx2.Minimum,x and y the mathematical expectation relation should satisfy(8)E∑i=1nλiYxi=EYx.Y function is defined as(9)Y=β0+Zx,where β0 is the mathematical expectation of Yx and Zx is the static random process with mean 0 and variance σ2=σ2x.The weighting factorλ is obtained and substituted into the y^x function to predict the value of the function at any x, thus greatly reducing the computational effort of numerical optimization.Figure1 shows the schematic diagram of the three microchannel cells of the microreactor, which is also the flow region of the multiphase fluid. The microchannel size and baffle size are 2.4 mm, the inlet and outlet widths are 5 mm and 3 mm, respectively, and each cell is connected by a narrow throat with a grid number of 120,000 after irrelevance verification. The geometrical dimensions of the throat width and baffle length are the research objectives to be further optimized.Figure 1
Schematic diagram of geometric structure.In order to quantitatively compare the effects of different geometries on the gas-liquid mixing effect and then find the optimal microchannel structure, the boundary coefficients are defined here:(10)σS=AσS×100%.Aσ is the area of the region between 0.1 and 0.9 for the volume fraction of the gas phase, and the larger the area of this part, the more adequate the intersection of gas and liquid phases. s is the total area of the three microchannel cells. The magnitude of the boundary coefficient is the ratio of the gas-liquid phase interface area to the total area when the fluid is flowing inside the microreactor, which can measure the mixing effect of the gas-liquid phases in the microchannels.


## 5. Microstructure Model Optimization Effect

Figure2 shows the gas-liquid phase distribution when calculating different widths and lengths at the same time, and the boundary coefficients are extracted using postprocessing software to compare their gas-liquid phase mixing effects.Figure 2
Gas-liquid phase distribution for different throat widths and baffle lengths att = 150 ms.By changing the throat size and baffle length, the geometry is optimized for the throat area and the channel area by changing the ratio of the area in front of the throat and the channel area.The improvement of these two structures is achieved by proxy model optimization. The structure of the proxy model for this problem can be summarized as follows: the maximum value of the relationship between the design variables throat sizex1, baffle length x2 and the objective function of the function y boundary coefficient is solved. The objective function of the design variables x1 and x2 with respect to the boundary coefficients is refined by collecting a certain number of simulation data points, and the optimal throat size and baffle length geometry is accurately predicted by this objective function.Figure3 shows the sensitivity analysis results of the two independent variables of throat size and baffle length. The blue part shows the results of the sensitivity analysis of the throat size, and the yellow part shows the percentage of sensitivity of the baffle length. The left side shows the global sensitivity analysis, i.e., in addition to these two independent variables, the effects of all the remaining unknown independent variables are also considered. In this analysis, the throat size and baffle length are the most important influencing factors of the boundary coefficient, which together account for 82% of the influencing factors, while the remaining variables account for only 18%. The reason for this result is that the throat size affects the flow of gas-phase fluid into the mixing enhancement domain through the throat area, while the baffle length affects the percentage of the channel area in the overall microreactor structure, with the baffle size remaining constant. As it was pointed out in the previous paper, a smaller throat area and a smaller channel area will be more favorable for the increase of the boundary coefficient, and these two variables are the main influencing variables for the throat area and the channel area, which will have a greater impact on the boundary coefficient than the other related variables.Figure 3
Sensitivity analysis results.The image on the right shows the comparison between the effect of throat size and baffle length only. The sensitivity of the throat size is 45%, and the sensitivity of the baffle length is 55%. It can be seen that the baffle angle has a greater influence on the boundary coefficient. This is also reflected in the data; as the baffle angle changes, the boundary coefficient will have a more obvious trend. The reason for this result is that the length of the baffle is the most direct variable affecting the area share of the access area, and the area share of the access area and the length of the baffle basically show a linear positive correlation growth relationship. The size of the throat has more influence on the flow of the gas phase from the throat into the mixing domain than the geometry, and its change does not produce significant geometrical changes on the area share of the throat and mixing domain.Figure4 shows the 3D kriging interpolation results for the three variables. The highest boundary coefficient point extracted is (1.4, 8), i.e., the throat size is 1.4 mm and the baffle length is 8 mm. This point is the boundary point of the two design variables, and it can be seen from the three-dimensional coordinate diagram that the throat size and baffle length show a negative linear relationship within the value range, and the boundary coefficient at the lower limit of the two obtains the highest value. There are two main reasons for this conclusion [21].(1)
The narrow throat size makes it easier for the gas flowing through the throat to separate from a larger bubble cluster into a group of small bubble clusters, which in turn forms a better gas-liquid mixing effect at the front of the baffle and increases the value of the boundary coefficient.(2)
The length of the baffle plate most directly reduces the area share of the channel area. The boundary coefficient in the channel zone is lower than that in the baffle front and corner zone, and the baffle length is the most direct factor affecting the area of the channel zone; therefore, the reduction of the baffle length will increase the boundary coefficient inside the microreactor very significantly.Figure 4
Three-dimensional kriging interpolation.The throat and baffle front have an important role in enhancing the gas-liquid phase mixing effect. Figure5 shows the results of numerical simulations and high-speed camera images. From Figure 5, it can be seen that the inside of the throat and the front of the baffle are the high-pressure areas of the microchannel unit, and the highest pressure is about 8000 Pa, which can greatly increase the mixing effect of the two phases. The gas phase will collide with the baffle wall after flowing through the throat, and the gas momentum direction will change after colliding with the baffle wall, and the large bubbles will be broken into smaller ones. The other part of the bubbles collides with the wall and rejoins with the new influx of fluid at the throat. This is the most important reason why the gas-liquid junction area in this region is significantly higher than the rest of the region.Figure 5
Experimental results and simulation results. (a) Experimental results. (b) Microchannel gas-liquid phase distribution diagram.


## 6. Conclusions

In this paper, a full linear regression of the SAR microchannel is performed using a proxy model. A series of numerical simulations using the CLSVOF multiphase flow model and the SST k-ω shear pressure transport turbulence model demonstrate outstanding computational accuracy in calculating the gas-liquid two-phase distribution and phase interface capture. The advantages and disadvantages of the level-set model, VOF model, and CLSVOF model are investigated and analyzed for the different multiphase flow models, i.e., the level-set model, the VOF model, and the CLSVOF model. Both of them have different degrees of computational drawbacks when applied individually. Therefore, the CLSVOF model, which combines the two computational models, has a greater computational advantage over the above two models for the study in this paper, which focuses on phase interface capture.



---

*Source: 6368018-2022-04-07.xml*