Modeling and Optimization of Gaseous Thermal Slip Flow in Rectangular Microducts Using a Particle Swarm Optimization Algorithm

In this study, pressure-driven flow in the slip regime is investigated in rectangular microducts. In this regime, the Knudsen number lies between 0.001 and 0.1. The duct aspect ratio is taken as 0 ≤ ε ≤ 1. Rarefaction effects are introduced through the boundary conditions. The dimensionless governing equations are solved numerically using MAPLE and MATLAB is used for artificial neural network modeling. Using a MAPLE numerical solution, the shear stress and heat transfer rate are obtained. The numerical solution can be validated for the special cases when there is no slip (continuum flow), ε = 0 (parallel plates) and ε = 1 (square microducts). An artificial neural network is used to develop separate models for the shear stress and heat transfer rate. Both physical quantities are optimized using a particle swarm optimization algorithm. Using these results, the optimum values of both physical quantities are obtained in the slip regime. It is shown that the optimal values ensue for the square microducts at the beginning of the slip regime.


Rarefied Gas Flows
Flows in microducts are found in microelectromechanical systems, nanotechnology applications, therapeutic and superhydrophobic microchannels, low-pressure environments, biochemical applications and cryogenics.The rarefaction effect can be found in microchannels and can be expressed in terms of the Knudsen number.In this case, the deviations from continuum behavior are smaller.Therefore, the Navier-Stokes equations can be employed with slip boundary conditions.The difference between a fully developed flow in a rectangular duct and in a microchannel is that a rectangular microduct needs 2D analysis, whereas a microchannel entail a one-dimensional analysis.
Generally, the Navier-Stokes equations can be solved with slip boundary conditions at the microduct walls [1,2].Previous studies [1,2] found that the solutions were validated by the experimental data in several microchannel flows.Otherwise, significant departures were observed between the numerical and experimental results without applying the slip conditions.Ameel et al. [3] confirmed the statement of Reference [1,2].Zade et al. [4] considered variable physical properties and investigated thermal features of developing flows in rectangular microchannels.They employed a collocated finite volume method and studied several channel dimensions for different values of slip parameters.They found that variable properties show substantial changes in flow behavior.Ghodoossi and E grican [5] determined the heat transfer rate using an integral method.They predicted thermal features for isothermal boundary conditions.They noticed that, for any aspect ratio, the rarefaction had obvious effects on the heat transfer.
The effects of rarefaction on the Nusselt numbers were investigated by Hettiarachchi et al. [6] for thermally developing flows.They noticed that these Nusselt numbers increased with the slip velocity, and decline with the thermal jump boundary condition.The slip flow forced convection in rectangular microchannels was also considered by Yu and Ameel [7].It was found that the velocity slip and aspect ratio enhance heat transfer.For rarefied gas flows, Renksizbulut et al. [8] observed exceptionally large reductions in the friction and heat transfer coefficients.They found that these coefficients are insensible to rarefaction effects in corner-dominated flows.Tamayol and Hooman [9] analyzed microchannels of polygonal, rectangular, and rhombic cross-sections and showed that the Poiseuille numbers decrease with increasing aspect ratios and Knudsen numbers for rectangular channels.Hooman [10] investigated forced convection in microducts using a superposition approach for both thermal boundary conditions and demonstrated that their results were valid for several cross-sections.Sadeghi and Saidi [11] considered the rarefaction effects in two different microgeometries and demonstrated that for both geometries, the Brinkman number reduced the heat transfer rate.In fact, the Brinkman number measures the viscous heating effects with reference to the heat conduction.They also developed a correlation for the Nusselt number.
Yovanovich and Khan [12] developed several slip flow models for microchannels of different cross-sections.The slip flow in circular and other cross-sections was investigated by a number of authors, including References [13][14][15][16].They developed different models for envisaging the friction factor under both developing and fully developed conditions.Duan and Muzychka [13,14] proposed a simple model to predict the Poiseuille numbers in these microchannels.Duan and Muzychka [15] considered slip flow in the entrance of circular and parallel plate microchannels and noticed that the Poisseuille number depends on the Knudsen number.Yovanovich and Khan [16] modeled the Poiseuille flow in long channels of different cross-sections.Ebert and Sparrow [17] analyzed abstemiously rarefied gas flows in different ducts and noticed that the velocity flattened and the friction lessened due to the slip parameter.
Baghani et al. [18] treated rarefaction effects in microducts of different cross-sections and obtained Nusselt numbers under isothermal boundary condition.They tabulated Nusselt numbers for the entire slip flow range of the Knudsen number.Wang [19] considered the slip flow under an isothermal boundary condition in different ducts.It was found that both hydrodynamic and thermal boundary conditions have substantial effects on both the Poiseuille and Nusselt numbers.Also, both parameters showed an opposite trend with an increasing velocity slip.Beskok and Karniadakis [20] considered internal rarefied gas flows and presented mathematical models for different regimes.To account for heat transfer from the surface, they suggested a new boundary condition and demonstrated that it is applicable in all flow regimes.Yovanovich and Khan [21] proposed compact models for isothermal, incompressible slip flows in different microchannels.They used the principle of superposition in the analysis and introduced new flow parameters.They compared their results with previous numerical results for different cross-sections and found them in a very good agreement.Klášterka et al. [22] obtained analytical and numerical solutions for the flow variables and demonstrated the effect of slip parameters on the friction and heat transfer coefficients.
The above literature survey reveals that there are only a few studies related to rectangular microducts.However, there has been no study where the numerical models are developed using an artificial neural network (ANN) and are optimized for the duct geometry.This fact inspires us to model our numerical results using ANN and to acquire optimized values of friction factors and Nusselt numbers using a particle swarm optimization algorithm.

Artificial Neural Networks
Artificial neural networks (ANNs) are mathematical models, which are inspired by biological nervous systems [23,24].The multilayer perceptron neural network (MLPNN) is one of the widely known feedforward neural networks.MLPNNs have three types of layers, which are the input layer, the hidden layers, and the output layer.Figure 1 shows the structure of MLPNN [25,26].

Artificial Neural Networks
Artificial neural networks (ANNs) are mathematical models, which are inspired by biological nervous systems [23,24].The multilayer perceptron neural network (MLPNN) is one of the widely known feedforward neural networks.MLPNNs have three types of layers, which are the input layer, the hidden layers, and the output layer.Figure 1 shows the structure of MLPNN [25,26].The activation function is defined as an output of a node or a neuron and can be written as: where x is an input value from the input layer.Supervised learning of ANNs is useful for improving the performance of ANN models.Several algorithms are used for the training to find the best parameters of the neural networks [27,28].The optimization algorithms are used to find the best values for inputs and output weights that are parameters of MLPNNs.In this study, the particle swarm optimization (PSO) algorithm is used to train the neural networks, because it is one of the most effective meta-heuristic algorithms used for optimization problems.Equation ( 2) is the sum of squared error (SSE) function, which was used as a fitness function for testing the performance of the ANNs [25].

SEE= ∑(𝐴𝑐𝑡𝑢𝑎𝑙 𝑜𝑢𝑡𝑝𝑢𝑡 − 𝑡𝑎𝑟𝑔𝑒𝑡 𝑜𝑢𝑡𝑝𝑢𝑡)
( where the   values are determined by Equation (3), while the   values are determined by the neural network (NN) output values.
where n, m, k,  , and ′ represent the number of neurons in the hidden layer, number of input data, number of neurons in input layer, number of neurons in output layer, input weight-which is between the input neurons and hidden neurons-and output weight, which is between the hidden neurons and output neurons, respectively.

Particle Swarm Optimization Algorithm
The principle of the particle swarm optimization (PSO) algorithm came from the ideas of the social behavior of bird flocks, as well as from the shoaling behavior of fish [27][28][29].All solution The activation function is defined as an output of a node or a neuron and can be written as: where x is an input value from the input layer.Supervised learning of ANNs is useful for improving the performance of ANN models.Several algorithms are used for the training to find the best parameters of the neural networks [27,28].The optimization algorithms are used to find the best values for inputs and output weights that are parameters of MLPNNs.In this study, the particle swarm optimization (PSO) algorithm is used to train the neural networks, because it is one of the most effective meta-heuristic algorithms used for optimization problems.Equation ( 2) is the sum of squared error (SSE) function, which was used as a fitness function for testing the performance of the ANNs [25].
where the target output values are determined by Equation (3), while the Actual output values are determined by the neural network (NN) output values.
where n, m, k, w jk , and w ki represent the number of neurons in the hidden layer, number of input data, number of neurons in input layer, number of neurons in output layer, input weight-which is between the input neurons and hidden neurons-and output weight, which is between the hidden neurons and output neurons, respectively.

Particle Swarm Optimization Algorithm
The principle of the particle swarm optimization (PSO) algorithm came from the ideas of the social behavior of bird flocks, as well as from the shoaling behavior of fish [27][28][29].All solution members of the algorithm are called particles.The particle flies by updating its velocity vectors v and position x, by using Equations ( 4) and (5), respectively [30].
where x i (t), v i (t), y i (t), and ŷk (t) represent the position at time t, the velocity at time t, the best personal position (p best ) at time t, and the global best position (g best ) of population at time t respectively, whereas c 1 and c 2 are the learning factors of (p best ) in interval between 0 and 2 and r 1 and r 2 are the random numbers in the interval 0 to 1.
It is important to note that the other particles follow the performance of the best particle [31][32][33][34].In addition, each particle keeps the best position that it has achieved so far.The flowchart of Figure 2 shows the steps of the present PSO algorithm [27,32,35].
Symmetry 2018, 10, x FOR PEER REVIEW 4 of 14 members of the algorithm are called particles.The particle flies by updating its velocity vectors v and position x, by using Equation (4) and Equation ( 5), respectively [30].
( + 1) =  () +  ( + 1) where  (),  (),  (), and  ∧ () represent the position at time t, the velocity at time t, the best personal position (  ) at time t, and the global best position  ) of population at time t respectively, whereas  and  are the learning factors of  ) in interval between 0 and 2 and r1 and r2 are the random numbers in the interval 0 to 1.
It is important to note that the other particles follow the performance of the best particle [31][32][33][34].In addition, each particle keeps the best position that it has achieved so far.The flowchart of Figure 2 shows the steps of the present PSO algorithm [27,32,35].Thus, the particles should move towards the best directions, and utilize useful information from other particles along with the best particles [30,35].
In this study, the algorithm was used to find the optimal values of the Nusselt number (Nu) and Poiseuille number (Po) corresponding to the aspect ratio and the Knudsen number in the slip flow regime.The number of the iterations was set as 50, with initial populations 60, and c1 = c2 = 0.05 [30].
To the best of our knowledge, most of the existing analytical solutions for slip flows are limited to fluid flow in rectangular ducts.There is a lack of literature related to heat transfer in microducts.Furthermore, there is no study related to ANN modeling and optimization of pressure drop and heat transfer in microducts.The aim of this paper was to cover these aspects.

Hydrodynamic Analysis
Let us consider a fully developed flow in a rectangular duct, as shown in Figure 3.The semimajor and -minor axes are taken as  and  along the x-and y-axes, respectively.The aspect ratio is / b a ε = .The gas is assumed to flow along the z-axis.The momentum equation for the flow in a microduct can be written as: Thus, the particles should move towards the best directions, and utilize useful information from other particles along with the best particles [30,35].
In this study, the algorithm was used to find the optimal values of the Nusselt number (Nu) and Poiseuille number (Po) corresponding to the aspect ratio and the Knudsen number in the slip flow regime.The number of the iterations was set as 50, with initial populations 60, and c 1 = c 2 = 0.05 [30].
To the best of our knowledge, most of the existing analytical solutions for slip flows are limited to fluid flow in rectangular ducts.There is a lack of literature related to heat transfer in microducts.Furthermore, there is no study related to ANN modeling and optimization of pressure drop and heat transfer in microducts.The aim of this paper was to cover these aspects.

Hydrodynamic Analysis
Let us consider a fully developed flow in a rectangular duct, as shown in Figure 3.The semi-major and -minor axes are taken as a and b along the x-and y-axes, respectively.The aspect ratio is ε = b/a.The gas is assumed to flow along the z-axis.The momentum equation for the flow in a microduct can be written as: which is a 2D equation for the axial velocity u, uniform gas viscosity µ and constant pressure gradient dP/dz.Due to double symmetry, only the OABCO is selected for the simulation.The velocity slip and symmetry boundary conditions for this problem can be written as: and: where σ v is the tangential-momentum accommodation coefficient, which lies between 0.9 and 1, and λ is the mean free path.If A is the area of cross-section and P is the perimeter of the duct, then the hydraulic diameter, D h can be written as: Symmetry 2018, 10, x FOR PEER REVIEW 5 of 14 which is a 2D equation for the axial velocity u , uniform gas viscosity μ and constant pressure gradient / dP dz .Due to double symmetry, only the OABCO is selected for the simulation.The velocity slip and symmetry boundary conditions for this problem can be written as: and: where v σ is the tangential-momentum accommodation coefficient, which lies between 0.9 and 1, and λ is the mean free path.If A is the area of cross-section and P is the perimeter of the duct, then the hydraulic diameter, h D can be written as: At this stage, it is convenient to define the Knudsen number (the measure of rarefaction effects) as follows: Following Ebert and Sparrow [17], Equations (6-8) can be non-dimensionalized as follows.Let: The dimensionless momentum equation can be written as: where the axial velocity u is normalized with  = | |.Using Equation (11), dimensionless boundary conditions can be written as: At this stage, it is convenient to define the Knudsen number (the measure of rarefaction effects) as follows: Following Ebert and Sparrow [17], Equations ( 6)-( 8) can be non-dimensionalized as follows. Let: The dimensionless momentum equation can be written as: where the axial velocity u is normalized with u o = b 2 µ dP dz .Using Equation (11), dimensionless boundary conditions can be written as: Symmetry 2019, 11, 488 6 of 13 The mean velocity in the duct is defined as: The Poiseuille number Po is defined as: where c f = τ w /0.5ρu 2 is the coefficient of friction and Re D h = uD h /ν is the Reynolds number, based on mean velocity.For a fully developed flow, the force balance on a control volume of width dz gives: Thus, the dimensionless shear stress can be expressed in terms of the Poiseuille number Po: where u/u 0 is the dimensionless average axial velocity in the microduct and can be determined from Equation ( 14).

Thermal Analysis
For fully thermally developed and steady flow in the rectangular microducts, the energy equation is given by: with temperature jump and symmetry boundary conditions: ∂T ∂x x=a (19) and: where Pr is the Prandtl number for the selected gas, σ t is the energy accommodation coefficient, and γ is the ratio of specific heats of gas.Using Equation ( 11), the energy balance in the flow direction can be written as: or: where Q is the total energy per unit length of the duct and is assumed to be constant.Using θ = (T − Tw)/(Q/k f ), Equation ( 22) can be written as: Symmetry 2019, 11, 488 7 of 13 with transformed boundary conditions [17]: For a constant temperature jump condition, the Nusselt number can be defined as:

Solution Procedure
The solution procedure comprises three steps.In the first step, the numerical values of the Nu and Po numbers are determined, corresponding to different values of ε and Kn by solving the governing Equations ( 12) and ( 23) by a finite difference method using MAPLE, and numerical data was obtained.In the second step, the back-error-propagation algorithm is employed to train a multi-layer perceptron neural network (MLPNN) using the neural network toolbox in MATLAB.There are two inputs, namely, the aspect ratio ε and the Knudsen number Kn.The corresponding outputs are the predicted values of the Nu and Po numbers, where each output is represented by a separate model.In each model, different input, hidden and output neurons are used.Finally, the results are optimized using PSO algorithm in order to find the optimal values of Nu and Po without gaining to the numerical values.Accordingly, the dimensionless partial differential equations will not use again in any stage of our work.As a result, the PSO algorithm uses neural network models to obtain the optimum values of Nu and Po.

Results and Discussion
The numerical values of Poiseuille numbers for forced convection in microducts were compared with the available data for different values of Kn in Tables 1-3.The comparison showed excellent agreement with the available data.The predicted values of the Nusselt numbers are also provided in Table 4 for different values of ε and Kn.To build the NN models, 70% of data was used as the training data, 15% as the testing data, and 15% as the validation data.In NN MATLAB tools, the back-error-propagation algorithm (an optimization algorithm) was used to determine the best models that have the optimal errors.These models were based on the mean squared error function of MSE = 0.00020768 and MSE = 0.00022112 for the two models, respectively, as shown in Figures 4 and 5.The figures show that the best validation performance, in terms of the mean squared error function, were found after 187 and 452 iterations, respectively.The correlation coefficients of the models are shown in Figures 6 and 7.As expected, the correlation coefficient of the models was greater than 0.95.In addition, Figures 6 and 7 show a significant convergence between the real values and the corresponding values that resulted from models.Accordingly, there was enough support to believe that the models were suitable to find the estimated values of the Nu and Po numbers.In addition, the coefficient of determination of the models was more than 0.99, which meant that both models could represent whole target data.
For optimization, the PSO algorithm was used to find the optimal values of Po and Nu, based on NN models.The optimal values of Po and Nu were achieved when ε = 1, and Kn = 0.001.The behavior of Po and Nu values in the PSO algorithm is shown in Figures 8 and 9.The optimal values of the curves were determined when the curves became stable.Accordingly, the optimal values for Nu and Po were found in iteration number 4 and 10, respectively.Note that, the data set values were converted into z-score values for use in the neural network toolbox in MATLAB.models was more than 0.99, which meant that both models could represent whole target data.
For optimization, the PSO algorithm was used to find the optimal values of Po and Nu, based on NN models.The optimal values of Po and Nu were achieved when 1 ε = , and 0.001 Kn = . The behavior of Po and Nu values in the PSO algorithm is shown in Figures 8 and 9.The optimal values of the curves were determined when the curves became stable.Accordingly, the optimal values for Nu and Po were found in iteration number 4 and 10, respectively.Note that, the data set values were converted into z-score values for use in the neural network toolbox in MATLAB.

Conclusions
In this study, the slip flow and heat transfer were investigated in rectangular microducts.The dimensionless governing partial differential equations with slip boundary conditions were solved numerically using MATLAB.The numerical data were used successfully for building ANN models.The outputs of the models were the corresponding Poiseuille Po and Nusselt Nu numbers to the input values  and .The main outcomes of the study were: 1.The small tolerance values for the mean squared error function values and a very close correlation coefficient of 1 in the NN models proved that the NN models were suitable to use for predicting the Poiseuille Po and Nusselt Nu numbers.2. Without going back to the numerical data, the proposed models were used with the PSO

Conclusions
In this study, the slip flow and heat transfer were investigated in rectangular microducts.The dimensionless governing partial differential equations with slip boundary conditions were solved numerically using MATLAB.The numerical data were used successfully for building ANN models.The outputs of the models were the corresponding Poiseuille Po and Nusselt Nu numbers to the input values ε and Kn.The main outcomes of the study were:

Figure 6 .Figure 6 . 14 Figure 7 .Figure 7 . 13 Figure 7 .
Figure 6.Correlation coefficient and the regression of the NN model values of Nu and the desired values.

Figure 8 .Figure 8 .
Figure 8. Best performance of the PSO algorithm in terms of Nu.

Figure 9 .
Figure 9. Best performance of the PSO algorithm in terms of Po.

Figure 9 .
Figure 9. Best performance of the PSO algorithm in terms of Po.

Table 1 .
Comparison of the Poiseuille numbers for forced convection when Kn = 0.1.

Table 2 .
Comparison of numerical and exact values of Po when Kn = 0.001.

Table 3 .
Comparison of numerical and exact values of Po when Kn = 0.1.

Table 4 .
Numerical values of Nu for different values of Kn.