STUDY OF THE FLOW AND HEAT TRANSFER OF A VISCOELASTIC FLUID USING HYBRID NEURAL NETWORK-PARTICLE SWARM OPTIMIZATION (HNNPSO)

Fluid flow and heat transfer of a second-order viscoelastic fluid in an axisymmetric channel with a porous wall for turbine cooling applications are studied. The nonlinear differential equations of the fluid flow and heat transfer arising from similarity solutions are computed employing a Hybrid Neural Network-Particle Swarm Optimization algorithm (HNNPSO). A trial function, satisfying the boundary conditions, as a possible solution for the governing equations is introduced. The trial functions incorporate a multi-layer perceptron neural network with adjustable parameters (the weights and biases). The Particle Swarm Optimization algorithm (PSO) is applied to find the adjustable parameters of the trial solution to satisfy the governing equations. Finally, comparisons are made between the results of the present method (HNNPSO) and the results of the fourth order Runge–Kutta method, finite difference method, and Variational Iteration Method. The results indicate that HNNPSO method conveniently produces a polynomial analytic solution with remarkable accuracy, and the accuracy of the solution improves as the number of neurons of the neural network increases.


INTRODUCTION
Viscoelastic fluids are a kind of non-Newtonian fluids which have both viscous and elastic characteristics. In these materials, the viscosity decreases as the shear/strain rate remains constant. Polymer solutions and melts like nylon, gelatin solution, and most of the industrial fluids are viscoelastic in nature [1].
Analysis of the behavior of the non-Newtonian fluid flows, especially the viscoelastic flows, have attracted much attention in the recent years due to the various industrial applications in different fields such as extrusion of plastics, lubrication and hot rolling. Generally, the interest in flow and heat transfer problem of non-Newtonian fluids has grown considerably [1,2]. The boundary layer flow and heat transfer over surfaces is an important issue with many practical engineering applications. Various aspect of flow and heat transfer over surfaces of channels have been addressed in recent literature works. For instance, the boundary layer heat transfer of viscoelastic nanofluids [3], non-Newtonian nano-fluids [4], micropolar fluids [5], and radiation effects [6] have been addressed recently.
In most cases, the governing equations of the fluid motion in the boundary layer are inherently nonlinear in nature that cannot be solved analytically, therefore, special techniques such as numerical techniques or advanced analytical methods to tackle such nonlinear equations are needed. In this way, different methods such as Homotopy Perturbation Method (HPM) [7], Parameterized Perturbation Method (PPM) [8], Homotopy Analysis Method (HAM) [9], Adomian Decomposition Method (ADM) [10] and Variational Iteration Method (VIM) [11] have been developed to solve differential equation. Although these methods are very powerful and can deal with various types of nonlinear differential equations; they have also encountered some limitations and disadvantages. For example, the Homotopy Perturbation Method (HPM) and Parameter Perturbation Method (PPM) depend on perturbation, linearization, or small parameters. The Homotopy Analysis Method (HAM) requires many harmonic terms to produce an accurate solution, which such a solution with too many analytic terms practically seems unworthy. VIM shows some disadvantages, namely, repeated computations and computations of unneeded terms, which consumes time and effort. Hence, a simple and powerful method to solve the nonlinear equations, which reduces these limitations and difficulties, is still a challenge.
Many researchers have studied the flow and heat transfer phenomena over open surfaces and channels. For instance, the MHD flows [12,13,14], flow over flat surfaces [15], natural and mixed convection flow and heat transfer [16][17][18], flow over a cylinder [19], and entropy generation [20] have been studied in past years. The non-Newtonian fluids are important in industrial applications as the behavior of many practical fluids cannot be described using the Newtonian model of fluids. In this regard, many researchers have studied the flow and heat transfer of non-Newtonian fluids. Among them, White and Metzner [21] have evaluated the constitutive equations for viscoelastic fluids. Debruge and Han [22] studied the heat transfer in a channel with a porous wall for turbine cooling application. Following the study of Debruge and Han [22], Kurtcebe and Erim [23] used the Perturbation Method (PM), and Shakeri et al. [24] applied Variational Iteration Method (VIM) to solve the heat transfer of non-Newtonian viscoelastic fluid flow in an axisymmetric channel with a porous wall. They [23,24] analyzed the effects of different parameters namely Reynolds number, injection Reynolds number, Prandtl number, and power-law index. The same problem was also solved by Esmaeilpour et al. [25] using the Homotopy Analysis Method (HAM) to obtain an expression for velocity and temperature fields.
The intelligence exhibited by machines or software is known as artificial intelligence. An Artificial Neural Network (ANN) is a new form of artificial intelligence used in machine learning that tries to mimic the learning process of the human brain. Artificial neural networks are powerful tools for functions approximation. This method has been successfully applied to a wide variety of applications [26]. Particle swarm optimization (PSO) is a swarmbased artificial intelligence technique that tries to simulate the social behavior of biological populations [27]. This heuristic method can be used to find approximate solutions to numeric maximization and minimization problems.
In recent years, some new artificial methods have been introduced to solve differential equations arising from engineering problems [28][29][30][31][32][33][34]. Noghrehabadi et al. [30] studied the buckling shape of a multi-wall carbon nanotube cantilever. The buckling of the nanotube follows a fourth-order non-linear boundary value differential equation. They utilized the power series with adjustable variables as a trial function for the governing equation, satisfying the boundary conditions. Afterward, Artificial Bee Colony (ABC) algorithm was employed to adjust the power series variables to satisfy the governing differential equation in the domain of the solution. Using the power series as a trial function and ABC as the training method, Noghrehabadi et al. [30] have successfully reported a semianalytic solution in the form of a power series for the buckling shape of the nanotube.
In another study, Noghrehabadi et al. [31] successfully utilized the Cuckoo search optimization algorithm as the training method to adjust the coefficients of a trial function for the buckling of a microactuator. Lagaris et al. [32] and Tawfiq and Hussein [33] have utilized the feedforward neural networks to obtain a solution for differential equations. Malek and Beidokhti [34] examined the Feed Forward Neural Networks (FFNN) to obtain a solution for high order differential equations. They utilized the Nelder-Mead optimization method to adjust the parameters of a trial solution, which was constructed using the feedforward neural networks.
The approaches utilized in [28][29][30][31][32][33][34] can be considered as hybrid methods. An advantage of the hybrid method is that it can provide an analytical approximation to a wide class of nonlinear equations without the computational difficulties of HAM. Also, there is no need for any perturbation, linearization, or small parameter as required in perturbation methods. The hybrid method is simple and powerful to solve different differential equations, while the results have excellent agreement with numerical methods. In the previous studies [28][29][30][31][32][33][34], the hybrid method was applied for a single high order differential equation. Although there are various practical systems of differential equations, there has been no attempt to develop the hybrid method to the case of the system of high order differential equations. Recently, Noghrehabadi et al. [28] utilized the hybrid ANN and PSO optimization method to solve the differential equation of nanofluids over a flat plate.
As mentioned, the literature review shows that the problem of heat transfer of a non-Newtonian fluid in a channel has been studied by previous researchers [22][23][24][25] using Perturbation Method (PM) [22], Variational Iteration Method (VIM) [24] and Homotopy Analysis Method (HAM) [25] due to its important practical applications and inherent non-linearity. The governing equations of heat transfer of a non-Newtonian fluid in a channel result in a system of high order non-linear differential equations, one high order boundary value equation for heat and one for flow.
The utilized method in the literature requires as mall value [22], or too many terms [23][24][25] to provide an accurate analytical solution. The present study tends to obtain an accurate analytical solution with only a few terms for the model of heat transfer introduced in [22][23][24][25]. Here, the hybrid method is employed to deal with a system of high-order boundary value differential equations (a fourth-order differential equation and a second-order differential equation). For this purpose, two feed-forward neural networks are introduced as the trial functions, and the Particle Swarm Optimization (PSO) method as the training method is employed to obtain a semi-analytical solution for the heat transfer and fluid flow of a non-Newtonian fluid in a channel.

MATHEMATICAL FORMULATION
Consider the flow of a non-Newtonian viscoelastic fluid on a turbine disc to cool the surfaces of the disc. Fig. 1 shows the physical model of the problem and the selected coordinate system. The r-axis is adopted parallel to the surface of the disk, and the z-axis is normal to the surface of the disk. The porous disc of the channel is placed at z=+L. The wall that coincides with the r-axis is heated externally; a non-Newtonian fluid is injected uniformly from the other perforated wall. The flow and heat transfer equations and the similarity solution are proposed by Kurtcebe and Erim [23] and followed by Shakeri et al. [24] and Esmaeilpour et al. [25] for the above problem description. For a steady, axisymmetric and two-dimensional non-Newtonian fluid flow and neglecting dissipation effect, the governing equations for flow and heat can be represented as [23]: Where f and qn are the fluid and heat variables, respectively and η is the domain variable. Here, K1 is the viscoelastic parameter, Re is the injection Reynolds number, Pr is the Prandtl number, which are defined based on the first viscosity characteristic (φ0). n is the order of heat equation which can be considered as a constant integer positive number. The corresponding boundary conditions for flow and heat are: As seen, the momentum equation is a fourth order and highly non-linear equation, and the heat equation is a second order coupled to the momentum equation. To find an analytical solution, the hybrid neural network and swarm optimization method are introduced and utilized to solve these equations.

Artificial Neural Networks (ANNS)
An Artificial Neural Network (ANN) is a new mathematical model used in machine learning that tries to simulate the structure of biological neural networks and the learning process of the human brain. The basic element of every artificial neural network is an artificial neuron, which is a simple mathematical function [35]. Using the science of computer programming, one can design a data structure that acts as a neuron. Then, the desired network can be introduced and trained by establishing a network of interconnected artificial neurons and utilizing a learning algorithm. ANNs are known as universal approximations. These networks can efficiently estimate and approximate functions [36].
In cases that include too many input data, one neuron alone is not enough to solve the problem. Accordingly, there should be a layer of neurons. In a multi-layer neural network, there is an input layer, which receives the input data, some hidden layers, which receive the information from the previous layer. Ultimately, there is an output layer to deliver the outcomes. A wide variety of learning algorithms have been developed. A supervised network needs a teacher to tell the network what the desired output should be. In contrast, an unsupervised network adapts purely in response to its inputs. Such networks can learn to pick out a structure in their input.
In a feed-forward neural network, the neurons are in successive layers and send their output to the forwarded layers. Thus, there is no feedback, and the output layer does not affect the other layer. The topology of the neural network is an important factor. Multi-layer perceptron which first was introduced by Rosenblatt [37], is a simple feed-forward neural network that involves one input layer, one or more hidden layers, and one output layer. Fig.2 shows a typical three-layered perceptron.  The ability to function approximation is the most important application of a multi-layered perceptron. Based on the Kolmogorov theorem, a three-layered perceptron with n(2n+1) nodes can approximate a continuous function of n variables [38,39]. Hence, it can be concluded that the accuracy of the function approximation is under the significant influence of the number of neurons in the hidden layer.
Considering an input of x, the network output can be evaluated by Where wi denotes a weight parameter from the input layer to ith neuron in the hidden layer, bi denotes a bias parameter for ith neuron in the hidden layer, and vi denotes a weight parameter from ith neuron in the hidden layer to output layer [35]. Here, p is a vector containing wi, vi and bi.
The multi-layered feed-forward neural networks are trainable. In this paper, a feed-forward neural network is used to form a non-linear function from the input data. Here, the ANN learns the function behavior from the theory of differential equations. During the learning process the parameters of ANN as adjusted in a way to satisfy the differential equations in some collocation point samples. An ANN with minimum deviation from the behavior of the ODE is the final solution and has learned the ODE adequately. Therefore, it can be concluded that an ANN that has learned the ODE is an analytic function that approximates the solution to the corresponding ODE.

Particle Swarm Optimization (PSO)
The Particle Swarm Optimization (PSO) is a population-based optimization technique proposed by Eberhart and Kennedy [27], inspired by the social behavior of biological populations. Like other swarm-based techniques, PSO starts with a population of random solutions and searches for the optimum solution by updating generations. In PSO, each particle has a velocity and a position that can move. Each particle adjusts its movement in the search space according to two pieces of information, which are the best position the particle has reached so far, which is called the personal best (pbest), and the best position obtained by any particle in the neighbors of the particle, which is called the local best (lbest). When a particle takes the whole swarm as its neighbor, the value is called the global best (gbest). In the PSO, each particle i has a position Zi and a velocity vi that is updated at each iteration: Where ω is the inertia weight, ⃗ is the personal best found by particle, ⃗ is the global best, r1 and r2 are random values between [0,1], and c1 and c2 are positive constant parameters called acceleration coefficients which determine the effect of the cognitive and social behavior of the particles in the velocity vector. The position of each particle is the sum of the present position and the distance that the particles have traveled. Hence, the new positions of the particles can be evaluated as: Some important factors in PSO are inertia weight (ω), particle size and maximum iteration number (t). Generally, the PSO algorithm can be arranged in five key steps:  Initialize random velocity and position vectors for all particles. Then, maximum velocity and maximum particle movement amplitude are specified.  Evaluate the fitness of each particle.  If the obtained solution is better than the previous local best solution, replace it. If the existing solution is better than the previous global best solution, then the existing solution is selected as the new gbest.  Change positions and velocities by computing the new position employing Eqs. (6) and (7).  Repeat steps 2 to 4 until the gbest reaches a satisfactory optimized solution.

Hybrid Method For A System Of Differential Equations
In the hybrid method, mentioned in [30][31][32][33][34], a trial function with adjustable coefficients is introduced. The trial function consists of the sum of two parts. The first part is meant to satisfy the boundary conditions and the second part, which is an FFNN here, is meant to satisfy the differential equation in the domain of the solution. This technique has been referred to as the collocation neural network in the study of Tawfiq and Hussein [33]. In this section, we will propose the formulation of the hybrid method for a system of boundary value differential equations consisting of two non-linear differential equations. Now, assume the general form of a system of differential equations as ( , , … , , , ) = 0 (8) Subject to ( , , … , , , ) = 0 (9) Where z=z1, z2…zj are the governing differential equation, and y=y1, y2, … yj are the differential equation variables which are a function of the independent variable of η (the solution in domain D that should be computed). Here, b=b1, b2 … bL are the boundary equations. The indexes of j and L denote the number of differential equations and boundary conditions respectively. Consider a trial function, YT (η, a, p), with adjustable parameters of a and p, that exactly satisfies the boundary conditions introduced by b in Eq. (9). It should be noted that YT = YT1, YT2 … YTj and a and p are vectors of adjustable coefficients for each trial function of YT1, YT2 … YTj. Considering YT (η, a, p) as a trial solution with adjustable parameters of a and p, the problem is transformed into a discretized form as: In our proposed method the trial solution of YT employs a polynomial (A(η, a)), containing the adjustable coefficients of a to satisfy the boundary conditions and a function containing a feed-forward neural network (G(η,N(η, p)) as: A(η, a)= A1(η, a1), A2(η, a2)… Aj(η, aj) and G(η, N(η, p))= G1(η, N1(η, p1)), G2(η, N2(η, p2))… Gj(η, Nj(η, pj)) where each of a1, a2 … aj and p1, p2 … pj are vectors containing adjustable parameters. Each Ni(η, p), such as N1(η, p1) is a single layer feed-forward neural network the same as the one that was introduced in Eq.
Now, a minimization technique such as PSO, which was introduced in the previous section, is required to adjust the coefficients of a and p.
At this end, it is worth mentioning that the present formulation extends the ideas proposed in the previous studies of [28][29][30][31][32][33][34] to the case of a system of differential equations. In addition, in the studies of [28][29][30][31][32][33][34] the function of A(η) which satisfies the boundary conditions is solely a function of the independent variable of η. However, in the present study, the function A is considered as a function free of the independent variable of η and the adjustable coefficient of a. In this way, the trial function would have a higher degree of freedom next to the boundary conditions to satisfy the governing equation. Hence, using few degrees of freedom in the form of adjustable coefficients of a would significantly increase the accuracy of the trial solution with only a few adjustable coefficients.
In the next section, the present method will be employed to obtain a solution for the governing equations of Eq. (1) and (2) subject to the boundary conditions of Eqs. (3) and (4).

RESULTS AND DISCUSSION
The hybrid method and PSO algorithm were coded in MATLAB 2009 [40]. Then, the consequent results of the presented method are compared with the results of the fourth order Runge-Kutta method [24] and the VIM method [24] in literature. Besides, the governing equations (Eqs. (1) and (2) subject to the boundary conditions in Eqs. (3) and (4)) were solved using a finite difference method [41]. In order to solve the equations using the present method, the coefficients of the trial functions for H=5 (five neurons) with sigmoid units in the hidden layer and m=6 equally spaced points inside the interval [0, 1] was determined by minimizing Eq. (15) using PSO. In all of this section, the results are reported for a neural network with five neurons; otherwise, the number of the utilized neurons will be stated. The following parameters of PSO were used for this problem: Inertia weight = 0.9; Acceleration factors = 2.5; Maximum Iteration=150 and Population = 100. Table 1 displays the results of the optimization of Eqs. (15) using PSO. The optimal adjustable parameters in trial function YT1(η, a, p) and YT2(η, a, p) are summarized in Table 1 for a case with K1=0.01, Re=0.5 and Pr=1 and n=0. The optimal adjustable parameter in the first part of the trial function YT1(η, a, p), i.e., a11, was found to be a11=-0.2308 while the obtained errors are and. According to the results of Table 1, YT1(η, a, p) and YT2(η, a, p) can be written as: It should be noted that YT1 represents the solution for variable f and YT2 represents the solution for variable qn in the domain of D. The results of the present are summarized in Table 2. This table shows the results computed using YT1 and YT2 consisting of three, four and five neurons. It is obvious that the accuracy of the solution improves as the number of neurons of the neural network increases. Next, the results for f(η) and qn(η) as a function of η have been shown in Table 3. A comparison between the hybrid method with five neurons and the numerical solutions (fourth order Runge-Kutta and the finite difference method) is reported in this table. In Table 3, %Error is introduced as %Error= where NM refers to the fourth order Runge-Kutta method. Results show that the hybrid method gives an analytical solution with high accuracy by utilizing only a few neurons. A comparison between the results of the present HNNPSO and the results reported in [24] is performed in Figs. 3 and 4 and Table 3. The finite difference method used in Table 3 was implemented using bvp4c functions of MATLAB. Comparing the present method with other solutions, the results for f(η) and f '(η) are shown in Figs. (4) and (5), respectively. The results for qn(η) are also plotted in Fig. (6). Figures show that the results of the present HNNPSO are in good agreement with the numerical solution (fourth order Runge-Kutta method [24]) and the VIM method [24].   Figs. (7) and (8) show the effect of various parameters namely Reynolds number (Re), dimension parameter (n), Prandtl number (Pr) and K1 parameter on the flow and heat transfer. The effect of different Re values on the velocity profile is shown in Fig. (7), while Fig. (8) shows the effect of various parameters on the temperature profile. Again, excellent agreement is seen between the results of the proposed method and the results of the numerical solution. Figs. (7) and (8) show that increasing Reynolds number leads to an increase in the curvature of the temperature profile and decrease of qn(η) values. Also, from Fig. (8) can be concluded that for a constant value of η, temperature increases if the power-law index decreases.

CONCLUSION
The main purpose of this study was to apply Hybrid Neural Network-Particle Swarm Optimization (HNNPSO) to nonlinear differential equation arising from the similarity solution of the flow and heat transfer of a second-order viscoelastic fluid in an axisymmetric channel with a porous wall. PSO was applied to minimize the error function and find the optimized values of the adjustable parameters of the trial function. The obtained solution is a closed-form solution and represents a remarkable accuracy in comparison with the numerical methods and the VIM method, and the accuracy improves as the number of neurons of the neural network increases. Using three neurons can provide acceptable results for most of graphical and engineering purpose while using five neurons (five non-linear sigmoid terms) can provide a relative error lower than 0.01%. Therefore, the present method can be utilized to obtain highly compact and accurate close-form analytical relations for boundary layer problems, arising in various engineering problems.
Moreover, in the present study a uniform collocation point is utilized. However, a non-uniform sample points can be utilized to increase the local accuracy of the close form solution, which can be subject of future studies.

ACKNOWLEDGMENTS
The second and third authors acknowledge Shahid Chamran University of Ahvaz for the support of the present work.