Continuous Genetic Algorithm as a Novel Solver for Stokes and Nonlinear Navier Stokes Problems

Abstract

The one-dimensional continuous genetic algorithm (CGA) previously developed by the principal author is extended and enhanced to deal with two-dimensional spaces in this paper. The enhanced CGA converts the partial differential equations into algebraic equations by replacing the derivatives appearing in the differential equation with their proper finite difference formula in 2D spaces. This optimization methodology is then applied for the solution of steady-state two-dimensional Stokes and nonlinear Navier Stokes problems. The main advantage of using CGA for the solution of partial differential equations is that the algorithm can be applied to linear and nonlinear equations without any modification in its structure. A comparison between the results obtained using the 2D CGA and the known Galerkin finite element method using COMSOL is presented in this paper. The results showed that CGA has an excellent accuracy as compared to other numerical solvers.

1. Introduction

Computational fluid dynamics (CFD), which has been commercialized for few decades, is the science of predicting the effect of different phenomena on the overall performance of the processes through solving sets of mathematical equations that govern these processes using numerical techniques. There are many CFD commercial codes available today, which have been originated by the 80s and 90s, such as Fluent, CFX, Fidap, Polyflow, Phoenix, Star CD, Flow 3D, ESI/CFDRC, and SCRYU [1].

The methodology of CFD is, briefly, applying numerical method (called discretization) to develop approximations of the governing equations of fluid mechanics in the fluid region of interest. This discretization involves the governing differential equations expressed algebraically and the collection of cells (the grid). These sets of equations are solved simultaneously for the flow field variables at each node or cell, which enormously contributes.

The discretized conservation equations are solved iteratively. A number of iterations are usually required to reach a converged solution. The convergence is reached when the changes in solution variables from one iteration to the next are negligible; in other words, the residuals provide a mechanism to help monitor this trend [1]. The major advantage of CFD is the relatively low cost that is continuing to decline as computers and software become more powerful and sophisticated. Moreover, CFD simulations can be executed in a short period of time. It has also the ability of simulating real conditions.

The principle methods currently used in CFD are finite difference method (FDM), finite element method (FEM), the Galerkin method, spectral method (SM), filter scheme methods (FSM), and boundary integral equation methods (BIEM) [1]. The first four methods mentioned above are applicable in principle to all types of fluid flow problems (hyperbolic, parabolic, elliptic, and mixed) whereas (BIEM) method has more restrictions and is usually valuable for elliptic type problems [1].

Fluid flows are encountered in the daily life applications including environmental hazards (air pollution, transport of contaminants, etc.), ventilation and air conditioning applications, combustion in automobile engines, complex flows in furnaces, heat exchangers, chemical reactors, and the human body processes (blood flow, breathing, drinking, etc.) [2]. These fluid flow applications are expressed mathematically in sets of partial differential equations that represent the conservation principle of mass, momentum, and energy.

The Navier Stokes equations are the basic governing equations for viscous fluid flow. These equations are obtained by applying Newton’s law of motion to a fluid element; alternatively these equations are called the momentum equations. These equations are usually supplemented by the mass conservation equation, also called continuity equation. A combination of both momentum equation and continuity equation is usually termed as Navier Stokes equations. The viscous incompressible Navier Stokes equations (INSE) were firmly established in the 19th century as a system of nonlinear partial differential equations that describe the motion of most fluids. As a result, finding a solution with reasonable accuracy of INSE set of equations has attracted the research interest of many scientists around the world [1].

Genetic algorithm (GA) was invented by Holland in the 1960s at the University of Michigan [3]. In contrast to evolution strategies and programming, Holland’s original objective was not to design algorithms that solve specific problems but rather to formally study the phenomenon of adaptation as it occurs in nature and to develop ways in which the mechanisms of natural adaptation might be imported into computer systems [4].

GA was recently reviewed by Goldberg [5], Davis [6], and others. In general, GA is implemented using computer simulations in which an optimization problem is specified. A genetic algorithm begins with a population of typically random parents in which members of a space of candidate solutions called individuals are represented using abstract representations called chromosomes [4]. In each generation, multiple individuals are randomly selected from the current population based upon some application of fitness and then using crossover and modified through mutation to form a new population. And finally it reaches the fit population with the exact one.

There is no rigorous definition of “genetic algorithm” accepted by all in the evolutionary computation community that differentiates GA from other evolutionary computation methods. However, it can be said that most methods called “GA” have at least the following elements in common: populations of chromosomes, selection according to fitness, crossover to produce new offspring, and random mutation of new offspring [7].

Continuous genetic algorithms (CGAs) were developed by Abo-Hammour as an efficient method for finding the global solutions of smooth functions [7]. In CGAs, the smoothness of the solution curve is achieved. Their novel development has opened a wide venue for different engineering and mathematical applications. In the numerical field, CGAs have been applied for the numerical solution of boundary value problems [8–13], Laplace equations [14], differential-algebraic systems [15], fuzzy differential equations [16], and the solution of Stokes and nonlinear Navier Stokes problems [17]. CGAs have been also applied in robotics [18, 19] and optimal control problems [20, 21].

When using GA in optimization problems, there are two vital points that should be considered: firstly, seeking the interrelation between parameters of optimization and secondly whether there is some restriction on the smoothness of the resulting curve or not [8]. In the case of uncorrelated parameters or nonsmooth solution curves, the conventional genetic algorithm (GA) performs well and there is no need of the CGA. On the other hand, if the parameters are interrelated or if the smoothness of the solution curve is a must, then the continuous genetic algorithm (CGA) is preferable [18].

In summary, CGA has many advantages compared with conventional GA especially when it is applied to problems where smooth solution curve is a need and when more than one parameter is used [7].(1)CGA needs smaller memory than what GA needs, because GA uses genotype and phenotype representations of the population’s individuals while CGA utilizes only the phenotype data. This makes CGA more suitable to problems with larger number of parameters.(2)Execution time for CGA is smaller than time needed by GA because there are no encoding/decoding processes in CGA.(3)The conventional GA cannot be used in applications where the optimal solution is required.The novel application of the (CGA) possesses several advantages: (i) it guarantees the smoothness of the solution curves, (ii) the results obtained using CGA are found to be in good agreement with the analytical solutions, and (iii) it can be applied to linear and nonlinear problems without any modification in the algorithm.

In this paper, the novel method based on continuous genetic algorithms is applied to solve two problems modeled by the steady state Navier Stokes equations in two-dimensional rectangular regions [17]. These problems varied in terms of linearity, that is, linear and nonlinear problems of Dirichlet boundary conditions. The main objective of this work is to solve algebraic equations that result from finite differencing of the partial differential equations (PDE) by CGA as a global optimization technique. The efficiency of CGA method is demonstrated by comparing its results to the results of other well-known methods. Furthermore, the exact solutions for both cases are already known, which make the judgment of the CGA method results possible. The usage of exact solutions here is only to reveal the efficiency of CGA in solving problems.

2. Theory of Continuous Genetic Algorithm

Continuous genetic algorithms (CGAs) were firstly developed to avoid the sharp jumps in the parameter values which can be resulted by using conventional genetic algorithm (GA) form in some optimization problems. That can be shown when the parameters of the optimization problem are correlated with each other or if the smoothness of the solution curve is a must [8]. Continuous genetic algorithm (CGA) mainly depends on the evolution of curves in one-dimensional space, surfaces in two-dimensional space, and volumes in three-dimensional space [18].

The CGAs are of global nature; consequently, the operators of CGA result in smooth transitions but the results of the traditional GA are a step-function-like jump in the parameter values [8]. As discussed by Abo-Hammour et al. [14], the solution of two-dimensional equation in the Cartesian coordinates was found using continuous genetic algorithm (CGA). Following Laplace equation other equations such as Stokes type can be discretized in the same way [8].

The central difference discretization form of two-dimensional equations is 𝑢_𝑖,𝑗=𝑔𝑢_𝑖+1,𝑗,𝑢_{𝑖−1,𝑗},𝑢_𝑖,𝑗+1,𝑢_{𝑖,𝑗−1}. (1) The nodal residue for node (𝑖,𝑗) is defined as 𝑟_𝑖,𝑗=𝑢_𝑖,𝑗−𝑔𝑢_𝑖+1,𝑗,𝑢_{𝑖−1,𝑗},𝑢_𝑖,𝑗+1,𝑢_{𝑖,𝑗−1}. (2) The overall residue for equation is the sum of squares of the nodal residues and is stated as 𝑅=^𝑛𝑥𝑚_𝑖=1^𝑦_𝑗=1𝑟²_𝑖,𝑗, (3) where 𝑛_𝑥 is the number of unknown nodes along the 𝑥-direction and 𝑚_𝑦 is the number of unknown nodes along the 𝑦-direction.

Then, an appropriate fitness function can be given as [7] 1𝐹=.1+𝑅 (4) This means that the larger the residues value the smaller the fitness one. For the exact solution 𝑅≈0 and 𝐹≈1.

The CGAs consist of general steps as discussed by [8] and are followed during solving problems and these steps are shown in Figure 1.

Figure 1 

Block diagram for GA procedure.

(1) Initialization. This step is just like an initial guess for the solution curve that modified functions can be used to randomly generate the initial population or parents [22]. In general, the initialization functions used in the algorithm should satisfy the following conditions.(1)It should be smooth from one side and should satisfy any problem-specific constraints, if such constraints exist, from the other side [7]. That is, any smooth function that is close enough to the expected solution curve and satisfies the boundary values can be used.(2)It should result in an information-rich initial population based on the population diversity, which is split into two main types.(a)Smooth Functions Diversity. This is the diversity due to using a mixture of smooth functions in the initialization phase instead of one function. The expected solution curves are not known and correspondingly a mixture of smooth functions that satisfy the given boundary values will be beneficial to result in a diverse initial population.(b)Parameters Diversity. This is the diversity due to the random generation of the initialization parameters, and in order to make the initial population as much diverse as we can, randomness should be there to remove any bias toward any solution.It is to be noted that the closer the initialization function to the final solution, the faster the convergence speed [8–13]. However, this convergence speed improvement is minor. This means that the initialization functions have minor effect on the convergence speed of the algorithm because usually the effect of the initial population dies after few tens of generations and the convergence speed after that is governed by the selection mechanism, crossover, and mutation operators.

It is also worth mentioning that the population size is kept constant throughout the evolution process of the algorithm. This population size affects the convergence speed, the average fitness, and the corresponding errors of CGA [7–15]. Small population sizes suffer from larger number of generations required for convergence and the probability of being trapped in local minima, while large population sizes suffer from larger number of fitness evaluations that means larger execution time. Previous studies and analysis performed showed that a compromise population size for similar problems with similar number of unknown nodes is about 500 individuals [7–12].

The initial population in this case will be divided into four equal segments; each of them generated using a different set of smooth functions to obtain a diverse initial population. The first part of the initial population is generated using two-dimensional Gaussian functions, 𝑛₁(𝑥,𝑦), the second segment of the initial population is generated using a correlated two-dimensional tangent hyperbolic function, 𝑛₂(𝑥,𝑦), and the third part 𝑛₃(𝑥,𝑦) is generated using an uncorrelated two-dimensional tangent hyperbolic, while the last segment 𝑛₄(𝑥,𝑦) is generated using a mixture of the above three functions. The four functions are given as follows [22]: 𝑛₁=𝑥,𝑦_𝑘𝑖=1𝐵_1,𝑖−1×exp21−𝑙_𝑖²×𝑥−𝑐_𝑥,𝑖²𝜎²_𝑥,𝑖−2𝑙_𝑖𝑥−𝑐_𝑥,𝑖𝑦−𝑐_𝑦,𝑖𝜎_𝑥,𝑖𝜎_𝑦,𝑖+𝑦−𝑐_𝑦,𝑖²𝜎²_𝑦,𝑖,𝑛₂=𝑥,𝑦_𝑘𝑖=1𝐵_2,𝑖×𝑞×0.51+hyptan^𝑥,𝑖×𝑥−𝑐_𝑥,𝑖𝜎_𝑥,𝑖+2𝑙_𝑖𝑥−𝑐_𝑥,𝑖𝑦−𝑐_𝑦,𝑖𝜎_𝑥,𝑖𝜎_𝑦,𝑖+𝑞_𝑦,𝑖𝑦−𝑐_𝑦,𝑖𝜎_𝑦,𝑖,𝑛₃𝑥,𝑦=𝐵₃×𝑞+0.51+hyptan^𝑥𝑥−𝑐_𝑥𝜎_𝑥𝑞×hyptan^𝑦𝑦−𝑐_𝑦𝜎_𝑦,𝑛₄𝑥,𝑦=𝑛₁𝑥,𝑦+𝑛₂𝑥,𝑦+𝑛₃𝑥,𝑦, (5)(6) where 𝑘 is a random in the range [1,𝑁/5] and 𝑁 is the number of unknown nodes in the grid. 𝐵₁, 𝐵₂ are a random number in the range [−2,2] and 𝐵₃ is a random number in the range [−0.25,0.25]. 𝑐_𝑥 is a random number in the range [𝑥_min,𝑥_max]. 𝜎_𝑥 is a random number in the range [0,𝑥_max−𝑥_min]. 𝑐_𝑦 is a random number in the range [𝑦_min,𝑦_max]. 𝜎_𝑦 is a random number in the range [0,𝑦_max−𝑦_min]. 𝑙 is a random number in the range [−0.5,0.5]. 𝑞_𝑥, 𝑞_𝑦 are random numbers in the range [0,2] and in 𝑛₄(𝑥,𝑦), 𝑘=1, for both 𝑛₁ and 𝑛₂.

The upper functions 𝑛(𝑥,𝑦) are then rescaled within the range [𝐼_min,𝐼_max] representing the expected minimum and maximum initial nodal values.

It is to be noted that the parameter “𝑙” appearing in 𝑛₁(𝑥,𝑦) and 𝑛₂(𝑥,𝑦) as given in (5) is chosen between −0.5 and +0.5 based on the second term of the argument of the exponential or the tangent hyperbolic function which is related to the dual effect of both coordinates (𝑥,𝑦). The coefficient of that term is “2𝑙” and in order to have a contribution for that term within the range [−1,1], then the range of the “𝑙” values should be within the range [−0.5,0.5].

(2) Fitness Evaluation. In this step, the measure of goodness of each individual in the population is calculated. Fitness will be applied to all populations until reaching the maximum criteria which is around one. This step is applied according to (2)–(4) [7].

(3) Selection. Selection process is to choose the individuals according to their relative fitness in order to enter the devotion process to create the new individuals for the next generation. This step ensures that the overall quality of the population increases from one generation to the next [7].

(4) Crossover. Crossover is a major operation that really empowers the GA. It operates on two randomly selected highly fitted individuals at a time and generates offspring by combining both parent individuals’ features [14].

In two-dimensional cases, crossover function will also be expanded into the two-dimensional following form [14]:𝑂₁𝑥,𝑦=𝑆𝑥,𝑦×𝐼₁𝑥,𝑦+1−𝑆𝑥,𝑦×𝐼₂𝑥,𝑦,𝑂₂𝑥,𝑦=1−𝑆𝑥,𝑦×𝐼₁𝑥,𝑦+𝑆𝑥,𝑦×𝐼₂𝑥,𝑦, (7a)(7b)where 𝑂₁, 𝑂₂ are the two children resulted from the crossover of the two parents 𝐼₁ and 𝐼₂. 𝑆 is the mixing function in the range [0,1].

Four types of mixing function can be used: the first three types are tangent hyperbolic-related functions, while the last one is of Gaussian shape [14]: 𝑆₁𝑞𝑥,𝑦=0.5×1+hyptan^𝑥𝑥−𝑐_𝑥𝜎_𝑥+𝑞_𝑦𝑦−𝑐_𝑦𝜎_𝑦,𝑆₂𝑞𝑥,𝑦=0.5×1+hyptan^𝑥𝑥−𝑐_𝑥𝜎_𝑥+𝑞_𝑥,𝑦𝑥−𝑐_𝑥𝑦−𝑐_𝑦𝜎_𝑥𝜎_𝑦+𝑞_𝑦𝑦−𝑐_𝑦𝜎_𝑦,𝑆₃𝑞𝑥,𝑦=0.5×1+hyptan^𝑥𝑥−𝑐_𝑥𝜎_𝑥𝑞×hyptan^𝑦𝑦−𝑐_𝑦𝜎_𝑦,𝑆₄=𝑥,𝑦_{𝑁𝑐𝑖=1}−1exp21−𝑙²_𝑖×𝑥−𝑐_𝑥,𝑖²𝜎²_𝑥,𝑖−2𝑙_𝑖𝑥−𝑐_𝑥,𝑖𝑦−𝑐_𝑦,𝑖𝜎_𝑥,𝑖𝜎_𝑦,𝑖+𝑦−𝑐_𝑦,𝑖²𝜎²_𝑦,𝑖, (8) where 𝑐_𝑥=𝑥_min+0.1×𝑥_max−𝑥_min,0.9×𝑥_max−𝑥_min,𝑐_𝑦=𝑦_min+0.1×𝑦_max−𝑦_min,0.9×𝑦_max−𝑦_min. (9)𝜎_𝑥, 𝜎_𝑦 are chosen such that the tangent hyperbolic function achieves its complete transition from −1 to 1 within the given ranges of both 𝑥 and 𝑦. 𝑞_𝑥, 𝑞_𝑦 have values +1, −1 randomly. 𝑞_𝑥,𝑦 is random number in the range [−0.2,0.2]. 𝑁_𝑐 represents the number of crossing points and it is chosen in the range [1,𝑁/20]. 𝑙 is a random number within the range [−0.1,0.1].

(5) Mutation. Mutation introduces random variations into the population. It is applied to a single chromosome only. It is usually performed with low probability. Genetic diversity of the population must be maintained so in order to guard against premature convergence, mutation is usually taken into consideration [14].

Two-dimensional mutation processes can be governed by the following formula [14]: 𝑂^𝑚𝑥,𝑦=𝑂𝑥,𝑦+𝑤×𝑀𝑥,𝑦, (10) where 𝑂(𝑥,𝑦) represents the child produced from the crossover process. 𝑂^𝑚(𝑥,𝑦) is the mutated child. 𝑀 is the mutation function. 𝑤 is a random number in the range [−𝑅_ave,𝑅_ave] and 𝑅_ave is average overall residue of the population. There are three types of mutation functions [7, 8]: 𝑀₁𝑥,𝑦=𝑚_dc1−1+exp21−𝑙_𝑖×𝑥−𝑐_𝑥²𝜎²_𝑥−2𝑙𝑥−𝑐_𝑥𝑦−𝑐_𝑦𝜎_𝑥𝜎_𝑦+𝑦−𝑐_𝑦²𝜎²_𝑦,𝑀₂𝑥,𝑦=𝑚_dc2×𝑞+0.51+hyptan^𝑥𝑥−𝑐_𝑥𝜎_𝑥+𝑞_𝑥,𝑦𝑥−𝑐_𝑥𝑦−𝑐_𝑦𝜎_𝑥𝜎_𝑦+𝑞_𝑦𝑦−𝑐_𝑦𝜎_𝑦,𝑀₃𝑥,𝑦=𝑚_dc3, (11) where 𝑐_𝑥 is a random number in the range [𝑥_min,𝑥_max]. 𝜎_𝑥 is a random number in the range [0,𝑥_max−𝑥_min]. 𝑐_𝑦 is a random number in the range [𝑦_min,𝑦_max]. 𝜎_𝑦 is a random number in the range [0,𝑦_max−𝑦_min]. 𝑞_𝑥, 𝑞_𝑦 have values +1, −1 randomly. 𝑞_𝑥,𝑦 is random number in the range [−0.2,0.2]. 𝑚_dc1, 𝑚_dc2 are some random dc offset values within the range [−0.25,0.75]. 𝑙 is a random number within the range [0,1]. 𝑚_dc3 is a random number within the range [−0.5,0.5].

There are three methods for selecting the value of 𝑐 that is usually used in the mutation process [8].(1)Random method: in this method 𝑐 is randomly selected as given in the crossover step.(2)Lamarckian method: this method depends on the value of residuals at certain mesh point that is the larger residual value at that point then the probability of choosing this point as a mutation center.(3)Deterministic manner: in this method the mesh point with the maximum residual is chosen as the mutation center.CGA also used “extinction and immigration” operator. This operator is applied when all individuals in the population are identical or when the improvement in the fitness value of the best individual over a certain number of generations is less than some threshold value. This means that no new information will be obtained through crossover process [7]. The CGA thus tends to stagnate; “extinction and immigration” operator is used to bypass this difficulty. This operator, as indicated by its name, consists of two stages; the first stage is the extinction process where all of the individuals in the current generation are removed except the best of generation individual [7]. The second stage is the mass-immigration process, where the extinct population is filled out again by generating 𝑁_𝑝−1 individuals to keep the population size fixed. The generated population is divided into two equal segments each of 𝑁_𝑝/2 size; the first segment, with 2≤𝑗≤𝑁_𝑝/2, is generated as in the initialization phase, while the other segment is generated by performing continuous mutations to the best-of-generation individual as given by the following formula [14]: 𝑂_𝑖𝑥,𝑦=𝑂₁𝑁𝑥,𝑦+𝐹×𝐺𝑥,𝑦,2≤𝑖≤^𝑝2, (12) where 𝑂_𝑖 is the 𝑖th parent generated using immigration operator. 𝑂₁ is the best of generation individual. 𝐹 is a random number in the range [−𝑅₁,𝑅₁], where 𝑅₁ is the overall residue of the best individual and 𝐺(𝑥,𝑦) is one of the following functions [14]: 𝐺₁𝑥,𝑦=𝑚_dc×𝑞+0.51+hyptan^𝑥𝑥−𝑐_𝑥𝜎_𝑥+𝑞_𝑥,𝑦𝑥−𝑐_𝑥𝑦−𝑐_𝑦𝜎_𝑥𝜎_𝑦+𝑞_𝑦𝑦−𝑐_𝑦𝜎_𝑦,𝐺₂𝑥,𝑦=𝑚_dc×𝑞+0.51+hyptan^𝑥𝑥−𝑐_𝑥𝜎_𝑥𝑞×hyptan^𝑦𝑦−𝑐_𝑦𝜎_𝑦,𝐺₃𝑥,𝑦=𝑚_dc−1+exp21−𝑙_𝑖×𝑥−𝑐_𝑥²𝜎²_𝑥−2𝑙𝑥−𝑐_𝑥𝑦−𝑐_𝑦𝜎_𝑥𝜎_𝑦+𝑦−𝑐_𝑦²𝜎²_𝑦𝐺₄𝑥,𝑦=𝑚_dc4. (13)𝑞_𝑥, 𝑞_𝑦 have values +1, −1 randomly. 𝑞_𝑥,𝑦 is random number in the range [−1,1]. 𝑚_dc are some random dc offset values within the range [−0.25,0.75]. 𝑙 is a random number within the range [−0.5,0.5]. 𝑚_dc4 is a random number within the range [0,1].

(6) Replacement. Applying the genetic operators to the parent’s population generates the offspring’s population and then the parents’ population is totally replaced by the offspring’s population. This is called nonoverlapping, generational, or replacement. In this step the “life cycle” of the population can be completed [7].

(7) Termination. When some convergence criterion is met, the GA is terminated. The convergence criteria are as follows: (i) the fitness of the best individual so far found exceeds a threshold value, (ii) the maximum number of generations is reached or the progress limit, and (iii) the improvement in the fitness value of the best member of the population over a specified number of generations is less than some predefined threshold. The best solution of the problem is found after the termination of the algorithm [7].

3. Numerical Results for Stokes Equations

The general forms of Stokes equations are shown in the following: −𝜕²𝑢𝜕𝑥²+𝜕²𝑢𝜕𝑦²=𝑓_𝑥,−𝜕²𝑣𝜕𝑥²+𝜕²𝑣𝜕𝑦²=𝑓_𝑦, (14) in which the right hand side is the shear stress term and the left hand side is the external force applied on the fluid, with continuity equation 𝜕𝑢+𝜕𝑥𝜕𝑣𝜕𝑦=0. (15)Dennis and Hudson [23] tried to solve the nonlinear form of Navier Stokes equations but with no pressure term as shown in the following equations: 𝑢𝜕𝑢𝜕𝑥+𝑣𝜕𝑢=𝜕𝜕𝑦²𝑢𝜕𝑥²+𝜕²𝑢𝜕𝑦²+𝑓_𝑥,𝑢𝜕𝑣𝜕𝑥+𝑣𝜕𝑣=𝜕𝜕𝑦²𝑣𝜕𝑥²+𝜕²𝑣𝜕𝑦²+𝑓_𝑦. (16)(17) Continuity equation is as follows: 𝜕𝑢+𝜕𝑥𝜕𝑣𝜕𝑦=0, (18) within the range of 𝑥=[0,1] and 𝑦=[0,1].

Wang et al. [24] tried to solve Stokes equation using finite element method within the region 𝑥=[0,1], 𝑦=[0,1]; the authors used the following boundary conditions: 𝑢()𝑢()𝑣()𝑣()0,𝑦=0,𝑢𝑥,0=0,1,𝑦=0,𝑢𝑥,1=0,0,𝑦=0,𝑣𝑥,0=0,1,𝑦=0,𝑣𝑥,1=0. (19) The exact solutions as given by Wang et al. [24] are as follows: ()𝑥(),()𝑦).𝑢=−2𝑥𝑦𝑥−1𝑦−1𝑥−12𝑦−1𝑣=2𝑥𝑦𝑥−1𝑦−1𝑦−1(2𝑥−1 (20) The right hand sides for (14) were given by Wang et al. [24] in the following forms: 𝑓_𝑥=8𝑥²(𝑥−1)²𝑦−1+4𝑥²2𝑦−1(𝑥−1)²+8𝑥²𝑦(𝑥−1)²+4𝑦2𝑦−1(𝑥−1)²𝑦−1+4𝑥²𝑦(),𝑓2𝑦−1𝑦−1+8𝑥𝑦2𝑥−22𝑦−1𝑦−1_𝑦=−8𝑦²(𝑥−1)𝑦−1²−4𝑦²(2𝑥−1)𝑦−1²−8𝑥𝑦²𝑦−1²−4𝑥(2𝑥−1)(𝑥−1)𝑦−1²−4𝑥𝑦²()()().2𝑥−1)(𝑥−1−8𝑥𝑦2𝑥−12𝑦−2𝑥−1 (21) This right hand sides are resulted by substituting the exact solutions shown in (20) into (14). The nodal values of the exact solutions are shown in Tables 1 and 2.

Table 1 

Actual 𝑢 nodal values for Stokes equations including boundary values.

Table 2 

Actual 𝑣 nodal values for Stokes equation including boundary values.

Moreover, the graphical representations of those actual values are shown in Figures 3(a) and 3(b).

The boundary values are known from the main problems but the other nodal values are unknown values of 𝑢 and must be solved. These nodes can be represented in the (10×10) mesh of 100 elements, as shown in Figure 2.

Figure 2 

Rectangular mesh with step size equal 0.1.

(a)

(b)

(c)

(d)

(e)

(f)

Figure 3 

For Stokes case. (a) Graphical representation for actual values of 𝑢. (b) Graphical representation for actual values of 𝑣. (c) Finite element solution for 𝑢 variable using COMSOL. (d) Finite element solution for 𝑣 variable using COMSOL. (e) Graphical representation for nodal 𝑢 values using CGA. (f) Graphical representation for nodal 𝑣 values using CGA.

With the same previous conditions COMSOL solution was obtained and the nodal values of 𝑢 and 𝑣 variables are shown in Tables 3 and 4. Moreover, the graphical representations for 𝑢 and 𝑣 COMSOL solutions are shown in Figures 3(c) and 3(d).

Table 3 

Nodal values of 𝑢 variable using COMSOL for Stokes case solution.

Table 4 

Nodal values of 𝑣 variable using COMSOL for Stokes case solution.

Appling (22) and using the errors between the exact and COMSOL values, it is found that 𝐿₂ norm error equals to 1.29×10⁻⁵ for 𝑢 variable and about 1.21×10⁻⁵ for 𝑣 variable with average time for solving the problem of about 0.045: 𝐿₂normerror=^𝑛_𝑖=1𝑢_exact−𝑢_numerical². (22)

3.1. Solution of Stokes Equations Using Continuous Genetic Algorithm

All steps of CGA are applied here and they are repeated till the convergence criterion is met (fitness = 0.99), where the best unknown values of nodes are reached. The central difference discretization form of Stokes equation was shown in (23) but for both 𝑢 and 𝑣 separately; also (24) and (25) are used in order to define the unknown residual nodal 𝑢 and 𝑣 values: 𝜕²𝑢𝜕𝑥²+𝜕²𝑢𝜕𝑦²+𝑓_𝑥=𝑢_{𝑖−1,𝑗}+𝑢_𝑖+1,𝑗+𝑢_{𝑖,𝑗−1}+𝑢_𝑖,𝑗+1−4𝑢_𝑖,𝑗ℎ²+𝑓_𝑥, (23) where ℎ is the step size in which ℎ_𝑥 is the step size in the 𝑥-direction and ℎ_𝑦 is the step size in the 𝑦-direction. And in this case ℎ_𝑥=0.1 and ℎ_𝑦=0.1, again the problem is defined over the rectangular region 𝑥=[0,1], 𝑦=[0,1].

Then, after rearrangement, one has 𝑢_𝑖,𝑗=𝑢_{𝑖−1,𝑗}+𝑢_𝑖+1,𝑗+𝑢_{𝑖,𝑗−1}+𝑢_𝑖,𝑗+1+𝑓_𝑥4×ℎ². (24) The nodal residue for node (𝑖,𝑗) is defined as 𝑟_𝑖,𝑗=𝑢_𝑖,𝑗−𝑢_{𝑖−1,𝑗}+𝑢_𝑖+1,𝑗+𝑢_{𝑖,𝑗−1}+𝑢_𝑖,𝑗+1+𝑓_𝑥4×ℎ². (25) Then, the overall residues were solved for (25) as the sum of squares of the nodal residues, as stated in (3); after that the appropriate fitness function resulted from (4).

The obtained solution for Stokes equation [Dirichlet boundary conditions] using continuous genetic algorithms (CGAs) are shown in Tables 5 and 6, including the nodal values of the boundaries which are equal to zero. And the graphical representations of these values are shown in Figures 3(e) and 3(f).

Table 5 

Nodal 𝑢 values including boundaries for Stokes equation using CGA.

Table 6 

Nodal 𝑣 values including boundaries for Stokes equation using CGA.

By calculating the errors between the actual nodal values and CGA one, the 𝐿₂ norm error is found to be equal 3.02×10⁻⁴ for 𝑢 variable and 1.08×10⁻⁴ for 𝑣 variable. The time needed to reach solution in CGA was about 7523 seconds.

Wang et al. [24] used finite element method with a mesh size of 8×8 for solving Stokes equations and then reported 5.05×10⁻⁴ as 𝐿₂ norm error of the solution.

Comparing the error values which resulted by applying CGA with the errors given by finite element method by Wang, it is found that the error for solution obtained by CGA is close to that error for solution obtained by finite element. On the other hand, 𝐿₂ norm error for using COMSOL is less, that was about 1.29×10⁻⁵ for 𝑢 variable and about 1.21×10⁻⁵ for 𝑣 values and even these errors are less than the errors of CGAs method.

The average number of generations was 1139, and the average residual value was about 6.339×10⁻⁴. Then, the average fitness according to (4) is very close to the maximum fitness value of 1.

The convergence state of nodal values of the nodes was presented along all generations such as an example for [𝑢(2,2), 𝑢(6,6), and 𝑢(6,10)] nodes excluding the boundary values as shown in Figures 4, 5, and 6, respectively. It is concluded from the figures that steady state is reached for 𝑢(2,2) nodes much earlier than that for 𝑢(6,6) and 𝑢(6,10). The reason for that is that the node (2,2) location is closer to the boundaries compared with the other nodes locations. For the same reason, the convergence for node (6,10) is faster than for (6,6).

Figure 4 

Nodal values of the point 𝑢(2,2) along all generations for Stokes and nonlinear Stokes cases.

Figure 5 

Nodal values of the point 𝑢(6,6) along all generations for Stokes and nonlinear Stokes cases.

Figure 6 

Nodal values of the point 𝑢(6,10) along all generations for Stokes and nonlinear Stokes cases.

The nodal residues in this case were less than 14.925×10⁻⁴ for 𝑢 values as shown in Table 7 and were less than 18.441×10⁻⁴ for 𝑣 values as shown in Table 8. It can be concluded from these values that the solution is close to the exact solution.

Table 7 

Residual 𝑢 values for all nodes for Stokes equation after using CGA (×10⁻⁴).

Table 8 

Residual 𝑣 values for all nodes for Stokes equation after using CGA (×10⁻⁴).

Maximum fitness was reached with smooth increasing toward the maximum value along all generations and this assures how much the solution of CGA is close to the true solution. This behavior is shown in Figure 7(a).

(a)

(b)

Figure 7 

Maximum fitness values along all generations for Stokes and nonlinear Stokes cases.

4. Numerical Results for Nonlinear Navier Stokes Equations

The following conditions were taken into consideration while solving the problem. Boundary conditions are as follows: 𝑢()𝑢()𝑣()𝑣0,𝑦=0,𝑢𝑥,0=0,1,𝑦=𝑦,𝑢𝑥,1=𝑥,0,𝑦=0,𝑣𝑥,0=0,1,𝑦=𝑦²,𝑣(𝑥,1)=𝑥². (26) The exact solutions are as follows: 𝑢=𝑥𝑦,𝑣=𝑥²𝑦². (27)The second terms in the right hand side of (16) and (17) are defined in (28) and they were obtained by substituting the exact solutions given in (27) into (16) and (17): 𝑓_𝑥=𝑥³𝑦²+𝑥𝑦²,𝑓_𝑦=2𝑥⁴𝑦³+2𝑥²𝑦³−2𝑥²−2𝑦². (28) The actual nodal values of both 𝑢 and 𝑣 variables are shown in Tables 9 and 10. The graphical representations of those actual values are shown in Figures 8(a) and 8(b).

Table 9 

Actual 𝑢 nodal values for nonlinear Navier Stokes equation including boundary values.

Table 10 

Actual 𝑣 nodal values for nonlinear Navier Stokes equation including boundary values.

(a)

(b)

(c)

(d)

(e)

(f)

Figure 8 

For nonlinear Navier Stokes. (a) Graphical representation for actual values of 𝑢. (b) Graphical representation for actual values of 𝑣. (c) Finite element solution for 𝑢 using COMSOL. (d) Finite element solution for 𝑣 using COMSOL. (e) CGA solution for 𝑢. (f) CGA solution for 𝑣.

Tables 11 and 12 show the obtained nodal values of 𝑢 and 𝑣 variables as a result of applying COMSOL. Moreover, Figures 8(c) and 8(d) show the graphical representations of these solution values.

Table 11 

Nodal values of 𝑢 using finite element method in COMSOL for the nonlinear case.

Table 12 

Nodal values of 𝑣 using finite element method in COMSOL for nonlinear case.

For solution obtained by COMSOL, the norm error was found to be about 1.6258 for 𝑢 variable values and about 0.5845 for 𝑣 variable. And the time needed to solve this problem was about 0.546 second. For the nonlinear Navier Stokes, the results of using COMSOL to solve 𝑢 variable for different sizes of meshes are shown in Table 13.

Table 13 

COMSOL results for different meshes.

4.1. Solution of Nonlinear Navier Stokes Using Advanced Continuous Genetic Algorithm

Firstly, the discretization of (16) is 𝑢_𝑖,𝑗𝑢_𝑖+1,𝑗−𝑢_{𝑖−1,𝑗}2ℎ_𝑥+𝑣_𝑖,𝑗𝑢_𝑖,𝑗+1−𝑢_{𝑖,𝑗−1}2ℎ_𝑦=𝑢_{𝑖−1,𝑗}−2𝑢_𝑖,𝑗+𝑢_𝑖+1,𝑗ℎ²_𝑥+𝑢_{𝑖,𝑗−1}−2𝑢_𝑖,𝑗+𝑢_{𝑖,𝑗−1}ℎ²_𝑦+𝑓_𝑥 (29) and is the same for the second equation. Then, for residues, one has 𝑟_1𝑖,𝑗=𝑢_𝑖,𝑗𝑢_𝑖+1,𝑗−𝑢_{𝑖−1,𝑗}2ℎ_𝑥+𝑣_𝑖,𝑗𝑢_𝑖,𝑗+1−𝑢_{𝑖,𝑗−1}2ℎ_𝑦−𝑢_{𝑖−1,𝑗}−2𝑢_𝑖,𝑗+𝑢_𝑖+1,𝑗ℎ²_𝑥+𝑢_{𝑖,𝑗−1}−2𝑢_𝑖,𝑗+𝑢_{𝑖,𝑗−1}ℎ²_𝑦−𝑓_𝑥, (30) Overall residues are the summation for both first and second equations residues: 𝑅1=^𝑛𝑥𝑚_𝑖=1^𝑦_𝑗=1𝑟²_1𝑖,𝑗,𝑅2=^𝑛𝑥𝑚_𝑖=1^𝑦_𝑗=1𝑟²_2𝑖,𝑗,𝑅=𝑅1+𝑅2. (31)

A solution using continuous genetic algorithm (CGA) for case 2 (nonlinear Navier Stokes equation [Dirichlet boundary conditions]) is shown in Tables 14 and 15, including the nodal values of the boundaries. And the graphical representation of these values is shown in Figures 8(e) and 8(f).

Table 14 

Nodal 𝑢 values including boundaries for nonlinear Navier Stokes equation using CGA.

Table 15 

Nodal 𝑣 values including boundaries for nonlinear Navier Stokes equation using CGA.

The results of applying CGA can be summarized in calculating the 𝐿₂ norm error that was found to be about 1.42×10⁻⁴ for 𝑢 variable values and about 9.76×10⁻⁵ for 𝑣 variable values. The time that CGA needed to solve such kind of problems was about 7456 seconds.

It can be concluded from the results that solving nonlinear equations using finite element method in COMSOL application have a low accuracy. On the other hand, the trend of accuracy for CGA method is the same for both linear and nonlinear cases.

The CGAs process for solving nonlinear Navier Stokes equations in visual basic tools about 7456 seconds, in which the average number of generations was 1127 generations. Moreover, the average residual value was about 5.75433×10⁻⁴. Then, the average fitness values according to the equation 𝐹=1/(1+𝑅)=0.9994 are very close to the maximum fitness value of unity.

The convergence state of nodal values of the nodes was presented along all generations such as an example for [𝑢(2,2), 𝑢(6,6), and 𝑢(6,10)] nodes excluding the boundaries values as shown in Figures 4, 5, and 6, respectively. It is concluded from the figures that steady state is reached for 𝑢(2,2) nodes much earlier than that for 𝑢(6,6) and 𝑢(6,10). The reason for that is that the node (2,2) location is much closer to the boundaries compared with the other nodes locations. For the same reason, the convergence for node (6,10) is faster than for (6,6).

The nodal residues in this case were less than 22.1059×10⁻⁴ for 𝑢 values as shown in Table 16 and were less than 25.686×10⁻⁴ for 𝑣 values as shown in Table 17. Maximum fitness was reached with smooth increasing toward the maximum value along all generations and this assures how much the solution of CGA is close to the true solution. This behavior is shown in Figure 7(b).

Table 16 

Residual 𝑢 values for all nodes for nonlinear Navier Stokes equation after using CGA (×10⁻⁴).

Table 17 

Residual 𝑣 values for all nodes for nonlinear Navier Stokes equation after using CGA (×10⁻⁴).

Table 18 summarizes the overall performance of the used solvers solutions for the selected cases of Navier Stokes equations. As shown in the table, the CGA method has good results in solving both linear and nonlinear problems. It is worth mentioning that in the Stokes case the accuracy of COMSOL was higher than that of CGA and Wang et al. results. However, in the nonlinear case, COMSOL results’ errors were higher than those of the CGA results.

Table 18 

Summary of all Navier Stokes cases results.

As a result, it can be concluded that COMSOL software has good predictions and small errors for the linear cases but not for nonlinear problems since its accuracy dramatically deceased. This makes CGA technique advantageous over the COMSOL technique because of its adaptability to the nonlinear nature of the problems with reasonable and consistent errors. However, it should be mentioned that the CGA method requires longer time periods than COMSOL in solving the problems. This is a new challenge for CGA, which could be overcome through the continuous development of the CGA algorithm or the computer hardware. However, it is worth mentioning that CGA has good accuracy even in the first run of solving such kind of problems.

5. Conclusions

In this paper, two-dimensional CGA was applied for the solution of steady-state two-dimensional Stokes and nonlinear Navier Stokes problems. The proposed CGA for the solution of partial differential equations has the following advantages: first, it does not require any modification while switching from the linear to the nonlinear case, only the discretization of the problem will be changed, while the solution procedure using CGA will be the same; second, the algorithm requires minimal knowledge about the aimed problem that is the discretized form of the partial differential equation in the residual form; third, the method is found to be simple, efficient, and attractive with great potential in engineering applications. The resulting solutions provided good accuracy for Stokes case and excellent accuracy as compared to other commercial packages for the nonlinear case.

Conflict of Interests

The authors declare that there is no conflict of interests regarding the publication of this paper.