Aerodynamic-Aeroacoustic Optimization of a Baseline Wing and Flap Conﬁguration

: Optimization design was widely used in the high-lift device design process, and the aeroacoustic reduction characteristic is an important objective of the optimization. The aerodynamic and aeroacoustic study on the baseline wing and ﬂap conﬁguration was performed numerically. In the current study, the three-dimensional Large Eddy Simulation (LES) equations coupled with dynamic Smagorinsky subgrid model and Ffowcs–William and Hawkings (FW–H) equation are employed to simulate the ﬂow ﬁelds and carry out acoustic analogy. The numerical results show reasonable agreement with the experimental data. Further, the particle swarm optimization algorithm coupled with the Kriging surrogate model was employed to determine optimum location of the ﬂap deposition. The Latin hypercube method is used for the generation of initial samples for optimization. In addition, the relationship between the design variables and the objective functions are obtained using the optimization sample points. The optimized maximum overall sound pressure level (OASPL) of far-ﬁeld noise decreases by 3.99 dB with a loss of lift-drag ratio (L/D) of less than 1%. Meanwhile, the optimized performances are in good and reasonable agreement with the numerical predictions. The ﬁndings provide suggestions for the low-noise and high-lift conﬁguration design and application in high-lift devices.


Introduction
With the rapid development of civil aviation, aeroacoustic problems were of wide concern among aviation manufacturing enterprises and researchers the world over. Strict noise airworthiness regulations are formulated by the International Civil Aviation Organization (ICAO) to limit the noise of civil aircraft. The research projects to reduce the noise level of civil aircraft were developed by the European Advisory Committee on Aeronautical Research (ACARE) [1] and the National Aeronautics and Space Administration (NASA) [2].
During takeoff and landing of civil aircraft, the flap lifting device is one of the main components that produce noise [3]. In the noise generated by flap lifting device, the side edge vortex shedding of the flap is the main noise source [4].
The research on the flap side edge noise source can be traced back to 1979, when Fink [5] and Schlinker [6] found that the flap side edge was an important aircraft noise source. Meecham et al. [7] confirmed the existence of severe surface pressure fluctuation in the flap tip region and demonstrated that the formation of pressure fluctuation is directly related to the coherent structures in the shear layer. Khorrami et al. [8] analyzed and showed that the unstable wave in the lateral shear layer is the main cause of noise by using the instability theory.
With the continuous progress of measurement technology, the understanding of flap side edge noise was deepened. Streett [9] and Radezrsky [10] et al. respectively carried out The benchmark three-dimensional wing and flap configuration is investigated in the current study, which was conducted in the same manner as the one studied experimentally by Aviation Industry Corporation of China Aerodynamics Research Institute (ARI). The experiment was performed in a 0.5 m anechoic wind tunnel, and the details about the wind tunnel were listed in Table 1.  As shown in Figure 2, the model was installed in the middle of the side plates. The height of the plate is 375 mm (same as the test section), and the side plate is 900 mm wide and 970 mm long. The 63-channel phased array was utilized for far field directivity analysis. It is multispiral, and it has a radius of 1m. In the experiment, the phased array was placed 1.8 m away from the wind tunnel centerline. The boundary conditions of the numerical simulation are obtained according to the experiment. The velocity of the freestream flow is U∞ = 50 m/s, the freestream temperature is T∞ = 300 K, and the static pressure is P∞ =101325 Pa. The periodic conditions are set to the spanwise boundaries. The spanwise length is set as 0.4 chord length of the flap, which is 40 mm.

Numerical Algorithm
The three-dimensional Large Eddy Simulation (LES) equations [17,18] coupled with dynamic Smagorinsky subgrid model are employed for numerical simulation of the aerodynamics on wing and flap configuration. The commercial CFD program ANSYS FLUENT 19.1 is used to simulate the flow field. LES method directly calculates the turbulent motion larger than the grid scale through the unsteady Navier-Stokes (NS) equations, while the As shown in Figure 2, the model was installed in the middle of the side plates. The height of the plate is 375 mm (same as the test section), and the side plate is 900 mm wide and 970 mm long. The 63-channel phased array was utilized for far field directivity analysis. It is multispiral, and it has a radius of 1m. In the experiment, the phased array was placed 1.8 m away from the wind tunnel centerline.  As shown in Figure 2, the model was installed in the middle of the side plates. The height of the plate is 375 mm (same as the test section), and the side plate is 900 mm wide and 970 mm long. The 63-channel phased array was utilized for far field directivity analysis. It is multispiral, and it has a radius of 1m. In the experiment, the phased array was placed 1.8 m away from the wind tunnel centerline. The boundary conditions of the numerical simulation are obtained according to the experiment. The velocity of the freestream flow is U∞ = 50 m/s, the freestream temperature is T∞ = 300 K, and the static pressure is P∞ =101325 Pa. The periodic conditions are set to the spanwise boundaries. The spanwise length is set as 0.4 chord length of the flap, which is 40 mm.

Numerical Algorithm
The three-dimensional Large Eddy Simulation (LES) equations [17,18] coupled with dynamic Smagorinsky subgrid model are employed for numerical simulation of the aerodynamics on wing and flap configuration. The commercial CFD program ANSYS FLUENT 19.1 is used to simulate the flow field. LES method directly calculates the turbulent motion larger than the grid scale through the unsteady Navier-Stokes (NS) equations, while the The boundary conditions of the numerical simulation are obtained according to the experiment. The velocity of the freestream flow is U ∞ = 50 m/s, the freestream temperature is T ∞ = 300 K, and the static pressure is P ∞ =101,325 Pa. The periodic conditions are set to the spanwise boundaries. The spanwise length is set as 0.4 chord length of the flap, which is 40 mm.

Numerical Algorithm
The three-dimensional Large Eddy Simulation (LES) equations [17,18] coupled with dynamic Smagorinsky subgrid model are employed for numerical simulation of the aerodynamics on wing and flap configuration. The commercial CFD program ANSYS FLUENT 19.1 is used to simulate the flow field. LES method directly calculates the turbulent motion larger than the grid scale through the unsteady Navier-Stokes (NS) equations, while the influence of small-scale vortices on large-scale vortices is reflected through subgrid model [19,20]. In LES method, two parts are divided for instantaneous variable ϕ, given by where ϕ is the large-scale average component, which is processed using a filter function ϕ = D ϕD(x, x )dx . D is the flow region, x is the spatial coordinate in the filtered largescale space, x' is the spatial coordinate in the actual flow region, and D (x, x') determines the scale of the solved vortex. ϕ represents the small-scale components. The unsteady NS equation is processed by the filter function D(x, x ), and the following can be obtained: ∂ρ ∂t where τ ij is the subgrid scale stress, referred to as SGS stress, which reflects the influence of small-scale vortex motion on the NS equation. According to Smagorinsky model, SGS stress was given by where µ t = (C s ∆) 2 S is the turbulent viscosity at the sublattice scale.
where ∆ x represents the grid size along the x-axis, and C S is Smagorinsky constant. The sound source term is calculated by unsteady CFD, and then the sound propagation characteristics are obtained by acoustic analogy method. The calculation methods of sound propagation characteristics mainly include Lighthill acoustic analogy equation [21], Curle's acoustic analogy equation [22], FW-H acoustic analogy equation [23], and Kirchhoff integral method [24]. The FW-H acoustic analogy equation is used in present paper.
The FW-H acoustic analogy method is based on the near-field flow solution, which is used as the sound source signal, and the aerodynamic noise at the far-field observation point is obtained by integrating the FW-H formula, given by where v n and u n are the normal velocity of the control surface and the normal velocity of the fluid, H(f ) is the Heaviside function, δ(f ) Dirac function, Appl. Sci. 2022, 12, 1063

of 19
The solution of FW-H equation can be obtained: The acoustic integral surface used in the integral formula can be any spatial surface containing solid. In this paper, the solid wing surface is taken as the integral surface, and the unsteady flow field information obtained by CFD simulation is taken as the input condition.

Numerical Validation
According to the experimental data performed by Bao [25] on pressure distribution, the numerical results obtained by CFD code are verified. Y+ is the main indicator for grid analysis, which plays an important role in turbulence modeling; it determines the height of the first cell normal to the wall. In this study, the height of the first cell is refined to obtain a Y+ value of approximately 1.0 to confirm the numerical accuracy. The X+ and Z+ estimations are 100 and 150 according to the length and width of the first cell, respectively. Similar grid scales can be referred to in the aerodynamic and aeroacoustics analysis of high-lift configuration using LES method [26]. The structural grids and the zoom-in view of the leading edge and trailing edge are shown in Figure 3. The parameters of mesh configuration are listed in Table 2. where vn and un are the normal velocity of the control surface and the normal velocity of the fluid, H(f) is the Heaviside function, δ(f) Dirac function, ( ) The solution of FW-H equation can be obtained: The acoustic integral surface used in the integral formula can be any spatial surface containing solid. In this paper, the solid wing surface is taken as the integral surface, and the unsteady flow field information obtained by CFD simulation is taken as the input condition.

Numerical Validation
According to the experimental data performed by Bao [25] on pressure distribution, the numerical results obtained by CFD code are verified. Y+ is the main indicator for grid analysis, which plays an important role in turbulence modeling; it determines the height of the first cell normal to the wall. In this study, the height of the first cell is refined to obtain a Y+ value of approximately 1.0 to confirm the numerical accuracy. The X+ and Z+ estimations are 100 and 150 according to the length and width of the first cell, respectively. Similar grid scales can be referred to in the aerodynamic and aeroacoustics analysis of high-lift configuration using LES method [26]. The structural grids and the zoom-in view of the leading edge and trailing edge are shown in Figure 3. The parameters of mesh configuration are listed in Table 2.     In the numerical simulation, the implicit dual time-stepping approach was used to solve the unsteady flow field. To predict noise with frequency below 10 kHz in noise analysis, and according to the maximum frequency calculation formula F = 1/2∆t, the time step ∆t was set to 5 × 10 −5 s.
Comprehensively considering the computational accuracy and efficiency of numerical simulation, according to the research experience of relevant researchers [27], the spanwise width of the model was set to 0.4 times of the flap chord length.
The porous surface has an important influence on the analysis of noise [28]. Turbulence and vortices in the wake region are mainly quadrupole sound sources, which have little effect on the far-field noise at low speed. In this paper, the influence of dipole sound source is mainly considered, and the noise integral surface is directly taken on the solid wall.
After the flow field is obtained by steady calculation, the mean values were obtained in last 5 s of the unsteady simulation process. The comparison of the mean pressure distributions along the surface between the numerical simulation and experimental results from the VALIANT database is presented in Figure 4. Figure 4a shows that there is a pressure rise at the upper trailing surface of the main wing, which is due to the block effect from the flap. Certain deviations could be observed here between the numerical results and experimental results, which could also be found in the research by Chalmers [29] and NUMECA [29]. Basically, the pressure coefficients on the surface of wing and flap from numerical simulations agree well with the experimental results.
Streamwise resolution along the wing 480 cells Streamwise resolution along the flap 150 cells Boundary layer resolution of the wing 60 cells Boundary layer resolution of the flap 60 cells Vertical resolution in the leading edge of the wing 60 cells In the numerical simulation, the implicit dual time-stepping approach was used to solve the unsteady flow field. To predict noise with frequency below 10 kHz in noise analysis, and according to the maximum frequency calculation formula F = 1/2Δt, the time step Δt was set to 5 × 10 −5 s.
Comprehensively considering the computational accuracy and efficiency of numerical simulation, according to the research experience of relevant researchers [27], the spanwise width of the model was set to 0.4 times of the flap chord length.
The porous surface has an important influence on the analysis of noise [28]. Turbulence and vortices in the wake region are mainly quadrupole sound sources, which have little effect on the far-field noise at low speed. In this paper, the influence of dipole sound source is mainly considered, and the noise integral surface is directly taken on the solid wall.
After the flow field is obtained by steady calculation, the mean values were obtained in last 5 s of the unsteady simulation process. The comparison of the mean pressure distributions along the surface between the numerical simulation and experimental results from the VALIANT database is presented in Figure 4. Figure 4a shows that there is a pressure rise at the upper trailing surface of the main wing, which is due to the block effect from the flap. Certain deviations could be observed here between the numerical results and experimental results, which could also be found in the research by Chalmers [29] and NUMECA [29]. Basically, the pressure coefficients on the surface of wing and flap from numerical simulations agree well with the experimental results.  The z-direction vorticity contours calculated by the numerical approach are shown in Figure 5. Small vortex structures form at the regions where the pressure coefficient is high, and the flow becomes more chaotic. The wake of the wing has a great influence on the upper surface of the flap, which makes the flow on the upper surface of the airfoil separate rapidly. Figure 6 shows the iso-surfaces of Q-criterion. The vortex structures on the upper and lower surfaces of the airfoil are significantly less than that on the flap, since the flow on the airfoil surface is less disturbed and there is basically no flow separation. The vortex structures formed on the flap surface are relatively smaller compared to those The z-direction vorticity contours calculated by the numerical approach are shown in Figure 5. Small vortex structures form at the regions where the pressure coefficient is high, and the flow becomes more chaotic. The wake of the wing has a great influence on the upper surface of the flap, which makes the flow on the upper surface of the airfoil separate rapidly. Figure 6 shows the iso-surfaces of Q-criterion. The vortex structures on the upper and lower surfaces of the airfoil are significantly less than that on the flap, since the flow on the airfoil surface is less disturbed and there is basically no flow separation. The vortex structures formed on the flap surface are relatively smaller compared to those over the main wing surface, indicating that the flap has more high-frequency noise sources relative to the main wing.
Microphones at 1.8 m from the middle of main wing on rotating support scan from 45 • to 135 • and 225 • to 315 • , while each set consists of 60 microphones with angle interval of 1.5 • . The location of the probes in the far field is arranged according to the experiment, as shown in Figure 7. over the main wing surface, indicating that the flap has more high-frequency noise sources relative to the main wing. Microphones at 1.8 m from the middle of main wing on rotating support scan from 45° to 135° and 225° to 315°, while each set consists of 60 microphones with angle interval of 1.5°. The location of the probes in the far field is arranged according to the experiment, as shown in Figure 7.  over the main wing surface, indicating that the flap has more high-frequency noise sources relative to the main wing.   over the main wing surface, indicating that the flap has more high-frequency noise sources relative to the main wing.   The overall sound pressure levels of the probes in the far field are shown in the figure below. According to Figure 8, the numerical simulation results agree well with experimental results, keeping the maximum error within 4 dB.  The overall sound pressure levels of the probes in the far field are shown in Figure 11. The experimental data only covers half of the far field noise (45~135°). The numerical predictions agree well with the experimental data, and the maximum error is kept within 4 dB. The overall sound pressure levels of the probes in the far field are shown in the figure below. According to Figure 8, the numerical simulation results agree well with experimental results, keeping the maximum error within 4 dB.  The overall sound pressure levels of the probes in the far field are shown in Figure 11. The experimental data only covers half of the far field noise (45~135°). The numerical predictions agree well with the experimental data, and the maximum error is kept within 4 dB. The overall sound pressure levels of the probes in the far field are shown in the figure below. According to Figure 8, the numerical simulation results agree well with experimental results, keeping the maximum error within 4 dB.  The overall sound pressure levels of the probes in the far field are shown in Figure 11. The experimental data only covers half of the far field noise (45~135°). The numerical predictions agree well with the experimental data, and the maximum error is kept within 4 dB. The overall sound pressure levels of the probes in the far field are shown in Figure 11. The experimental data only covers half of the far field noise (45~135 • ). The numerical predictions agree well with the experimental data, and the maximum error is kept within 4 dB. Appl. Sci. 2022, 12, x FOR PEER REVIEW 9 of 19 Figure 11. Sound pressure levels of far field probes.

Optimization Design Process and Results
The surrogate-based optimization strategy consists of design of experiments, surrogate models and optimization algorithm [30]. The high-precision CAA/CFD platform coupled with surrogate model can effectively reduce the computational cost and improve the optimization efficiency [31]. The overall framework of the optimization process is shown in Figure 12.

Design Variables and Optimization Objectives
In current study, the optimization objective is to reduce the noise in the far field while the aerodynamic characteristics are maintained.  Table 3 shows the design space.

Optimization Design Process and Results
The surrogate-based optimization strategy consists of design of experiments, surrogate models and optimization algorithm [30]. The high-precision CAA/CFD platform coupled with surrogate model can effectively reduce the computational cost and improve the optimization efficiency [31]. The overall framework of the optimization process is shown in Figure 12.

Optimization Design Process and Results
The surrogate-based optimization strategy consists of design of experiments, surrogate models and optimization algorithm [30]. The high-precision CAA/CFD platform coupled with surrogate model can effectively reduce the computational cost and improve the optimization efficiency [31]. The overall framework of the optimization process is shown in Figure 12.

Design Variables and Optimization Objectives
In current study, the optimization objective is to reduce the noise in the far field while the aerodynamic characteristics are maintained.  Table 3 shows the design space.

Design Variables and Optimization Objectives
In current study, the optimization objective is to reduce the noise in the far field while the aerodynamic characteristics are maintained. The optimization variables are the flap position and the flap deflection. The flap position parameters should be included in selecting design variables. Coordinates of the flap leading point (x_flap, y_flap) and angle of the flap (angle_flap) are chosen as the design variables. Table 3 shows the design space. The maximum overall sound pressure level (OASPL) of far-field noise is taken into consideration when selecting optimization objectives, as shown in Table 4. The lift-drag ratio of the initial shape is 20.014. The aerodynamic loss is less than 1% after optimization, so the constraint condition is set to the lift-drag ratio of 19.814, as shown in Table 5. Table 5. Optimization constraint condition settings.

Parameter Lower Limit Upper Limit
Lift-Drag Ratio (L/D) 19.814 N/A

Design of Experiments
The experimental design method is a mathematical method about how to arrange experiments scientifically and reasonably [32]. It is a technology of arranging experiments economically and scientifically on the basis of probability theory and mathematical statistics [33]. In the present study, according to surrogate method and the computational cost, 19 samples of the three design variables are generated by Latin hypercube method. The distribution of sample points is shown in Table 6 and Figure 13. The design variables of sample points are uniformly distributed in the design space.  The relationship between design variable and objective function can be clarified by the scatter diagram between design variable and objective function [34,35]. Related tables and correlation maps can reflect the relationship between the two variables [36], but it is difficult to accurately show the degree of correlation between the two variables. The correlation coefficient is used to reflect the degree of linear correlation between two variables, which takes the form as The scatter plots and correlation coefficients of three design variables (x_flap, y_flap, angle_flap) and lift-drag ratio (L/D) are given in Figure 14. In the scatter plots, the values of x-axis and y-axis represent the variation range of the two variables to be analyzed. The value of the correlation coefficient indicates the degree of linear correlation between the two variables, the closer the value of the correlation coefficient is to 1, the positive correlation between the two variables is shown. On the contrary, the closer the value of the correlation coefficient is to −1, the more the negative correlation between the two variables is shown. If the correlation coefficient is close to 0, the two variables have no linear correlation. In Figure 14a,b, the correlation coefficients are 0.2564 and 0.3190. The two design variables of x_flap and y_flap have positive correlation with lift-drag ratio (L/D). In Figure  14c, the correlation coefficient is 0.9862, indicating that the correlation between angle_flap and constraints is very strong. The relationship between design variable and objective function can be clarified by the scatter diagram between design variable and objective function [34,35]. Related tables and correlation maps can reflect the relationship between the two variables [36], but it is difficult to accurately show the degree of correlation between the two variables. The correlation coefficient is used to reflect the degree of linear correlation between two variables, which takes the form as

Var[X]Var[Y]
In which, Cov(X,Y) is the covariance of X and Y, Var[X] and Var[Y] are the variances of X and Y, respectively.
The scatter plots and correlation coefficients of three design variables (x_flap, y_flap, angle_flap) and lift-drag ratio (L/D) are given in Figure 14. In the scatter plots, the values of x-axis and y-axis represent the variation range of the two variables to be analyzed. The value of the correlation coefficient indicates the degree of linear correlation between the two variables, the closer the value of the correlation coefficient is to 1, the positive correlation between the two variables is shown. On the contrary, the closer the value of the correlation coefficient is to −1, the more the negative correlation between the two variables is shown. If the correlation coefficient is close to 0, the two variables have no linear correlation. In Figure 14a

The Surrogate Model
Kriging model [37] is an interpolation model for spatial modeling and prediction of random process/random field according to covariance function. Kriging model can obtain more accurate global optimal results than quadratic RSM model [38]. The model can be divided into regression model and correlation model. The regression model is the global approximation in the design space, and the correlation model can reflect the design spatial distribution structure of variables, which has an impact on the prediction accuracy of the model.
The Kriging model as a type of surrogate model can effectively establish a global mapping between design variables and optimization objectives, with lower computational cost. The Kriging surrogate model was widely used in aerodynamic design optimization of high-speed train [39] and high-lift airfoil design optimization [40].

Optimization Algorithm
The particle swarm optimization algorithm (PSO) [41,42] is a famous design optimization algorithm that has good global optimization performance. Figure 16 shows the basic flow of a modified PSO algorithm by the authors [43]: first, the position and velocity of particles in the initial particle swarm are obtained by initializing the particle swarm. After that, the nondominated rank of the population is sorted by particle fitness. Then, the crowding distance of each particle and niche count can be obtained. The better progeny populations could be obtained, and the population could be updated.

The Surrogate Model
Kriging model [37] is an interpolation model for spatial modeling and prediction of random process/random field according to covariance function. Kriging model can obtain more accurate global optimal results than quadratic RSM model [38]. The model can be divided into regression model and correlation model. The regression model is the global approximation in the design space, and the correlation model can reflect the design spatial distribution structure of variables, which has an impact on the prediction accuracy of the model.
The Kriging model as a type of surrogate model can effectively establish a global mapping between design variables and optimization objectives, with lower computational cost. The Kriging surrogate model was widely used in aerodynamic design optimization of high-speed train [39] and high-lift airfoil design optimization [40].

Optimization Algorithm
The particle swarm optimization algorithm (PSO) [41,42] is a famous design optimization algorithm that has good global optimization performance. Figure 16 shows the basic flow of a modified PSO algorithm by the authors [43]: first, the position and velocity of particles in the initial particle swarm are obtained by initializing the particle swarm. After that, the nondominated rank of the population is sorted by particle fitness. Then, the crowding distance of each particle and niche count can be obtained. The better progeny populations could be obtained, and the population could be updated.

The Surrogate Model
Kriging model [37] is an interpolation model for spatial modeling and prediction of random process/random field according to covariance function. Kriging model can obtain more accurate global optimal results than quadratic RSM model [38]. The model can be divided into regression model and correlation model. The regression model is the global approximation in the design space, and the correlation model can reflect the design spatial distribution structure of variables, which has an impact on the prediction accuracy of the model.
The Kriging model as a type of surrogate model can effectively establish a global mapping between design variables and optimization objectives, with lower computational cost. The Kriging surrogate model was widely used in aerodynamic design optimization of high-speed train [39] and high-lift airfoil design optimization [40].

Optimization Algorithm
The particle swarm optimization algorithm (PSO) [41,42] is a famous design optimization algorithm that has good global optimization performance. Figure 16 shows the basic flow of a modified PSO algorithm by the authors [43]: first, the position and velocity of particles in the initial particle swarm are obtained by initializing the particle swarm. After that, the nondominated rank of the population is sorted by particle fitness. Then, the crowding distance of each particle and niche count can be obtained. The better progeny populations could be obtained, and the population could be updated. Appl. Sci. 2022, 12, x FOR PEER REVIEW 13 of 19 Figure 16. Flow chart of modified particle swarm optimization (PSO) algorithm.

Optimization Results
After adding sample points for three times, the optimization process is finished according to the progressive optimization strategy. And the evolutionary process of optimization objectives is shown in Figure 17. Overall sound pressure level (OASPL) values are dispersed, and the average value is larger in the initial population. With the optimization, the population of the OASPL improves, and the average value tends to converge. Therefore, the optimization design is carried out smoothly. Comparing further the optimization in the second and third rounds, there is little change in the optimal population of the two rounds.

Optimization Results
After adding sample points for three times, the optimization process is finished according to the progressive optimization strategy. And the evolutionary process of optimization objectives is shown in Figure 17. Overall sound pressure level (OASPL) values are dispersed, and the average value is larger in the initial population. With the optimization, the population of the OASPL improves, and the average value tends to converge. Therefore, the optimization design is carried out smoothly. Comparing further the optimization in the second and third rounds, there is little change in the optimal population of the two rounds.

Optimization Results
After adding sample points for three times, the optimization process is finished according to the progressive optimization strategy. And the evolutionary process of optimization objectives is shown in Figure 17. Overall sound pressure level (OASPL) values are dispersed, and the average value is larger in the initial population. With the optimization, the population of the OASPL improves, and the average value tends to converge. Therefore, the optimization design is carried out smoothly. Comparing further the optimization in the second and third rounds, there is little change in the optimal population of the two rounds.     Figure 19. With the optimization proceeding, the evolutionary tendency of the population becomes smaller and converges.
After the three rounds of optimization design, the configuration parameters of the optimized shape are compared with those of the initial shape, as shown in Table 7. With the constraint of L/D ratio, there is little change on the flap deflection angle. After the change of the position of the flap, the width of the channel between the flap and the wing changes, so the effect of reducing noise is realized.
The optimization results are validated by CFD/CAA, as shown in Table 8. The errors of OASPL value and lift-drag ratio value between the two models are small, which verifies the accuracy of the surrogate model and the correctness of the optimization method.  Figure 18 shows the variation of the constraints in the optimization process. The liftdrag ratio values are dispersed and the mean value is small in the initial population. With the second round of optimization, most of the lift-drag ratio values are larger than the constraints and the mean value tends to converge.  Figure 19. With the optimization proceeding, the evolutionary tendency of the population becomes smaller and converges. After the three rounds of optimization design, the configuration parameters of the optimized shape are compared with those of the initial shape, as shown in Table 7. With the constraint of L/D ratio, there is little change on the flap deflection angle. After the change of the position of the flap, the width of the channel between the flap and the wing changes, so the effect of reducing noise is realized. The optimization results are validated by CFD/CAA, as shown in Table 8. The errors of OASPL value and lift-drag ratio value between the two models are small, which verifies the accuracy of the surrogate model and the correctness of the optimization method.  Figure 20 gives the converging process of the optimization objectives in the threeround optimization design. Comparing the results of aerodynamic performance and aerodynamic noise between the initial shape and the optimized shape, the maximum far-   Figure 20 gives the converging process of the optimization objectives in the three-round optimization design. Comparing the results of aerodynamic performance and aerodynamic noise between the initial shape and the optimized shape, the maximum far-field total sound pressure level value of the optimized shape decreases by 3.99 dB, and the lift-drag ratio value increases by 0.8%. The comparison of the far-field OASPLs before and after the optimization is given in Figure 21. After the three-round optimization design, the optimization results tend to converge and satisfy the constraints.

Discussions
Aiming at the change of flow field caused by the shape variation, the optimal design results are discussed. Figure 22 shows the comparison between the initial shape and the optimal shape. In the x-axis direction, the flap leading edge of the optimal shape is close to the leading edge of the main wing, and the flap leading edge of the optimized shape is downward, which widens the channel between the flap and the main wing, and the deflection angle of the optimal shape is almost unchanged compared with that of the initial shape.

Discussions
Aiming at the change of flow field caused by the shape variation, the optimal design results are discussed. Figure 22 shows the comparison between the initial shape and the optimal shape. In the x-axis direction, the flap leading edge of the optimal shape is close to the leading edge of the main wing, and the flap leading edge of the optimized shape is downward, which widens the channel between the flap and the main wing, and the deflection angle of the optimal shape is almost unchanged compared with that of the initial shape.

Discussions
Aiming at the change of flow field caused by the shape variation, the optimal design results are discussed. Figure 22 shows the comparison between the initial shape and the optimal shape. In the x-axis direction, the flap leading edge of the optimal shape is close to the leading edge of the main wing, and the flap leading edge of the optimized shape is downward, which widens the channel between the flap and the main wing, and the deflection angle of the optimal shape is almost unchanged compared with that of the initial shape.
The pressure coefficients on the airfoil surface before and after optimization are shown in Figure 23. For the upper surface of the main wing, the pressure coefficients of the optimal airfoil are almost the same as that of the original airfoil. At the lower trailing edge of the main wing, the pressure coefficients of the optimized airfoil decrease compared with the initial airfoil, and the area of the inverse pressure gradient decreases. Thus, the interference between main wing and flap is reduced. The flap has a limited influence on the wing, and there is no change on the pressure in front of the wing.

Discussions
Aiming at the change of flow field caused by the shape variation, the optimal design results are discussed. Figure 22 shows the comparison between the initial shape and the optimal shape. In the x-axis direction, the flap leading edge of the optimal shape is close to the leading edge of the main wing, and the flap leading edge of the optimized shape is downward, which widens the channel between the flap and the main wing, and the deflection angle of the optimal shape is almost unchanged compared with that of the initial shape. The pressure coefficients on the airfoil surface before and after optimization are shown in Figure 23. For the upper surface of the main wing, the pressure coefficients of the optimal airfoil are almost the same as that of the original airfoil. At the lower trailing edge of the main wing, the pressure coefficients of the optimized airfoil decrease compared with the initial airfoil, and the area of the inverse pressure gradient decreases. Thus, the interference between main wing and flap is reduced. The flap has a limited influence on the wing, and there is no change on the pressure in front of the wing. The z-direction vorticity contours and the iso-surfaces of Q-criterion (Q = 5000) of the initial shape and the optimal shape are shown in Figures 24 and 25. The difference between the flow field of the optimal shape and that of the initial shape is mainly reflected in the gap between the main wing and the flap. After optimization, the area of the inverse pressure gradient region on the main wing surface is reduced, and more air flow can pass by the passage between the main wing and the flap. The passage between the main wing and the flap becomes wider and the blocking effect of the flow is reduced and the separation area at the tail of the main wing is significantly reduced. Consequently, the aerodynamic noise is also reduced.  The z-direction vorticity contours and the iso-surfaces of Q-criterion (Q = 5000) of the initial shape and the optimal shape are shown in Figures 24 and 25. The difference between the flow field of the optimal shape and that of the initial shape is mainly reflected in the gap between the main wing and the flap. After optimization, the area of the inverse pressure gradient region on the main wing surface is reduced, and more air flow can pass by the passage between the main wing and the flap. The passage between the main wing and the flap becomes wider and the blocking effect of the flow is reduced and the separation area at the tail of the main wing is significantly reduced. Consequently, the aerodynamic noise is also reduced.
Appl. Sci. 2022, 12, x FOR PEER REVIEW 16 of 19 Figure 22. Comparison between initial shape and optimal shape.
The pressure coefficients on the airfoil surface before and after optimization are shown in Figure 23. For the upper surface of the main wing, the pressure coefficients of the optimal airfoil are almost the same as that of the original airfoil. At the lower trailing edge of the main wing, the pressure coefficients of the optimized airfoil decrease compared with the initial airfoil, and the area of the inverse pressure gradient decreases. Thus, the interference between main wing and flap is reduced. The flap has a limited influence on the wing, and there is no change on the pressure in front of the wing. The z-direction vorticity contours and the iso-surfaces of Q-criterion (Q = 5000) of the initial shape and the optimal shape are shown in Figures 24 and 25. The difference between the flow field of the optimal shape and that of the initial shape is mainly reflected in the gap between the main wing and the flap. After optimization, the area of the inverse pressure gradient region on the main wing surface is reduced, and more air flow can pass by the passage between the main wing and the flap. The passage between the main wing and the flap becomes wider and the blocking effect of the flow is reduced and the separation area at the tail of the main wing is significantly reduced. Consequently, the aerodynamic noise is also reduced.  (a) (b) Figure 24. z-direction vorticity contours: (a) initial shape; (b) optimal shape.

Conclusions
In this paper, the particle swarm optimization algorithm combined with Kriging surrogate model is applied