Aerodynamic Inverse Design of Transonic Compressor Cascades with Stabilizing Elastic Surface Algorithm

The upgraded elastic surface algorithm (UESA) is a physical inverse design method that was recently developed for a compressor cascade with double-circular-arc blades. In this method, the blade walls are modeled as elastic Timoshenko beams that smoothly deform because of the difference between the target and current pressure distributions. Nevertheless, the UESA is completely unstable for a compressor cascade with an intense normal shock, which causes a divergence due to the high pressure difference near the shock and the displacement of shock during the geometry corrections. In this study, the UESA was stabilized for the inverse design of a compressor cascade with normal shock, with no geometrical filtration. In the new version of this method, a distribution for the elastic modulus along the Timoshenko beam was chosen to increase its stiffness near the normal shock and to control the high deformations and oscillations in this region. Furthermore, to prevent surface oscillations, nodes need to be constrained to move perpendicularly to the chord line. With these modifications, the instability and oscillation were removed through the shape modification process. Two design cases were examined to evaluate the method for a transonic cascade with normal shock. The method was also capable of finding a physical pressure distribution that was nearest to the target one.


Introduction
Considerable efforts have been made to optimize the systems that include heat transfer or fluid flow in external and internal flows. The modeling limitations and computational costs are the major challenges when using design methods. Direct methods, inverse design methods, and different optimization algorithms have been used for optimal shape design. Direct methods are based on geometry corrections using feedback from a flow analysis code, and rely heavily on the previous designer's experience, which makes them time-consuming and inefficient because of the trial-and-error process [1].
Many studies have recently used direct numerical optimization for the aerodynamic shape design. Various optimization algorithms such as adjoint method, neural networks, evolutional algorithm or machine learning are coupled with a flow solver to reach the best cost function [2]. The optimization methods are often time-consuming or mathematically complicated with a high computational cost [3]. However, the freedom in selecting the multi-objective cost functions and applying the constraints are the benefits of the optimization methods [4,5].
In contrast to the direct methods, inverse design methods calculate the geometry corresponding to a given flow variable distribution on the boundary (e.g., wall target pressure distribution) [6]. By improving the aforementioned flow parameter distributions, the geometry profile can be inversely designed to increase the aerodynamic load, decrease the drag force, or decrease the flow losses [7]. Inverse design methods only calculate the unknown boundaries corresponding to the distribution of the target flow parameter along the boundaries, and they cannot obtain the optimal geometry [8]. From the viewpoint of engineering application, the inverse problem technique is desirable because it requires a much smaller number of flow simulations than the numerical optimization does [9]. In fact, inverse design methods make a direct relationship between geometry corrections and target flow parameter, which causes the faster convergence rate.
Coupled (non-iterative) and decoupled (iterative) methods are the two main types of inverse shape design methods. Also, indirect transpiration techniques and stream-functionas-coordinate trends are two classifications of coupled methods. In the formulation of coupled methods, the geometry coordinates are implicitly integrated with flow variables (such as pressure and velocity) [10,11]. Therefore, by solving the governing equations, the geometry shape is directly obtained along with flow variables. However, the complication of viscous flow equations and the complicated 3D geometries are the main challenges for using coupled design methods [12,13].
Stanitz [14] transformed the physical space into a computational space and developed an inverse design form of the Laplace equation in 2D potential flows. In this method, the desired velocity distribution is the target function, which is defined along the geometry boundaries. In addition, a coupled flow-geometry equation was generated so that the coordinates of equipotential lines, streamlines, and the values of flow angles were obtained by solving the coupled equation with suitable boundary conditions.
Zannetti et al. [15] converted the physical space into the computational space to present the inverse design form of 2D and axisymmetric Euler equations. Ashrafizadeh et al. [16,17] introduced the method of proper closure equations, which was a new finite-elementvolume CFD technique for the incompressible flow simulations on co-located grids. This method used appropriate physical equations to impose flow variables at integration points.
In iterative inverse design methods, the equations for flow simulation and geometry modification are solved separately. In fact, at each shape modification step, a series of consecutive problems are solved to modify the geometry. The main advantage of the iterative methods is that the flow analysis codes, such as commercial or open-source flow analysis software, which are developed in different regimes and complicated geometries, can easily be used to compute the required flow parameters on the target boundaries [18,19]. Iterative inverse shape design methods have always been the most common methods in applied shape-design problems [6].
There are several inverse design methods [20][21][22], which are capable of determining the proper aerodynamic shape to satisfy a target pressure distribution on its surface. Most of these methods require the development of new complex mathematical formulations and the accompanying new solvers, which are time-consuming and not suitable for industrial applications. Therefore, iterative inverse design methods were introduced to reduce the code development and to be integrated with any reliable CFD solver. In the following, several different types of iterative inverse design methods such as permeable wall, streamline curvature, flexible membrane or MGM, integral equations, average swirling velocity, potential flow-based, FSA, BSA and ESA methods are introduced.
In the permeable wall method [23], it is assumed that the geometry wall is porous and as a result, a normal penetration flux is defined on the wall. The cause of this flux is the difference between the target and current pressure distributions at each step. As the current geometry approaches the target geometry, the current pressure distribution approaches the target pressure distribution and consequently the penetration flux will be zero. The purpose of this method was to remove the velocity component normal to the permeable wall. Using this method, Demeulenaere et al. [23] performed the inverse design of axial turbine and axial compressor blades, assuming inviscid flow and solving three-dimensional Euler equations. They redesigned the transonic rotor blades of an axial compressor and the low-aspect-ratio blades of an axial turbine.
Veress et al. [24] applied the permeable wall method to redesign multi-stage radial compressors with the aim of achieving a more uniform Mach number distribution. In addition, based on this method and the use of pressure distribution as the blade loading, Dang et al. [25] presented their three-dimensional inverse design method by solving the viscous flow. The focus of their work was on the design of a transonic compressor blade with lateral suction on the blade. They showed that the use of this inverse design method reduces the blade suction mass flow.
Barger [26] proposed the streamline curvature method that relates the change in surface curvature to the velocity change. Since then, a number of methods have been developed based on this idea. This idea was used by Campbell [27] for the full potential equations, Bell [28] for the Euler equations and Malone [29] for the Navier-Stokes equations.
The flexible membrane method, which is an iterative mathematically based inverse design technique, was initially introduced by Garabedian and McFadden [30,31] as GM method. Then, Malone et al. [32][33][34] developed this method to the MGM (Malone-Garabedian-McFadden) method. The MGM method models the aerodynamic body wall as a membrane with a mathematical formulation that deforms owing to the difference between the desired and current pressures. Thus, the membrane is modeled by a partial differential equation (PDE), which converts to an ordinary differential equation (ODE) upon spatial discretization. The ODE has constant coefficients, and a force function along with the surface coordinate and its first and second derivatives determine the ODE coefficients. Therefore, this membrane model will be applicable. Moreover, the MGM method has been developed for the inverse design of wings, rotors and other 3D aerodynamic geometries [34].
Dulikravich [35,36] upgraded the MGM method using a solution for analytical Fourier series. The upgraded method provided appropriate results for subsonic and transonic, 2D and 3D flow regimes. In Dulikravich's design method, the aerodynamic surface was considered as a flexible membrane (the same as MGM method), so that the membrane characteristics and coefficients were determined by analytical Fourier series solution. The membrane was deformed under the difference between the desired and current pressure distributions, until the desired pressure distribution is achieved. A benefit of this method is that the shape modification solver is completely independent of the flow solver. Therefore, any proven flow solver (e.g., potential flow, Euler, Navier-Stokes, etc.) could be conjunct with the code [36].
Hazarika [37] developed an iterative inverse design method based on MGM for designing body-alone and wing-body configurations directly in the physical plane. This method was a general-purpose analysis-and-configuration-independent design scheme. In fact, it was not restricted to a specific geometry and despite using a full potential flow solver, the inverse method could be conjunct with either Euler and Navier-Stokes solvers.
Matsushima et al. [9] presented an iterative inverse design method using integral equations for design of airfoils and wings. This residual-correction inverse method was based on the small disturbance equation of velocity potential as the approximation of the perturbation of the Navier-Stokes equations. The inverse design method was applied in subsonic, transonic and supersonic external flow regimes. In the formulation of these inverse problems, the partial differential equations, which govern the flow field, were mathematically converted to an integral equation system.
The use of average swirling velocity distribution as a loading parameter is common in the inverse design algorithms, especially in turbomachines. The idea was first proposed by Hawthorne [38]. The advantage of this method in the design of turbomachines is that the required power of the machine is related to the rotor circulation [39].
Borges [40] and Zangeneh [41] used the average swirling velocity to design radial and mixed-flow impellers, assuming potential flow. Borges designed a low-speed radial turbine and made suggestions for choosing the optimal average swirling velocity distribution.
Zangeneh et al. [42] performed the inverse design of a centrifugal compressor with seven main blades and sevensplitter blades. In this design, the loading parameter was the meridian derivative of the average swirling velocity of the blade, which was directly related to the blade circulation. The blade thickness distribution was one of the input data and the geometry of the meridional plane remained unchanged. The results showed that the pressure ratio of the modified geometry increased by 4.5%.
In addition, Tiow et al. [43] redesigned the NASA rotor 67 transonic fan using a threedimensional inverse design method and applied the average swirling velocity distribution of the flow as a target function. They focused on the near-blade shock, the interaction of shock and boundary layer, and the tip leakage, and applied a three-dimensional viscous CFD solver based on the finite volume method.
G.S. et al. [44] performed the inverse design of the airfoil using the vorticity element method. They carried out the aerodynamic calculations for the inviscid potential flow, approximated the airfoil surface with vorticity sheets, obtained the velocity and pressure, and minimized the difference between the target pressure and the obtained pressure using a genetic algorithm.
Kamoun et al. [45] performed the inverse design for all the cross-sections of wind turbine blades, and obtained the velocity field at any point on the airfoil wall using the single element method, which is based on the distributions of sources and vortices on the airfoil wall. Finding the velocity field, the vortex distribution on the airfoil wall was determined. In the next step, the resulting vortex distribution was compared with the induced vortex distribution (determined by the designer) and the geometry was modified if the values were different. The above solution process was repeated until the difference in values became insignificant. Finally, according to the relationship between the lift force and the circulation, the amount of lift force was determined. They finally reduced the drag coefficient.
The Flexible String Algorithm (FSA), a physical iterative inverse design method for internal flows, was introduced by Nili-Ahmadabadi et al. [10]. In the FSA, the channel wall is modeled as a flexible string. The difference between desired and calculated pressure distributions deforms the wall until the desired pressure distribution is achieved. The convergence rate and stability of the inverse method are controlled by the string mass. The FSA was implemented for viscous incompressible [46] and inviscid compressible [47] internal flows.
The Ball-Spine Algorithm (BSA) was another iterative inverse method that was presented by Nili-Ahmadabadi et al. [7,48] in order to design the meridional plane of a centrifugal compressor in the quasi-three-dimensional space. The meridional plane walls were envisioned as a chain of virtual balls, which can move freely along certain directions named as spines. At each shape modification step, the difference between the desired and calculated pressure distributions is applied to the set of balls as a force vector that frequently deforms the wall. The BSA was developed for the inverse design of subsonic and transonic external-flow 2D airfoils [49].
Samadi Vaghefi et al. [50] incorporated the BSA into the FLUENT software [51] using user-defined functions (UDF) for the inverse design of an axisymmetric duct. Hoghoughi et al. [52] applied the BSA for the aerodynamic design of the nozzle of a wind tunnel to achieve a uniform airflow with minimum turbulence intensity and flow angle. Shumal et al. [53] developed the BSA to enhance the performance of an axisymmetric 90-degree bend duct with the swirling viscous flow regime. Mayeli et al. [54] used a proper physical-sense law in order to develop the BSA for inverse heat conduction problems. Hesami et al. [55] presented a new shape design method for the numerical solution of inverse heat convection problems (IHCPs) in nanofluids based on the BSA inverse method.
Madadi et al. [56] utilized the BSA for the quasi-3D inverse design of an S-shaped duct. In addition, they employed the BSA for the design of 2D cascade blades with sharp [57,58] and blunt [59] leading edges, and for the quasi-3D inverse design of the axialflow compressor rotor and stator blades [60]. They redesigned the rotor blade of an axial compressor with this method, so that the rotor loading coefficient was increased, which led to an increase of the compressor outlet pressure by 10%. Nili-Ahmadabadi et al. [61] developed a hybrid aerodynamic design method for curved diffusers using the genetic algorithm and BSA inverse design method. Kariminia et al. [62,63] developed a full 3D version of BSA for the aerodynamic design of 3D curved ducts.
Safari et al. [64] presented a novel iterative physical method for the inverse design of subsonic and transonic external-flow airfoils, which was called the Elastic Surface Algorithm (ESA). The ESA convergence rate in transonic flow regimes was much higher than that of BSA and MGM. The ESA models the airfoil wall as an elastic curved beam, which deforms owing to the difference between the desired (target) and calculated (current) pressure distributions in each modification step. Therefore, the inverse design process is implemented on a fluid-solid interaction basis. The innovation of this method is that the convergence rate and stability are controlled by physical characteristics such as the thickness and elasticity modulus of the beam. This property increases the robustness and simplicity of the ESA in comparison to mathematical-based methods such as MGM.
Noorsalehi et al. [65] applied the ESA for the inverse design of subsonic and transonic axial-flow compressor blades with separation. They applied the shear stress distribution as the target function to reduce the size and length of the separation zone. Noorsalehi et al. [66] developed the ESA for the inverse design in external separated flows by considering a linear combination of normalized pressure and shear stress distribution as the target flow parameter. They upgraded the ESA for separated flows by removing the geometrical filtrations and defining dynamic spines instead of the vertical spines. Nasrazadani et al. [67] upgraded the Elastic Surface Algorithm for the inverse design of 2D cascade blades, which have sharp leading and trailing edges and sharp pressure gradients near the stagnation point. In the upgraded ESA, the blade wall deformation is modeled with the deflection curve of a Timoshenko beam, which is differentiable and smooth at all the nodes. Therefore, unlike the original ESA, there is no need to use geometric smoothing for the fracture removal in the geometry profile after each modification step. Although the original ESA restricts the nodal displacements in the direction perpendicular to the chord line, the upgraded ESA eliminates this constraint to enhance the curved beam deformation based on the large deformation equations. Moreover, Nasrazadani et al. [67] extracted the optimized values of width, thickness and elasticity modulus of the curved beam to achieve a higher convergence rate of the upgraded ESA.
The original UESA had some limitations for transonic flow regimes with shock occurrence. When the shock on a transonic cascade is intensified due to increasing Mach number, there is a greater difference between the target and current pressure distributions near the location of shock occurrence, which leads to some oscillations of the elastic beam in the shape modification process. The movement of normal shock during the geometry corrections also causes oscillations and instabilities, which cause divergence of the inverse design.
The main goal of the current research is improving the Upgraded Elastic Surface Algorithm (UESA) for the compressor cascade in a viscous transonic flow regime with normal shock that can achieve severe pressure gradients in the shock region without geometry filtration. A distribution for the elastic modulus along a Timoshenko beam was chosen to increase the stiffness of the beam near the shock and to control the unstable deformation in this region.
In the following, the original UESA and the CFD solver are completely described. The steps of improving the UESA for transonic compressor blade cascade, especially the elastic modulus distribution, are explained. Then, the UESA is validated for the blade cascade with high subsonic and transonic flow regimes. Three design cases of transonic blade cascade with UESA are presented and the performance characteristics are compared before and after inverse design.

Original UESA
In the UESA, a cascade of sharp-edged blades is assumed as the primary (initial) geometry, and its walls are modeled as elastic curved Timoshenko beams. Then, an arbitrary CFD solver simulates the fluid flow over the 2D blades in order to compute the wall pressure distributions. The difference between the desired (target) and computed (current) pressure distributions is exerted to the blade walls, which are considered as elastic beams. The elastic beams deform until the internal stresses neutralize the exerted force vectors. After wall deformation at each step, the beam internal stresses should set to zero so that the elastic curved beams will have no residual stress at the beginning of the next modification step. This process is essential for convergence to the desired blades. Specifically, nonlinear finite-element equations based on Timoshenko's theory are applied to calculate the beam deformations [67,68]. The basic similar flowchart of the ESA and UESA is shown in Figure 1 in which the words TPD, CPD and ESM stand for "Target Pressure Distribution", "Current Pressure Distribution" and "Elastic Surface Method", respectively. The UESA is used for the inverse design of 2D cascade blades, which have sharp leading and trailing edges and sharp pressure gradients near the stagnation point. In the UESA, the blade wall deformation is modeled with the deflection curve of a Timoshenko beam, which is differentiable and smooth at all the nodes. Therefore, unlike the original ESA, there is no need to use geometric smoothing for the fracture removal in the geometry profile after each modification step.
Although the original ESA restricts the nodal displacements in the direction perpendicular to the chord line, the UESA eliminates this constraint to enhance the curved beam deformation based on the large deformation equations. Moreover, the optimized values of width, thickness and elasticity modulus of the curved beam are extracted to achieve a higher convergence rate of the upgraded ESA [67].
The following are the convergence criteria for the UESA method ( Figure 1) and beam finite element numerical code, respectively: In Equation (2), three degrees of freedom in the x, y, and θ directions are considered for each node. The symbol of ϕ shows an arbitrary degree of freedom. If the total number of curved beam nodes is n, the total number of degrees of freedom is 3n. If the spine constraint is used for node displacement, there are two degrees of freedom.

Numerical Solver for Fluid Flow
In the present research, the fluid flow passing over the blades was simulated using a compressible RANS (Reynolds-averaged Navier-Stokes) flow solver based on the conservation equations of mass, momentum and energy (Equations (3)-(5)). The flow over the blades was assumed two-dimensional, steady, compressible, and viscous.

∂ρ ∂t
where u i is the Reynolds-averaged velocity component with average sign omitted,ú i is the fluctuating velocity component, δ ij is the component of Kroneck tensor, Pr t is turbulent Prandtl number and µ e f f is the sum of turbulent viscosity and molecular viscosity, which is called effective viscosity. Moreover, −ρú iúj can be calculated by Boussinesq's hypothesis as follows: The turbulence model was k-ω-SST and the compressibility effects were considered for this model [69]. The conservation equations of k (turbulence kinetic energy) and ω (specific dissipation rate) are described in Equations (7) and (8).
where µ t (eddy viscosity or turbulent viscosity) is described as follows: The discretization of the equations was achieved using the pressure-based implicit finite-volume method with a coupled solver. The second-order upwind discretization scheme was applied to interpolate the variables stored at the cell centers to the faces. The gradients and derivatives inside the cells were computed from the Least-Squares Cell-Based relationship. The SIMPLE algorithm was used for a combination of continuity and momentum conservation equations (Equations (3) and (4)) in order to derive an equation for a pressure correction. In addition, the MUSCL (Monotonic Upstream-centered Scheme for Conservation Laws) scheme [70] was implemented for shock capturing scheme.
The non-slip boundary condition was used for the blade wall. Other solution and boundary conditions are presented in Table 1.

Improvement of the UESA for Compressor Cascade with Normal Shock
As described in the introduction, the original UESA had some limitations for transonic flow regimes with shock occurrence. When the shock on a transonic cascade is intensified due to increasing Mach number, there is more difference between the target and current pressure distributions near the occurrence location of shock, which leads to some oscillations of the elastic beam in the shape modification process. The movement of normal shock during the geometry corrections also causes oscillations and instabilities, which cause divergence of the inverse design. Two modifications are needed for the convergence of the inverse design problem to control these oscillations and instabilities in a transonic compressor cascade with high Mach number.

Constrained Movement of the Wall Nodes along the Spine
In the original UESA, to increase the convergence of the inverse design method, the wall nodes can move along both the X and Y directions. To obtain a unique solution for the inverse design problem, dynamic spines were defined. In the present work, for controlling the oscillations and instabilities and convergence of the inverse design problem, it is necessary to define constant spines for the node displacement direction through the shape modification process. In other words, the wall nodes can move only along the direction perpendicular to the chord line of the blade.

Using a Distribution of Elasticity Modulus along the Beam
In the new version of this method, a distribution for the elastic modulus along the Timoshenko beam is chosen to increase the stiffness of the beam near the shock and to control the high deformation in this region. FORTRAN code for the finite element equations of the beam deformation solves the matrix equation of KU = F. The force F is due to the difference between the target and the current pressure distribution, which is applied as a vector to the beam nodes. For each beam element, there are two degrees of freedom along the Y-axis and the rotation around the Z-axis (θ). The force vector at each node has one component in the Y direction with zero moment [68].
The beam stiffness matrix K is produced by super-positioning the stiffness matrices for all elements. The beam displacement matrix U is generated by applying the force matrix to the beam in the presence of the beam stiffness matrix. Thus, the pressure difference in one area of the beam affects the displacement of the beam in other locations. The distribution of E along the beam is used to solve this problem and to create more dependency between the pressure difference in one area of the beam and the displacement in the same area. Therefore, the node displacement decreases when increasing the elastic modulus near the locations of shock occurrence with a high difference between target and current pressure.
The decrease in average elastic modulus (Ē) is proportional to the decrease in the difference between the current and target pressure, which improves the convergence rate with no divergence and instability in the inverse design. The dimensionless parameter Q is thus defined according to Equation (10), and its value is constant and equal to Q in the first geometry correction during the inverse design process.
In Equation (10),Ē is the average elastic modulus of the beam, and ∆P int is calculated from Equation (11), in which P target ds * and P iteration ds * are integrals of the target and current pressure distributions at each geometry correction, respectively. The domain of these integrals is the dimensionless length of the airfoil (pressure or suction) wall. The continuous reduction ofĒ at the final geometry corrections causes the geometry profile to oscillate and be unstable. To prevent this instability, the reduction ofĒ was stopped after 100 reduction steps in the pressure difference.
In To validate an inverse design method, a target geometry is assumed, and the corresponding wall pressure distribution is considered as the target pressure. The validation of the inverse design method has been correctly performed if the arbitrary geometry used as an initial guess converges to the target geometry through the inverse design procedure. The original UESA [67] is validated for a cascade of DCA (Double-Circular-Arc) blades using the conditions in Table 1. The inlet Mach number is 0.87, and a structured grid is used for the numerical domain. The RANS equations are solved for each shape modification. Figures 2 and 3 show the generated grids and the contours of the Mach number for the cascade of DCA blades as the target geometry. The maximum local Mach number is 1.03. This means that the flow regime is mainly subsonic through the passage of the cascade.
A thin blade with 70% of the thickness of the target DCA blade is considered as the initial geometry, and AOA (Angle of Attack) is set to one degree. Figure 4 shows the procedure of geometry correction on both the suction and pressure surfaces of a cascade of sharp-edged blades. Although there is a big difference between the initial guess and target geometry, the initial guess geometry precisely converges to the target geometry after 150 iteration steps by the original UESA. It should be noted that for plotting the process of geometry correction, the blade is rotated around the origin so that the blade chord coincides with the X-axis. Then, the coordinates of rotated blade are normalized versus the chord length in order to achieve X * and Y * .    Figure 5 shows the residuals of the pressure difference during the geometry corrections in a high gradient of pressure difference without any profile fluctuation. Although the geometry completely converges to the target geometry, there are some fluctuations in the residuals of the inverse design due to a small region of transonic flow on the suction side of the blade. Later, it will be shown that these fluctuations will increase and lead to divergence of the original UESA result when increasing the inlet Mach number. Figure 6 shows how the current pressure distribution after 150 iteration steps coincides with the target pressure on both the suction and pressure surfaces. Figure 7 shows how the final pressure distributions on the suction and pressure side of the blade coincide with the target pressure distributions.

Validation of UESA for the Cascade with Transonic Flow (Mach Inlet = 0.96)
To have a transonic flow regime through the DCA cascade, the inlet stagnation pressure and outlet static pressure are set in a way that the inlet Mach number equals 0.96. With these conditions, the local Mach number exceeds 1 over the suction surface. Figure 8 shows the contours of the Mach number for the DCA cascade as the target geometry. The local Mach number reaches 1.34 on the suction surface. In the transonic regime, the pressure distributions of suction and pressure surfaces are strongly influenced by each other, which causes oscillations of the geometry and pressure distribution during the geometry corrections. In addition, the displacement of the shock occurrence region during the different geometry corrections increases the oscillations and instabilities. Therefore, for convergence of the inverse design process, two modifications are needed:

1.
The geometry correction of the pressure surface and suction surface should be independently carried out to control the oscillations. In other words, when the pressure surface is being corrected, the suction surface should be fixed, and vice versa. According to our experience with numerical inverse design in transonic regimes, after five iterations for the geometry correction of the pressure surface, five iterations for the geometry correction of the suction surface should be applied.

2.
By applying the constraint of the spine for nodal displacements (∆X = 0), the oscillations of the flexible wall through the inverse design process are significantly controlled.
After these two modifications, an inverse design was carried out as a validation case, in which a 70% scale of the target geometry was considered as the initial guess geometry. Figure 9 shows the procedure of the geometry correction on both the suction and pressure surfaces of the sharp-edged blade cascade by the UESA. The initial guess geometry converged to the target geometry after 175 iteration steps.  Figure 10 shows how the pressure distribution of the initial guess has coincided with the target pressure distribution on both the suction and pressure surfaces of the blade. This coincidence indicates flexibility of the method, especially near the leading edge with high pressure gradient and near the oblique shock with a Mach regime change. Figure 11 shows the residual of the pressure differences through the geometry corrections. There is no oscillation in the residual after 175 geometry corrections.  To show the effects of the modifications on the original UESA, the previous inverse design was repeated without using the constraint of ∆X = 0. Figure 12 shows that the residual of the pressure differences has high oscillations during the geometry corrections, while the oscillations were significantly removed by applying the spine constraint (∆X = 0), according to Figure 11. As mentioned before, the difference between the target and current pressures in subsonic flow is less than that in transonic flow. Thus, three degrees of freedom (in the x, y, and θ directions) are considered for each node displacement to increase the convergence rate of inverse design process in the subsonic flow. However, in transonic flow, due to the increase of pressure differences, especially near the shock occurrence region, this approach causes oscillations and even buckling in the geometry modification process. Therefore, reducing the degrees of freedom by one and applying the spine constraint (∆X = 0) eliminate these oscillations and reduce the number of geometry corrections. In all validation and design cases studied here, computations were performed using a PC with 2.3 GHz Intel ® Core™ i7 CPU and 8 GB RAM. CPU times, number of shape modifications (iteration steps) and the required residual for the convergence of flow solver have been shown in Table 2 for both cascade validation cases (Sections 2.4 and 2.5).

Design of Sharp-Edged Blades in Transonic Flow
In the validation of an inverse design method such as UESA (Sections 2.4 and 2.5), a known geometry is assumed as the target geometry, and the corresponding pressure distribution is considered as the target pressure distribution. An arbitrary geometry is then considered as the initially guessed geometry. The inverse design method is correctly verified if the initial pressure distribution converges to the target pressure distribution, and the initially guessed geometry converges to the target geometry.
However, for the blade aerodynamic design, the existing geometry is considered as the initial guess and its corresponding pressure distribution is modified and corrected (due to the designer's experience) to be considered as the target pressure distribution. During the inverse design process, the modified pressure distribution should be achieved. This design process is performed for a certain incidence angle (blade design condition). Finally, the performance characteristic curves of the initial and modified blades will be obtained by using a high-fidelity numerical solution in various incidence angles to compare the performance of initial and modified blades (Section 3.2).

Initial Design with Fixed Blade Stagger Angle
The design of the target pressure distribution along a transonic blade surface is very important. In this step, the designer modifies the pressure distribution of the current blade based on some simple concepts in fluid mechanics and considers the modified pressure distribution as the target pressure distribution. In transonic flow regimes, a supersonic flow region occurs in a part of the passage, which causes a shock on the blades. Displacement of the shock location through geometry corrections can cause the inverse design process to diverge, especially at higher inlet Mach numbers. Furthermore, it is hard to find a feasible pressure distribution for the blades, which increases the probability of divergence. To show this problem, an initial design for only the suction side of the DCA blade was carried out based on the details in Table 3, and the following constraints.

1.
Two degrees of freedom (Y, θ) were considered for node displacement, which means the nodes can only move in the spine direction (∆X = 0).

2.
The leading and trailing edges of the blade were fixed. Therefore, the blade stagger angle was fixed. The results of the initial design are shown in Figures 13 and 14. The plot of the pressure difference residual for this case in Figure 13 shows divergence after 70 iteration steps. However, because of the spine constraint for nodal displacements, there is no oscillation in the residual diagram. Figure 14 compares the current (initial), design (target), and final pressure distributions after 80 geometry corrections. According to this figure, the final pressure distribution is moving away from the design pressure distribution.

Second Design with Elastic Modulus Distribution and Fixed Blade Stagger Angle
The FORTRAN code for the finite element equations of the beam deformation solves the matrix equation of KU = F. The force F is due to the difference between the target and the current pressure distributions, which is applied as a vector to the curved beam nodes. For each beam node, there are three degrees of freedom along the X-axis, the Y-axis, and the rotation around the Z-axis (θ). The force vector at each node has two components in the X and Y directions, with zero moment. The stiffness matrix of the curved beam (K) is produced by super-positioning the stiffness matrices for all elements. By applying the force vector to the beam in the presence of the beam stiffness matrix, the beam displacement vector (U) is generated, so the pressure difference in one area of the beam affects the displacement of the beam in other locations. The distribution of E along the curved beam is used to solve this problem and to create more dependency between the pressure difference in one area of the beam and the displacement in the same area. For the convergence of this inverse design process, the following constraints were used:

1.
Two degrees of freedom (Y, θ) were considered for node displacement, which means the nodes could only move in the spine direction (∆X = 0).

2.
The leading and trailing edges of the blade were fixed. Therefore, the blade stagger angle was fixed.

3.
A distribution of the beam elastic modulus along the blade surface was implemented.
One of the constraints for the convergence of UESA is to use a distribution of the beam elastic modulus (E) prior to the shock region to design in the transonic flow regime. Due to the displacement of the shock location in different geometry corrections and the non-physical design pressure distribution, increasing E after the leading edge and then reducing it step by step along the blade surface prevents the design process from diverging. According to the finite element equations for the beam deformation, the elastic modulus can be considered a local property. The unique advantage of this method is that despite the beam being continuous, there is a chance of applying a local constraint to each part of the beam. Three different initial distributions of the elastic modulus ( Figure 15) were applied separately for the inverse design of the blade suction side. Other constraints are given above, and the design (target) pressure distribution is the same as in Section 3.1.1 ( Figure 14). According to Figure 15, the constant elastic modulus along the curved beam causes the inverse design process to diverge (Section 3.1.1). Applying the first and second elastic modulus distributions (in a stepwise descending manner) causes the convergence of the design process after 70 and 167 iteration steps, respectively. The first distribution has fewer iterations and a higher convergence rate, so it was used for the second design process. Figure 16 illustrates the current (initial), design (target), and final pressure distributions on the blade suction side for the second design case. The spine constraint for nodal displacements (∆X = 0) and the first elastic modulus distribution along the blade surface were considered. The geometry correction was applied to only the suction surface. The target pressure distribution was chosen so that the shock moves downstream and its intensity decreases. Accordingly, the final pressure converges to a distribution that is close to the target pressure distribution. The target pressure distribution is not exactly achieved because it is not physical for this fixed stagger angle. Actually, in the UESA inverse design method for a blade with normal shock, if the blade stagger angle is fixed, only the physical target pressure distributions will completely be achieved. However, if the blade stagger angle can change during the inverse design process, the target pressure distribution will be achieved anyways, which is the subject of Section 3.1.3. Figure 16. Current, design, and final pressure distributions versus non-dimensionalized X*coordinate, for the second design of the suction side of a transonic blade by UESA with the spine constraint and E distribution applied. Figure 17 shows the convergence history of the residual of pressure differences for this design case in a transonic flow regime. According to this figure, the residual of pressure differences almost reaches a constant value after about 70 iteration steps, which indicates the convergence of the inverse design process. The final pressure distribution in Figure 16 is a purely physical pressure distribution for the fixed stagger angle and can be used as a modified design pressure distribution in the transonic flow regime. Therefore, another inverse design process was performed with the modified design pressure distribution, as shown in Figures 18 and 19.    Figure 18 compares the current, final, and modified design pressure distributions on the suction side. According to this figure, the final pressure converges to the modified pressure distribution, which is a physical pressure distribution. Figure 19 illustrates the designed blade based on the modified design pressure distribution.

Third Design with Elastic Modulus Distribution and Low Variable Blade Stagger Angle
The target pressure distribution for the third design process is chosen to decrease the total pressure loss and to increase the blade loading (lift coefficient) in addition to moving the shock downstream (Figure 20a). The wall pressure before the shock for the third design is higher than that for the second design. The higher pressure before the shock weakens the shock, which leads to a loss reduction. However, the displacement of the shock is almost the same as in the previous design cases. Here are the constraints of the third design:

1.
Two degrees of freedom (Y, θ) are considered for node displacement, which means the nodes can only move in the spine direction (∆X = 0).

2.
Although the trailing edge of the blade is fixed, but the leading edge can freely move along the direction perpendicular to the blade chord. Therefore, the blade stagger angle is variable.

3.
A distribution of the beam elastic modulus along the blade surface is implemented. Unlike the previous design cases, the leading edge node in the third design case can freely move along the direction perpendicular to the blade chord, so the blade stagger angle can change slightly (Figure 20b). Figure 21 shows the leading edge constraint for the third design case, in which the leading edge node can freely move along a vertical groove. Figure 22 shows the geometry corrections in the third design process. The number of iteration steps is 110, and the leading edge is vertically displaced through the design process. The initial stepwise descending distribution of the elastic modulus along the suction side is shown in Figure 23.   CPU times, number of shape modifications (iteration steps) and the required residual for the convergence of flow solver have been shown in Table 4 for the last two design cases (Sections 3.1.2 and 3.1.3).

Comparison of Performance Curves
The designs of DCA blade cascade in Section 3.1 were carried out for a certain incidence angle (blade design condition). After that, the performance characteristic curves (including loss and deviation angle) of the initial and modified blades are obtained by using a high-fidelity numerical solution in various incidence angles, to completely compare the performance of initial and modified blades.
Theoretically, the errors in the CFD solution related to the grid must disappear with increases of mesh resolution. To obtain the performance characteristic curves of the designed blades and the current blade, a high-fidelity numerical solution with a finer grid is needed to calculate the accurate values of the flow losses and deviation angles and to determine the influence of the mesh size on the solution. It is also necessary to carry out a large number of high-fidelity numerical solutions with a wide range of inlet flow angles to calculate the loss and deviation angle in a wide range of blade incidence angles. For the grid study, the current blade cascade was numerically simulated with different grid sizes, and the flow losses were obtained. The loss in the 2D cascade was calculated from the following equation: where P 01 is the inlet total pressure, P 02 is the outlet total pressure, and W is the difference between the inlet and outlet total pressures. In addition, ρ and V 1 are the values of density and velocity in the cascade inlet, respectively. Figure 24 shows how the losses reach an asymptotic value as the number of cells increases.  Figure 25 shows the grid with 18,000 cells, which was chosen for the high-fidelity numerical solution of the current and designed blade cascades. The y+ value is less than five for this grid, which is proper for using the k-ω-SST turbulence model. The contours of the velocity for the current and two design cases are compared in Figures 26-28. For both designed blade cascades, the shock moves toward the trailing edge. The maximum Mach numbers in the third design case are lower than in the current geometry and second design case. The lowest maximum Mach number occurs in the third design case.    To obtain the performance curve of a blade cascade, inlet total pressure and outlet static pressure with periodic boundary conditions were applied to the numerical solutions. The inlet boundary condition for the inlet flow angle was also changed for each numerical solution to calculate the total pressure loss and the lift and drag coefficients at each incidence angle. Figures 29-32 compare the performance curves for the current blade cascade and each design case. A total of 13 inlet flow angles from 44 to 60 degrees were considered for calculating the performance curves.     Figure 29 shows the variations of the deviation flow angle and loss coefficient versus the inlet flow angle. The design point was placed at 80% of the maximum deviation angle. According to this figure, the design point is equivalent to the inlet flow angle of 53 degrees. Compared to the current case, there is a loss reduction for the inlet flow angles greater than 52.5 degrees and a loss increase for the inlet flow angles less than 52.5 degrees. A 2% loss reduction at the design point is observed for the second designed blade relative to the current blade. Figure 30 compares the variations of the lift and drag coefficient for the current and second designed blade cascades. According to this figure, as the inlet flow angle increases, the lift coefficient increases before reaching the stall point, and the drag coefficient significantly increases when approaching the stall point. In the second design case, the lift coefficient increases for flow angles more than 51 degrees compared to the current blade. For inlet flow angles from 53 to 57 degrees, the drag coefficient also decreases for the second design case. A 2% load increment (relative increase in the lift coefficient) at the design point is observed for the second designed blade relative to the current blade. Figure 31 compares the variations of the loss coefficient versus the inlet flow angle for the current and third designed blade cascades. The loss coefficient for the third designed blade relative to the current blade decreases at inlet flow angles of 52-55 degrees. For the third designed blade, there is a 7% decrease in the loss coefficient at the design point relative to the current blade. Figure 32 compares the variations of the lift and drag coefficients versus the inlet flow angle for the current and third designed blade cascades. The lift coefficient of the third designed blade increases for the inlet flow angles more than 52.5 degrees compared to the current blade cascade, indicating an increase in lift coefficient in this area. The drag coefficient for the third designed blade relative to the current blade decreases for inlet flow angles of 51.5-56.5 degrees. There is a 3.5% increment in the lift coefficient at the design point for the third designed blade (relative to the current blade).

Conclusions
In this study, the upgraded Elastic Surface Algorithm (UESA) was improved for the inverse design of a sharp-edged blade cascade in transonic compressors with normal shock. The displacement of the shock occurrence region during the geometry correction and the dependency of the pressure distributions along the pressure and suction surfaces, especially in transonic cascades, cause oscillations and instabilities for the geometry correction in the UESA method. To overcome these instabilities, the geometry corrections for the pressure and suction surfaces were carried out, independently. Furthermore, the constraint of the spine for node displacement (∆X = 0) was applied, and a distribution for the elastic modulus along the Timoshenko beam was chosen to increase its stiffness near the normal shock to control the large deformations and oscillations, especially near the shock. Finally, the target pressure distribution for a design case was chosen to decrease the total pressure loss and to increase the blade lift coefficient in addition to moving the shock downstream. The performance characteristics of the designed cascade were obtained by CFD analysis in design and off-design conditions and compared with those of the initial geometry, that showed the performance has been improved.