© R. El Maani et al., published by EDP Sciences, 2020

1 Introduction

The current need for faster and more accurate methods for calculating flow fields around technical interest configurations directed the rapid progress of CFD. In the past, CFD was the method of choice in the design of many industrial components and processes in which fluid or gas flows play a major role with many commercial packages available for modeling flow in or around objects. Features and details that are difficult, expensive or impossible to measure or visualize experimentally are shown due to the computer simulations, such as transition from laminar to turbulent flow that plays an important role in determining the flow features and in quantifying the airfoil performance. The creation of the geometry and the meshes with a preprocessor constitute the first step in modeling a problem and a successfully generated mesh for the domain geometry takes usually the majority of time spent on a CFD project in the industry allowing a compromise between desired accuracy and solution cost. Then the governing equations of the problem are solved and the solution are computed and monitored using the solver and the physical models. Afterwards, the results are examined and saved and if necessary, revisions to the numerical or physical model parameters are considered.

In the preliminary design of any airplane, the geometrical dimensions, with selection of the principal components, are analyzed with the aerodynamics characteristics to find the suitable combination. The provided design of the airfoils shapes aims to have, for given flight conditions, high lift values at low drag. More parameters, such taper ratio or aspect ratio, affect the overall lift (L) and drag (D) of the wing. The Reynolds number is also very important for airfoil performance. This number bounds the maximum thickness-to-chord ratio, beyond which point the airfoil will have unacceptable performance. It also determines the achievable section maximum lift coefficient and lift-to-drag ratio. Airfoils provide two-dimensional lift, drag and pitch momentum, which is equivalent to the characteristics of a section of an infinite span wing. Real wings, the wing with finite span, behave quite differently. Optimizing these ratios and forces presents a very usefull improvement in the aircraft wings aerodynamic performances.

By definition, optimization aims to obtain the best possible performance of a model, which is formulated in mathematically as the minimization of a function or a set of functions at the same time. This refers to single-objective optimization and multiobjective optimization, respectively. A multiobjective optimization problem (MOP) involves more than one objective function. The task of multiobjective optimization is not to find an optimal solution corresponding to each objective function but to find a set of solutions called Pareto-optimal front [1–3]. The multiplicity of principles involved in real world applications, depending on several variables of a given model, makes the use of native MOP techniques crucial and handling of these participating variables appropriately is the key for successful optimization.

A satisfactory solution under highly nonlinear and complex constraints is conditionned by optimizing simultaneously several objectives. The multiobjective nature of the problem is given by the existence of contradictory objectives meaning that the improvement of one objective implies the deterioration of another [4]. Multiobjective optimization then consists in seeking all the solutions which correspond to the best compromises between the objectives of security to be maximized and of cost to be minimized. The design problem we face therefore involves several contradictory objectives, which requires the use of an optimization method capable of managing the multiobjective nature of the problem. We are more particularly interested in the multiobjective optimization algorithm by Evolutionary Algorithm techniques (EA), which are generally used to find the global optimum when the objective function of an optimization problem is indistinguishable and not linear [5,6]. EAs have been used for problems of dynamic analysis of chemical processes [7], phase stability and equilibrium calculations for reactive systems [8], communication applications [9], dynamic optimization of biochemical processes [10], mechanical design and many other engineering problems.

It is very important, in an EA, to protect a population's diversity for the population's ability to iteratively support its development. An EA tries to develop gradually, through a “trial individual”, an individual with a better fitness value using a combination of a choosen existing individuals due to various genetic operators. The trial individual with a better fitness value replaces the original one in the next-generation population. EAs radically differ from one another based on their strategies for generating trial individuals. Because these strategies have a considerable effect on their problem-solving success and speed. It is believed that recombination, crossover, mutation, selection and adaptation [11–13], constitute the genetic diversity of a population. Many EAs are based on basic genetic rules, such as the covariance matrix adaptation evolution strategy (CMAES) [14], genetic algorithms [11,15] and the differential evolution algorithm (DE) [16–18].

NSGA-II [11], SPEA-2 [19] and NNIA [20] are examples of multiobjective optimization algorithms that can process a set of solutions simultaneously in a single run and for which an EA presents an effective tool. Recently, a new evolution algorithm, called backtracking search algorithm (BSA), was proposed and applied in several numerical single-objective optimization problems by Civicioglu [21]. It has been used in the synthesis of concentric circular antenna arrays by Guney et al. [22] and the parameter identification of hyperchaotic system [23]. In this algorithm, information obtained from past generations is used to search for better fitness solutions. It has good convergence and it is very effective in solving numerical optimization problems, compared with the performances of six widely used EA algorithms: JDE, ABC, CMAES, PSO, SADE and CLPSO [21], due to its unique mechanism for generating a trial individual. To generate trial individuals, BSA uses selection, mutation and crossover. In contrast with many genetic algorithms such as DE and its derivatives, for each target individual from individuals of a randomly chosen previous generation, BSA uses only one direction individual based on a random mutation strategy. BSA uses a non-uniform crossover strategy that is more complex than the crossover strategies used in many genetic algorithms [24].

Our objective, in this paper, is to use the backtracking search algorithm for solving multiobjective optimization problems (BSAMO). BSA is extended to deal with multiobjective design problems using the fast non-dominated sorting procedure and the crowding distance. The application of the proposed algorithm to real aerodynamic problems can be considered as a multi-disciplinary design optimization (MDO). The MDO is a field of engineering that focuses on the use of numerical optimization for the design of systems that involve a number of disciplines or subsystems [25–27].

This paper is organized as follows. In Section 3 we introduce the mesh morphing tool used to achieve the motion of the mesh at any specific control points and Section 4 presents some basic concepts of the multiobjective optimization and the Backtracking Search Algorithm for Multiobjective Optimization (BSAMO). Section 5 presents the numerical simulation consisting first of computing the transonic, compressible flow over a NACA0012 airfoil compared to published experimental data, then the coupled solution procedure for the CFD optimization process was presented with the multiobjective optimization results. Finally, the main conclusions are drawn in Section 6.

2 Governing equations

2.1 Fluid mechanics

We are generally interested in the properties of the flow at certain locations in the fluid domain. Therefore, to describe the fluid flows, the Eulerian formulation is generally used. We limit ourselves to the case of Newtonian fluids which are linear viscous isotropic fluids and the most important for practical applications. The Cauchy T stress tensor of these fluids is characterized by the following material law:(1)where p is the pressure, μ the dynamic viscosity, v_i the velocity vector with respect to Cartesian coordinate x_i, and the Kronecker symbol δ_ij. On the spatial fluid mechanics domain, the Navier-Stokes equations of incompressible flows may be written as:(2)where f, v, and ρ are the external force, velocity, and the density, respectively, and σ the stress tensor is defined as:(3)

Here p is the pressure, I is the identity tensor, μ is the dynamic viscosity, and ε(v) is the strain-rate tensor given by:(4)

The boundary conditions associated to the fluid domain are:(5)where can be a known velocity profile at the boundary.

2.2 Turbulence model

Turbulence models is a modeling of the terms related to fluctuations in statistical unstable Navier-Stokes equations by introducing mean and fluctuating quantities which form the Reynolds Averaged Navier-Stokes Equations (RANS), due to the average statistical procedure used to obtain these equations the turbulence models based on the RANS equations are called statistical turbulence models. In our study, the k − ω SST model will be used [28].

2.2.1 K-Omega SST model

According to Menter [29], the model k − ω of transport by shear stress (SST) combines the two models k − ω and k − ε, two equations are solved, one for the specific dissipation ω and the other for the kinetic energy of turbulence k, with a robust and precise formulation in the region close to the wall with the freestream independence in the far field.

The SST k − ω model has a similar form to the standard k − ω model:(6) (7)where D_ω represents the cross-diffusion term. G_ω is given by:(8)

represents the production of turbulence kinetic energy, and is defined as:(9)

The SST k − ω model is based on both the standard k − ω model and the standard k − ε model. To blend these two models together, the standard k − ε model has been transformed into equations based on k and ω, which leads to the introduction of a cross-diffusion term D_ω is defined as:(10)

Model constants are:(11)

3 Mesh morpher

For typical engineering problems, the shape sensitivity field can have smoothness properties that are not adequate to define a shape modification. Mesh morphing technology is used here for two-dimensional system in a two-fold role. The first role is as a smoother for the surface sensitivity field. The second role is to provide smooth distortions not only of the boundary mesh, but also the interior mesh. This approach is very appealing since it functions for arbitrary mesh cell types [30,31].

A Cartesian or cylindrical region is defined such that it encompasses all or a part of the problem domain. Only the mesh nodes that fall within this region may move as a result of the morphing operation. A local (u, v) coordinate system is defined, where u ∈ [0, 1] and v ∈ [0, 1]. A regular array of N_u × N_v control points is then distributed in the control volume. The motion of the mesh at any point in the domain is determined by the movement of the control points.

In general, an optimization problem may include both free-form minimization or maximization of the observable(s) and constraints that are localized in space. This requires an approach that can provide a deformation field that is well-behaved and consistent with manufacturability requirements, while still allowing locally sharper deformations where required to satisfy the imposed constraints.

To achieve this, two coincident sets of control points are created. Bernstein polynomials and B-splines are used, respectively, to map the control point motions to the computational mesh nodes.

The ith Bernstein polynomial of degree ℓ is:(12)

The ith B-spline of degree p is the non-periodic B-spline that uses a knot-vector, {t_i}, where:(13)The value of the B-spline is defined recursively as:(14)

Depending on the degree of the B-spline, the function is zero on some portion of the unit interval. This is in contrast to the Bernstein polynomials that are strictly non-zero everywhere on the unit interval.

Let denote the ith coordinate of the vth mesh node, and let denote the change in the position of the node due to the morphing process. Let denote the displacement of the (j, k)th control point associated with the Bernstein polynomials in the ith coordinate direction, and let denote the displacement of the (j, k)th control point associated with the B-splines. If (u ^v , v ^v ) denotes the position of the vth mesh node, then the displacement of the node is defined by the superposition:(15)

The control-point movement for the Bernstein polynomials controls the large-scale smooth deformation, while the control-point movement for the B-splines controls fine-scale motions.

4 Multiobjective optimization problem

A multiobjective optimization problem deals with more than one objective function and these goals are independent of each other. The task of multiobjective optimization is not to find an optimal solution corresponding to each objective function but to find a set of solutions called Pareto-optimal front.

Let consider the following multiobjective optimization problem (MOP):(16)where m is the objective functions number, f is the concurrent objective functions vector to optimise, x = (x₁, … , x_n) ∈ Ω is the D-dimensional decision space where each decision variable x_i is bounded by lower and upper limits x_li ≤ x_i ≤ x_ui for i = 1, … , D, g_i(x) are l inequality constraints and h_k(x) are p equality constraints.

Let us first introduce some basic concepts of the multiobjective optimization [32]:

• Pareto dominance: Let a = (a₁, ..., a_n) and b = (b₁, ..., b_n) be two vectors. The vector a dominates the vector b if and only if a is partially less than b, i.e.:

(17)

• Pareto-optimal solution: A solution a ∈ Ω is said to be Pareto-optimal if there is no solution b ∈ Ω which dominates a, i.e.

(18)

• Pareto-optimal set: The set Ps of all Pareto-optimal solutions, as defined by:

(19)

• Pareto front: The Pareto front, denoted P_F, is the image of Pareto optimal set Ps in the objective space and is defined as follows:

(20)

4.1 Multiobjective optimization based Backtracking Search Algorithm

In this paper, BSA, which was initially developed for single-objective optimization problems only, is extended to handle MOPs. This development maintain the BSA's spirit and allows a good efficiency with great benefit from new capacities of exploration of the search space and of an adapted search direction matrix due to the adaptation of BSA's main strategies and algorithm simplicity. The main development should concern the fitness values assigned to individuals in the population.

Algorithm 1 gives the proposed pseudocode of BSAMO (Backtracking Search Algorithm for Multiobjective Optimization). It integrates BSA's mutation and crossover operators presented in [24] and the fast non-dominated sorting and the crowding distance of Deb et al. [11].

Some explanations for this last item are given in the following subsections.

In the elaboration of the fast non-dominated sorting procedure, for every solution, and before creating the first non-dominated front and initialising it with all solutions having zero as domination count, the number of solutions which dominate the solution p, the domination count n_p, and the set of solutions that the solution p dominates S_p are calculated. Then, for each solution p with n_p = 0, each member q of its set S_p is visited, and its domination count is reduced by one. The second non-dominated front is then created as the union of all individuals belonging to Q, which is a separate list of any member with a domination count equal to zero. The procedure is repeated for subsequent fronts (F₃, F₄, etc.) until all individuals are assigned their ranks. The fitness is set to a level number, lower numbers correspond to higher fitness (F₁ is the best). The fast non-dominated sorting procedure was developed in the framework of NSGA-II [11].

The average distance between two individuals located on either side of the given particular solution along each objective is called the crowding distance. It is an estimation of the perimeter of cuboid formed using the nearest neighbours as mentioned in Figure 1. This metric represents half of the perimeter of the cuboid encompassing the solution i. The sum of individual crowding distance values corresponding to each objective gives the overall crowding distance value [11].

Based on the normalised values of objectives, the computation of the crowding distance is given by Algorithm 2, where and are the maximum and minimum values of the mth objective function, respectively.

is a non-dominated set, [i]_m is the mth objective value of the ith individual in , n is the number of elements of , and sort (, m) is sorting of the individuals of according to the mth objective. The theoretical aspect of this algorithm is developed in [16,24].

Fig. 1 Crowding distance of individual i.

5 Numerical simulation

5.1 Problem statement

In this work, a numerical simulation of the NACA0012 airfoil was evaluated two different Mach numbers, M = 0.5 and M = 0.7, same with the reliable experimental data from T.J. Coakley (1987) [33], in order to validate the present simulation. The flow can be described as compressible for this Mach numbers and the numerical simulation was conducted using FLUENT with a segregated implicit solver. Table 1 illustrates the parameters used in this numerical simulation.

In the appropriate boundary conditions of the fluid domain we will be setting “pressure far-field” for the inlet and the outlet where the gauge pressure, the Mach number, the temperature, and the components of the flow direction are given, and for the airfoil lower and upper a “wall” boundary condition is set and the velocity there will be 0.

The airfoil profile meshes were created using structured meshes, which consist of a variety of quadrilateral elements. The resolution of the mesh was greater in regions where greater computational accuracy was needed, such as the region close to the airfoil (Fig. 2). Generally, a numerical solution becomes more accurate as more nodes are used and the appropriate number of nodes can be determined by increasing the number of nodes until the mesh is sufficiently fine so that further refinement does not change the results.

5.2 CFD and validation results

This symmetrical airfoil is widely used for testing different computing methods and is considered as the main source of experimental data used here for comparison, and collected in the Experimental Data Base for Computer Program Assessment (AGARD-AR-138, 1979) for checking computational fluid dynamic methods. The data on the flow over the airfoil are obtained in several wind tunnels, but, according to Hoist (1987) [34], it is best to take the data from Harris's tests in the Langley 8-Foot Transonic Pressure Tunnel [35]. Computations were performed at a Reynolds number of Re = 9.10⁶. Comparison of the computed results with Harris's data for the pressure coefficient are shown in Figures 3 and 4 at M =  0.5 and M =  0.7, respectively, and at an angle of attack α = 1.55^∘. Table 2 pressents a comparision of the lift and drag coefficients with the experimental data.

Figures 3 and 4 show experimental and computed pressure and forces distributions at a Mach number of 0.5 and 0.7 and with angle of attack of 1.55°. In the first case The computations are in good agreement with experimental data but a significant difference can be seen at the upper airfoil near the chordwise location of x/c = 0.15. The computed lift and drag coefficients for those cases are compared with experimental normal force and drag coefficients in Table 2. In the second case (M = 0.7) there is a small amount of separation, and the computed and experimental pressure distributions are in close agreement with one another. The computations do, however, indicate the presence of a shock wave, near the chordwise location of x/c = 0.15, which causes the separations at the upper airfoil ; it is not apparent in the experimental pressure distrlbutions. It is found that the second case gives lift and drag coefficients that are slightly the same with the experimental values. Changing the turbulence model can also give significant variations in the predicted separations [33].

Fig. 2 Boundary conditions.

Fig. 3 Pressure coefficient along airfoil surfaces for Mach 0.5.

Fig. 4 Pressure coefficient along airfoil surfaces for Mach 0.7.

Fig. 5 Mesh morphing control points.

Fig. 6 Multiobjective optimization and FLUENT coupling process flowchart.

5.3 Optimization

The objectives of this optimization are to obtain a minimal drag coefficient of the airfoil with a maximal lift coefficient which leads to a maximum aerodynamic finesse ratio, implicitly given as below:(21)

It is required in all optimization problems to identify parameters that can be modified in order to reach the optimized solution. It is the geometry, in the case of the mesh morpher, that must be parameterized. The large variety of shapes available in engineering applications involves complications of the geometric parameterization for general shapes used in CFD. In order to minimize such complications in FLUENT simulation, the problem of shape parameterization is reduced to a problem of the parameterization of changes in the geometry.

The next essential requirement for mesh morphing is a tool that can smoothly alter the shape, irrespective of the underlying mesh topology. In FLUENT, this is accomplished using a free form deformation technique. This technique manipulates designated deformation regions via displacements applied to a set of control points. The mesh region that is to be deformed is defined by a “box” that is, a rectangle for 2D cases, and the control points must be located within the box. The displacements of the control points are the result of user-defined motions (each of which involves a parameter value and other directional settings) and these displacements are then applied to the mesh as a smooth deformation by either interpolating the displacement based on radial basis functions or using the tensor product of Bernstein polynomials [30].

The shape of the airfoil is defined by a spline. The shape parameterization is done with one parameter “Param” that prescribed the relative ranges of motion of the four select control points according to Figure 5. Two other input parameters are used for the optimization process that are the Mach number and the incidence as illustrated in Table 3.

Fig. 7 Pareto solutions of the NACA0012 airfoil.

5.3.1 FLUENT and MOP coupling process

The flowchart in Figure 6 shows the multiobjective optimization process of CFD problems. The optimization work was entirely automated via a main script written in MATLAB. First of all, starting the Multiobjective optimization process begins with the script generation of an initial population using a uniform random distribution function. Then multiple sets of decision variables (i.e., multiple geometries) are submitted and wait for each individual objective functions values before proceeding further. For changing the input parameters to compute the aerodynamic responses taken as objective values of the new shape, MATLAB calls FLUENT software and then it returns back those values to the optimization process for executing the multiobjective algorithm and displaying the optimal results. The coupling MATLAB/FLUENT is compiled using FLUENT as a Server session that is very similar to a conventional standalone FLUENT session, but with the addition of an Internet Inter-ORB Protocol (IIOP) interface exposed that accepts connections from suitable client applications, which is MATLAB is this paper. This capability is different from batch mode operation in that commands can be issued to the running session at any time rather than just from a predefined journal file. This allows solution steering and other manipulations without exiting from the FLUENT session.

Fig. 8 Pressure contour of an optimized airfoil.

5.3.2 Optimization results

The multiobjective shape optimization in this part is carried out for the airfoil NACA0012 by using BSAMO and NSGA-II algorithms coupled with a CFD solver according to the optimization process presented in Figure 6. For a population size and a maximum number of generations set respectively at 200 and 10, the Pareto fronts obtained for both algorithms (BSAMO and NSGA-II) are shown in Figure 7 and the results from Table 3 marked in bold present the best objective functions optimal values comparison found by the two algorithms for the airfoil.

Figure 8 displays the pressure contour of an optimized airfoil and non symmetrical compared with the original NACA0012 given in Figure 5. The freeform mesh morpher in FLUENT provides a powerful tool for arbitrary smooth changes in the geometry without being limited by a constrained parameterised geometry.

The results from Figure 7 and Table 4 indicate that the expansibility of the Pareto fronts distribution by BSAMO and NSGA-II is slightly the same when solving this optimization problem with a good quality distribution and approximation. Also for the two multi-physics simulations it obtains better optimal results for the objective functions.

6 Conclusion

In this paper a workbench project, including FLUENT, has been used to compute the transonic, compressible flow over a NACA0012 airfoil. The implicit density based solver with solution steering was employed and the computed results have been compared to published experimental data and good agreement was achieved. A case comparison has been carried out within CFD Post to compare the pressure fields at Mach 0.5 and Mach 0.7. Then the drag minimization and lift maximization study was carried out for the airfoil NACA0012 by using BSAMO and NSGA-II coupled with FLUENT. The more challenging Navier-Stokes computation was carried out to demonstrate the reliability and robustness of multiobjective optimization algorithms based on a main script written in MATLAB and its successful coupling with the CFD software. The free stream Mach number Ma, the angle of attack and the mesh morpher parameter were allowed to vary during the course of the optimization process. The convergence history of the computation was noted after about 2199 CFD calls, that the fitness value reaches its converged value. It should be mentioned here that the averaged fitness corresponds to the entire members in the generation and the maximum fitness corresponds to the best member in each generation.

The obtained results, compared to the NSGA-II algorithm, indicate that BSAMO is generally successful and competitive in dealing with complex and multi-physics systems. Therefore, we conclude that BSAMO offers a potential alternative solution for solving MOPs.

References

All Tables

All Figures

	Fig. 1 Crowding distance of individual i.
In the text

	Fig. 2 Boundary conditions.
In the text

	Fig. 3 Pressure coefficient along airfoil surfaces for Mach 0.5.
In the text

	Fig. 4 Pressure coefficient along airfoil surfaces for Mach 0.7.
In the text

	Fig. 5 Mesh morphing control points.
In the text

	Fig. 6 Multiobjective optimization and FLUENT coupling process flowchart.
In the text

	Fig. 7 Pareto solutions of the NACA0012 airfoil.
In the text

	Fig. 8 Pressure contour of an optimized airfoil.
In the text