Maximum Thickness Location Selection of High Subsonic Axial Compressor Airfoils and Its Effect on Aerodynamic Performance

Solidity and camber angle are key parameters with a primary effect on airfoil diffusion. Maximum thickness location has a considerable impact on blade loading distribution. This paper investigates correlations of maximum thickness location, solidity, and camber angle with airfoil performance to choose maximum thickness location quickly for compressor airfoils with different diffusion. The effects of maximum thickness location, solidity, and camber angle on incidence characteristics are discussed based on abundant two-dimensional cascade cases computed through numerical methods. Models of minimum loss incidence, total pressure loss coefficient, diffusion factor, and static pressure rise coefficient are established to describe correlations quantitatively. Based on models, dependence maps of total pressure loss coefficient, diffusion factor, and static pressure rise coefficient are drawn and total loss variation brought by maximum thickness location is analyzed. The study shows that the preferred selection of maximum thickness location can be the most forward one with no serious shock loss. Then, the choice maps of optimal maximum thickness location on different design conditions are presented. The optimal maximum thickness locates at 20–35% chord length. Finally, a database of optimal cases which can meet different loading requirements is provided as a tool for designers to choose geometrical parameters.


INTRODUCTION
With the increase of stage loading, the area coverage of boundary layer separation enlarges and disordered vortex structures become more complicated within the flow passage (Sun et al., 2018;Hu et al., 2021;Ju et al., 2020). To achieve good aerodynamic performance, employing effective flow control approaches is necessary both in rotors and stators for the highly loaded aero-engine compressor . The low-reaction design method is proposed to simplify the sophisticated flow control strategy (Qiang et al., 2008;Sun et al., 2019;Sun et al., 2020;Sun et al., 2021). For a highly loaded low-reaction stage, both the inlet Mach number and the static pressure rise of the stator significantly increase. In this case, the high subsonic stator airfoil design with high flow turning which can achieve low losslevel, increase efficiency and static pressure rise capability becomes a huge challenge (Zhang et al., 2019).
When inlet Mach number increases to a certain degree, supersonic flow appears at local region of suction surface and brings sharp rise of loss (Shi and Ji, 2021). It is called supercritical flow condition when incoming flow is high subsonic and transonic flow exists in the passage. With shock-free and controlled diffusion features, controlled-diffusion airfoil (CDA) introduced in the early 1980s was originally designed for supercritical cascades (Dunker et al., 1984;Niederdrenk et al., 1987;Elazar and Shreeve, 1990;Steinert et al., 1991;Jiang et al., 2015). As CDA blades were primarily developed for aero-engines but not the optimum solution for heavy-duty gas turbines, Köller et al. (2000) and Kusters et al. (2000) designed a high subsonic compressor airfoil family (Ma ≤ 0.8, θ≤ 30°, e/b 25%). The optimization results demonstrated that the well-controlled front loading airfoil obtained good performance through promoting early transition without strong laminar separation. Besides, it is conducive to reducing the adverse pressure gradient and the loading at the rear region of suction side. Sieverding et al. (2004) transferred and adapted the preceding design system to industrial compressors and the optimized airfoils (Ma 0.5, e/b 25%) were generated. The results indicated that as the maximum thickness location moved forward, the peak velocity position on the suction side migrated forwards, and stabilized boundary layer is the reason why front-loading airfoils present good performance. Wang et al. (2020) investigated the evolution of vortex structures and separated flow transition process in frontloaded airfoil through large eddy simulations. As the laminar separation was suppressed on the suction surface of the frontloaded airfoil, the associated vortex dynamics was weakened. Hence, the high level of Reynolds shear stresses recedes.
Although previous conclusions for front-loading airfoils are effective, most optimization research objects are dispersive and lack of revealing the universal rules of parameter selection. When facing different design requirements, it is difficult to achieve a quick front-loading blade profile design. Moreover, these questions must be raised: What kind of flow problems will occur when the maximum thickness position moves forward? Is there a limit to the forward movement of the maximum thickness position? What is the effect of cascade key parameters on the influence brought by forward movement of maximum thickness position? The systematic study of cascade geometrical parameters is helpful to reflect the influence rules of key geometrical parameters on airfoil aerodynamic performance and reveal the flow evolution mechanism. In this situation, it is necessary to build up performance prediction models for supercritical airfoils, which include maximum thickness location and key geometrical parameters. The results of models contribute to revealing the generality of maximum thickness location selection and providing a tool of blade profile selection for mean-line and through-flow design.
For high-load aero-engine compressor design, solidity and camber angle are key geometrical parameters in preliminary design, which have a remarkable influence on flow turning capacity, blade loading, and working range (Howell, 1945;Zweifel, 1945;Carter, 1950;Lieblein et al., 1953;Johnsen and Bullock, 1965;Britsch, 1979). An increase of solidity could achieve high flow turning, but increases the number of blades and profile loss (Sans et al., 2014). To meet the needs of empirical solidity inputs in mean-line and through-flow design, performance prediction models including solidity which can be used in multistage axial flow compressor are needed (Larosiliere et al., 2002;Horlock and Denton, 2005;Burberi et al., 2020). Sans et al. (2014) updated the correlations for controlled diffusion blades and extended their application to the high Mach number flow regimes. Bruna and Cravero chose inlet flow angle, inlet Mach number, AVDR, Reynolds number, and solidity as model parameters to study profile loss (Ma1 0.5-0.84, σ 0.625-2, θ 37.4°-43.4°) (Bruna et al., 2006). Xu and Du established performance characteristics prediction models for the curved blade containing solidity and camber angle in NACA65 compressor cascade (Ma1 0.2-0.6, σ 1.2-1.8, θ 30°-60°) (Xu et al., 2018), and integrated aspect ratio into the models in the follow-up research (Xu et al., 2019).
As an attempt to provide a reference for compressor preliminary design, this paper investigates rules of maximum thickness location selection matching cascade key parameters. Several camber angles are selected to investigate different typical flow problems for the generality of the conclusion. The in-house program introduced in reference (Li, 2007) is used for the generation of profiles. The numerical calculation carried out in abundant two-dimensional linear cascades has been successfully applied for performance evaluations and flow phenomena detections. The evolution of performance parameters (ω o , D o , and Cp 2o ) at minimum loss incidence versus geometrical parameters (σ, e/b, and θ) is studied. The flow mechanism leading to the performance differences is analyzed based on the calculated performances and flow details. Then performance parameters are described in function of geometrical parameters. Based on polynomial models, dependence maps ofω o , D o , and Cp 2o and the effect of e/b on them are discussed. In addition, choice maps of optimal e/b at selected θ, D o , or Cp 2o are presented respectively. Finally, a database of all the optimal e/b is provided for designers with different demands as a tool to choose parameters.

Investigation Programs
The investigated blade profile is used in the stator of the second stage of a three-stage compressor. The curve of the camber line is described by a quartic polynomial and normalized by axial chord length, as shown in Figure 1A. It is well controlled to arrange the blade loading chordwise distribution and to alter the deviation angle. Two cubic polynomial curves connected at maximum thickness location are used to describe the thickness distribution, as shown in Figure 1B. The d 1 , d 2 are the leading edge and trailing edge radius. The variables K 1 , K 2 are the rate of thickness changing along the chord at leading edge and trailing edge. The e is the location of maximum thickness, while e/b is normalized by chord.
Geometrical parameters and flow conditions of the cascade are shown in Table 1. Maximum thickness and its position are normalized by chord length. Chord length, maximum thickness and inlet curvature of camber line has been optimized in the previous design of blade profile. Cascade geometrical parameters (σ, e/b, and θ) are changed in the ranges which cover the usually used values in industry applications. Effects of e/b and θ on profiles are shown in Figure 2. Figure 3 shows the axial distribution of the thickness and suction side curvature of the blade at θ 40°. As shown in the figure, the curvature of the suction side has a step, whose position is obviously related to e/b. The suction side curvature increases before the step but reduces after it with the forward migration of e/b.

Numerical Method and Validation
The structured grid is generated by the IGG/AutoGrid5 module of the NUMECA software package. As shown in Figure 4, a mixed O-and H-block is used for the computation domain. The inlet of the computational domain is located 100% chord length ahead of the leadingedge. The outlet of the computational domain is set 100% chord length behind the trailing-edge with a buffer zone. The measurement plane location is 40% chord length behind the trailing-edge for complete wake development. The thickness of the first grid layer is 10 −6 m, and 5 × 10 −6 m near the hub and shroud, which can insure y + <3 at solid wall in all cases and satisfy the solving demands of the selected turbulence model in this paper.
With such a topology structure, a grid independence study is carried out with different grid nodes number, including seven meshes with respectively from about 45 to 140 thousands cells. As shown in Figure 5, when further refining the mesh grid with grid number beyond 100 thousands, the results almost have no modification. Therefore, the spanwise node number adopted is 5 in this study, and the total grid node number is about 109 thousands.
The code solver FINE ™ EURANUS is used to solve the Reynolds-Averaged Navier-Stokes equations with finite volume form. The governing equation for the relative velocities in the rotating frame can be expressed as: where U is the vector of the conservative variables; F I → is inviscid flux vector; and F V → is viscous flux vector. Their variable forms are given by:  where the superscripts "-" and "∼" are the time average value and density weighted average value, respectively; w i is the x i component of the relative velocity; E is the total energy; p* is the total pressure; δ ij is the Kronecker number; τ ij is the Reynolds stress; and q i is the heat flux component. The contributions of the Coriolis and centrifugal forces are contained in the source item vector Q, which is given by: where ω is the angular velocity of the relative frame of reference; r is the radius; and w is the relative velocity. (NUMECA International, 2015). The MOGA profile in literature (Song and Ng, 2006) is selected to validate the calculation methods of the twodimensional cascade. Several turbulence models were adopted to show the effect of turbulence models on the simulation accuracy at Ma 0.77 in Figure 6A. It can be seen that the trend of numerical simulation is basically consistent with experiments. The Spalart-Allmaras turbulence model is employed, as it provides high robustness and good compromise in accuracy for the prediction of complex flows (Spalart and Allmaras, 1994). The simulations are fully turbulent without boundary layer transition effects.
Besides, the space discrete form is TVD of the second-order upwind scheme, and Min-Mod flux limiter is adopted to enhance the scheme stability. The temporal discretization scheme is fourstage Runge-Kutta scheme, with multigrid solver and implicit residual smoothing method to accelerate the residual convergence.
Inlet boundary conditions are prescribed with total pressure (430 kPa), total temperature (440 K) and inflow angle. Free-slip and adiabatic conditions are imposed at the solid walls. Translational periodic boundary condition is used in the pitchwise direction. In the computations, the mass-flow specified at the outlet is modified to achieve the desired inlet Mach number. Non-reflecting boundary condition is used at inlet/outlet planes.

RESULTS AND DISCUSSION
Based on Lieblein's definition (Lieblein, 1955), the diffusion factor given in Eq. 4 is used to assess the loading of blades. In Eq. 4,v 1 ,v 2 ,Δv u and σ are inlet velocity, outlet velocity, tangential velocity variation, and solidity. Total pressure loss coefficient (ω) and static pressure rise coefficient (Cp 2 ) are selected to evaluate airfoil aerodynamic performance. The definitions of ω and Cp 2 are described in Eq. 5 and Eq. 6. Here, P p 1 andP 1 are total and static pressure at the cascade inlet, whileP p 2 andP 2 are total pressure and static pressure at the outlet. Mass-weighted averaging is used.
The effect of maximum thickness location (e/b), solidity (σ) and camber angle (θ) on airfoil performance are discussed, including minimum loss incidence (i o ), diffusion factor (D o ), total pressure loss coefficient (ω o ), and static pressure rise coefficient (Cp 2o ). Minimum loss incidence (i o ) refers to the incidence with minimum loss within the incidence range, while the other three performance parameters are obtained in this case. Figure 7 shows incidence characteristics of total pressure loss coefficient and diffusion factor from negative to positive stall at different solidities (Ma 1 0.8). It can be seen that as solidity increases, total pressure loss coefficient curves move up and rightwards, and diffusion factor curves migrate rightwards. Obviously, i o , ω o and D o increase with solidity. Figure 8 shows incidence characteristics of total pressure loss coefficient and diffusion factor with different camber angles (Ma 1 0.8). It can be seen that as camber angle increases, incidence range shrinks, while both total pressure loss coefficient and diffusion factor curves move up, which means ω o and D o increase but i o changes slightly. Figure 9 shows incidence characteristics of total pressure loss coefficient and diffusion factor with different maximum thickness locations (Ma 1 0.8). It can be seen that when maximum thickness location moves upstream, positive incidence range increases and loss reduces, except for It can be found that solidity (σ), camber angle (θ), and maximum thickness location (e/b) have impacts on twodimensional cascade incidence characteristics. The trend of optimum incidence performance varying with cascade geometrical parameters will be discussed in detail afterward.

Minimum Loss Incidence
Supplementary Figure S10 shows i o versus σ with different e/b and θ. It can be seen that as σ increases, i o rises nearly linearly, except for the cases at e/b 0.2, θ ≤ 40°. The trend is in accord with the results in reference (Sans et al., 2014). Supplementary Figure S11 shows i o versus θ with different σ and e/b. It shows that the evolvement of i o versus θ is quadratic almost for all e/b and σ. In most cases, i o decreases slightly with an increasing θ, except for cases at e/b 0.2, σ ≥ 1.46. In these cases, i o decreases more significantly, which agrees with Supplementary Figure  S10A. Supplementary Figure S12 shows i o versus e/b with different σ and θ. It illustrates that the evolution of i o with e/b is quadratic, reflecting in two different patterns. In the cases with big θ and small σ, i o increases firstly and then decreases while e/b is enhanced. In the cases with small θ and big σ, i o decreases as e/b increases.

Total Pressure Loss Coefficient
Supplementary Figure S13 shows ω o versus σ with different e/b and θ. It illustrates that ω o almost increases linearly with σ. It is because that a rise in a number of blades attracts linear increases of profile loss. Supplementary Figure S14 shows ω o versus θ with different e/b and σ. As θ rises, ω o decreases firstly and then increases, which can also be observed in Supplementary Figure  S13. Supplementary Figure S15 shows ω o versus e/b with different σ and θ. It indicates that ω o decreases firstly and then increases with the rise of e/b, which presents quadratic curves except for some low-load cases. The minimum of ω o occurs at e/b 0.2-0.3 and alters with σ and θ.

Diffusion Factor
Supplementary Figure S16 shows D o versus solidity with different e/b and θ. As solidity is increasing, D o increases and then keeps constant with a quadratic trend. The influence of σ is very limited. It is because that the rise of i o with σ causes an increase of blade load correspondingly and compensates for the Frontiers in Energy Research | www.frontiersin.org November 2021 | Volume 9 | Article 791542 5 decrease of blade load caused by σ. Another observation is that D o increases faster for lower cambers at e/b 0.2, which is different from the reference (Sans et al., 2014). It is because that the e/b is different from that of the literature (Sans et al., 2014), while the observation in literature (Sans et al., 2014) agrees with the cases at e/b ≥ 0.3.
Supplementary Figure S17

Static Pressure Rise Coefficient
Supplementary Figures S19-21 show Cp 2o versus σ, θ, and e/b respectively. The trends of Cp 2o evolution resemble these of D o in Supplementary Figures S16-18. Besides, it can be noticed that Cp 2o is influenced more significantly than D o by σ.

Flow Mechanism Analysis of Normal Airfoils
This section aims at revealing how e/b affects the loss-level of the blade and transonic flow of airfoils with normal camber angle (θ ≤ 40°) without flow separation. Analysis indicates that the related flow mechanism is the same at θ 35°and θ 40°. For brevity, the cases with θ 40°are selected to discuss here. As shown in Supplementary Figure S22, the e/b can cause a significant influence on velocity distribution at blade surface and load of the blade. With the upstream propagation of e/b, the position of suction peak velocity moves forward, and peak velocity increases, which result in the front-loaded velocity distribution. When e/b moves forward to a certain extent, the flow after the suction peak decelerates sharply. Besides, suction peak velocity increases with the reduction of solidity, which means a rise in blade loading.
As shown in Supplementary Figure S23, the flow separation does not occur in moderate-load blade passage. It can be observed through the flow field that the acceleration in the front portion of the passage amplifies with the upstream propagation of e/b. Shock wave generates when e/b moves forward to some extent. It indicates that the forward migration of e/b induces shock. In conjunction with Figure 3B, it can be attributed to that in the flow acceleration region before the shock (before 20% chord), the forward migration of e/b causes the increase of the curvature of the suction side, which is conducive to the flow acceleration. By comparing the flow field at different solidities, it can be found that the shock is easy to occur at low σ, which is caused by the increased suction peak velocity. In general, the forward migration of e/b and reduction of σ induces the generation of shock.
Supplementary Figure S24 illustrates the velocity gradient on the blade surface, which reflects viscous shear in the boundary layer. It can be seen that in the case without serious shock, the velocity gradient reduces with the forward migration of e/b in the most portion (30-80% chord), which indicates the weaker viscous shear. It is the reason why the small e/b reduces total pressure loss in the case without shock. By comparing with Figure 3B, it can be found that in the deceleration area after the suction peak velocity, the reduction of suction side curvature caused by the upstream propagation of e/b contributes to the lower gradient of velocity. In conjunction with analysis of flow field, how the e/b alteration influences loss of the blade is demonstrated: in the flow acceleration region before 20% chord, the forward migration of e/b induces the generation of shock; while in the deceleration region after the suction peak velocity, the upstream propagation of e/b contributes to the weaker viscous shear.

Flow Mechanism Analysis of High-Load Airfoils
This section illustrates how e/b affects the loss of the blade and transonic flow of high-load airfoils with higher camber angle (θ ≥ 44°) with flow separation. For brevity, the cases with θ 46°are selected to discuss for the clearer difference. As shown by Supplementary Figure S25, the trend of velocity distribution on the blade surface and load of the blade caused by e/b and σ are the same as normal airfoils (θ ≤ 44°, as shown in Supplementary Figure S22). The case e/b 0.31 is the optimal e/b at θ 46°, σ 1.25. It should be mentioned that the rise of the camber angle generates a slight drop in suction peak velocity.
Supplementary Figure S26 illustrates the flow field at θ 46°. It can be observed that compared with normal airfoils, flow separation exists at the rear portion of the blade (behind 80% chord). The scale of separation shrinks with the upstream migration of e/b. Besides, the effects of e/b and σ on  acceleration flow are the same as that in normal airfoils. It should be noted that the strength of supersonic flow at θ 46°is weaker than that of θ 40°, which can be attributed to the lower peak velocity in Supplementary Figure S25. The influence of e/b on boundary layer behavior is used to explain the reason for loss-level. It can be observed in Supplementary Figure  S27 that the influence of e/b on velocity gradient and viscous shear agrees with the trend reflected in Supplementary Figure S24. Especially, reverse flow occurs in the area where the velocity gradient is close to zero. At small solidity (σ 1.25), the forward migration of e/b delays the onset of boundary layer separation.
Supplementary Figure S28 reflects the boundary layer growth before separation. It can be noticed that low e/b contributes to the stabilized boundary layer, which reflects in the delay of the onset of rapid growth of boundary layer. The total pressure distribution of blade wake shown in Supplementary Figure S29 also confirms that the high total pressure loss region of blade wake significantly shrinks with the forward migration of e/b. In general, small e/b conduces to the drop of the loss (caused by viscous friction) and the separation scale of high-load airfoils, although contributes to the generation of shock.

Models
Based on the discussion above, it can be inferred that i o , ω o , D o and Cp 2o change with σ, θ, and e/b, and their correlations can be described by such quadratic multiple regression model as Eq. 7.
The variableX x i denotes cascade geometrical parameter (σ, e/b, or θ). The variableNdenotes the number of input parameters (3). The variableFdenotes performance parameters (i o , ω o , D o , or Cp 2o ). Variables a、b、c denote model coefficients. The regression models are verified by the coefficient of determination R 2 ≥ 0.96.
Supplementary Figure S30 shows the comparison of results obtained by models and CFD calculations, in which the horizontal axis and vertical axis represent CFD results and model predicting results. The red lines imply that the results predicted by models are equal to CFD calculations, while blue dots are close to them.
Supplementary Figure S31 shows errors of the model predicting results compared with CFD calculations. The error of minimum loss incidence (Δi o ) is defined as

Loss Variation on Isosurfaces of Camber Angle
The established models help to obtain numerous continuous data in the investigation sample space and to explore the correlations maps of parameters. The distribution of minimum loss (ω o ) is given in Supplementary Figure S32 when θ is constant. Three axes represent e/b, θ, and σ, while the blue color represents low ω o . It can be seen in Supplementary Figure S32 that ω o rises sharply when θ increases over 44°. For constant θ, the low loss area exists at low σ and upstream e/b. Because the camber angle is determined in premier compressor design, it is necessary to explore the correlations of e/b and σ at selected θ.
Quantitative analysis of total loss variation brought by e/b is necessary in order to evaluate the impact of e/b. It can be found in Supplementary Figure S32 that for constant θ and σ, the maximum loss is mostly at e/b 0.6. The total loss variation (δω o ) is defined as δω o ω o,e/b 0.6 −ω o ω o,e/b 0.6 . While ω o is the local total pressure loss coefficient, ω o,e/b 0.6 is the total pressure loss coefficient at e/b 0.6 with the same σ and θ.
Supplementary Figure S33 displays the distribution of total loss variation (δω o ) caused by the alteration of e/b (θ 40°), as well as contours of D o and Cp 2o . Three axes of Supplementary Figure S33 represent σ, e/b, and ω o , while the red color represents high total loss variation. For each σ, there is a maximum δω o (namely, the minimum ω o ) case, whose e/b is called as optimal e/b (e/b opt ). Such cases compose the black line named as (δω o ) max, θ 40°.
In Supplementary Figure S33, contours of D o and Cp 2o at constant θ are displayed under the contour of δω o . The solid black lines and dashed ones represent contours of D o and Cp 2o respectively. It can be seen that the upstream move of e/b leads to an increase of D o and Cp 2o , which means that the optimum e/b reduces loss and increases diffusion.
All the cases with the maximum of δω o composing the surface of e/ b opt for selected θ and σ are given in Supplementary Figure S34. The points on the surface produce lower loss than those with the same θ and σ but off the surface. It can be observed that e/b opt locates at 25-35% chord. This observation includes but is not limited to the reference (Koller et al., 1999;Kusters et al., 2000;Sieverding et al., 2004).
The distribution of optimal e/b is illustrated in Supplementary  Figure S35, while δω o of e/b opt is shown in Supplementary  Figure S36. It can be observed that e/b opt moves downstream when θ or σ decreases and the shock generates more easily. It can achieve more than a 10% reduction of total pressure loss at θ ≥ 44°. Thus, compared with normal-load airfoils, smaller e/b is preferred to produce excellent aerodynamic performance in high-load cascades, while total loss variation caused by e/b opt is more significant.
The cases with e/b opt (θ 40°or 46°, σ 1.25 or 2.05) in Supplementary Figure S35 are selected to verify the accuracy of e/b opt predicted by the model. Their velocity distribution and boundary layer behaviors are depicted in details in Supplementary Figures S22-29. In conjunction with the previous analysis, it can be found that the e/b opt is the most forward e/b without serious shock loss. Therefore, the preferred selection of e/b for the high subsonic airfoil can be suggested as the most forward e/b with no serious shock loss.

Loss Variation on Isosurfaces of Diffusion Factor
Because diffusion factor indicates blade load, it is necessary to explore correlations of parameters at different In Supplementary Figure S38, contours of θ and Cp 2o are drawn as solid black lines and dashed ones under the contour of δω o . It can be seen that the upstream move of e/b leads to θ and Cp 2o reduction. It means that e/b opt achieves lower ω o and Cp 2o and maintains D o with lower θ than e/b 0.6 at the same σ.
All the cases with e/b opt at selected σ and D o among the sample space are given in Supplementary Figure S39. The points on the surface produce lower loss than those with the same D o and σ but off the surface. It can be observed that the e/b opt locates at a 20-35% chord.

Loss Variation on Isosurfaces of Static Pressure Rise Coefficient
The situation where Cp 2o is constant is discussed in this part. Because Cp 2o rises with the increase of D o , the isosurfaces of Cp 2o in Supplementary Figure S40 are similar to those of D o in Supplementary Figure S37. Solidity has a significant effect on Cp 2o , which causes the difference between isosurfaces of D o and Cp 2o . Since highest loss occurs at e/b 0.6 for selected Cp 2o and σ, δω o is defined as δω o ω o,e/b 0.6 −ω o ω o,e/b 0.6 , in which ω o is the local total pressure loss coefficient and ω o,e/b 0.6 is the total pressure loss coefficient at e/b 0.6 with the same σ and Cp 2o .
Supplementary Figure S41 displays δω o and ω o versus σ and e/b when Cp 2o is constant, as well as contours of D o and θ. The black line which is named as (δω o ) max, Cp 2o 0.4 consists of the cases with maximum δω o at each σ. At selected Cp 2o , the optimal e/b moves forward as σ increases.
In Supplementary Figure S41, the solid black lines and dashed ones under the isosurface of Cp 2o are contours of θ and D o . Because θ reduces and D o increases as e/b decreases, e/b opt achieves lower loss and higher D o with lower θ than e/b 0.6 at the same Cp 2o and σ.
The cases with e/b opt are given in Supplementary Figure S42 and produce lower loss than those with the same Cp 2o and σ but off the surface. The e/b opt locates at a 20-35% chord.

Database of Optimal Airfoils
All the optimal e/b that can meet the three different kinds of design requirements (selected θ, D o , or Cp 2o ) are contained by a database, whose flow chart is shown in Supplementary Figure S43. The design requirement is regarded as a restriction, and correlation of optimal e/b, ω o , D o , Cp 2o , and θ with σ on the line of (δω o ) max are depicted in maps. According to needed performance parameters (such as ω o or Cp 2o ), σ can be chosen appropriately and input. Then corresponding e/b opt and performance parameters (ω o , D o , Cp 2o ) predicted by models are provided as the output.
The program interface of the database is shown in Supplementary  Figure S44. When the θ input is 44°, maps of e/b opt , ω o , D o and Cp 2o versus σ are shown below the θ input box in Supplementary Figure  S44A. It can be seen that as σ increases, e/b opt moves forward, and ω o , D o , Cp 2o rise. When the desired σ (such as 1.3) is input, corresponding e/b opt , ω o , D o and Cp 2o are indicated by red stars on the curves, and their values are displayed on the right side. Small σ (for example, σ 1.25) can be chosen to reduce ω o , when e/b opt is at 32% chord length. Large σ (for example, σ 2.05) is available to improve D o and Cp 2o , when e/b opt is at 23% chord length.
The database is also useful when the design goal is to achieve a certain D o (for example, D o 0.45) or Cp 2o (for example, Cp 2o 0.4). When D o is selected, the effects of σ on e/b opt , ω o , θ and Cp 2o are given below the D o input box in Supplementary Figure S44B. It can be seen that as σ increases, e/b opt moves forward, the corresponding θ decreases, while ω o and Cp 2o rise. When Cp 2o is selected, the effects of σ on e/b opt , ω o , θ and D o are presented below the Cp 2o input box in Supplementary Figure S44C. As σ increases, e/b opt moves forward, corresponding θ and D o decrease, and ω o rises. Therefore, the predetermined Cp 2o can be achieved at σ 1.25 with the lowest ω o , when e/b opt is at 33% chord length.

CONCLUSION
This paper aims to investigate the influence of geometrical parameters (e/b, θ, σ) on the aerodynamic performance (D o , ω o , Cp 2o ) of high subsonic stator airfoil. Numerical investigation indicates that optimization of e/b can lead to a considerable reduction of loss level at minimum loss incidence condition. The key conclusions are as follows: 1) Substantial two-dimensional cascade cases are simulated through validated numerical methods. It is indicated that D o , ω o , and Cp 2o versus σ are linear, but quadratic versus e/b. D o and Cp 2o versus θ are linear, while ω o versus θ is quadratic. 2) In the flow acceleration region before the 20% chord, the forward migration of e/b causes the increase of the curvature of the suction side. It is conducive to the flow acceleration and induces the generation of shock. While in the deceleration region after the suction peak velocity, the reduction of suction side curvature caused by the upstream propagation of e/b contributes to the drop of loss and separation scale. Therefore, the preferred selection of e/b can be suggested as the most forward e/b with no serious shock loss. 3) Quadratic multiple regression models of performance parameters are established. The models help to obtain the optimal e/b surfaces at selected θ, D o or Cp 2o . The results show that when e/b opt locates at 20-35% chord, δω o is maximum. As θ, D o , or Cp 2o increase, δω o becomes more significant. So, the selection of e/b has an important effect on high-load blade profile. 4) A database is established to meet design requests, which provides optimal e/b according to the selected θ, D o , Cp 2o and σ. When θ is selected, maps of performance parameters (ω o , D o , Cp 2o ) versus σ are depicted, according to which σ can be appropriately chosen. When selected σ is input, corresponding optimal e/b and performance parameters predicted by models are presented as the output. The database also helps when the design request is selected D o or Cp 2o .

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding author.

AUTHOR CONTRIBUTIONS
CT provides innovation, numerical simulation experiment, model establishment and writing for this article. XD provides financial support and guidance about the content of the article. JD and ZW contributes similar to XD. YL has engaged in some technical work.