An Inverse Design Method for Cascades for Low-Reynolds Number Flow

Inverse design of cascades for low Reynolds number can be applied for the aerodynamic development of fans and compressors. The present contribution describes a complete design procedure by taking into account the transition from laminar to turbulent boundary layer flow. A shape factor distribution is prescribed along the suction surface of the blades. The inverse boundary layer calculation is performed by the application of a finite difference method. On the pressure side the velocity distribution is prescribed in such a way that the given flow angles in front of and behind the cascade are realized. An inverse calculation based on potential theory is applied in order to determine the geometry of the cascade. At the end of the present contribution a cascade is designed by the described inverse design procedure and the flow is simulated by the application of CFD.


Introduction
The design of highly loaded blades for compressors and fans is still a challenging task for engineers.Particularly, it is of interest to reduce the flow losses and noise.This is the reason why design engineers are looking for design tools, which make it possible to take into account these requirements of fluid machines.There are three design strategies: i use of design data based on measurement and/or experience, ii design by fast recalculation, and iii inverse design.
The first mentioned method, for example, using NACA65-data 1 , works quickly and robustly, but is restricted to special profile classes and therefore cannot be used to reach an optimal design.A more advanced design method is the design by fast recalculation.The geometry of the cascade is parameterized and the flow values are calculated quickly on PCs by the application of more or less simple flow models.Nevertheless, the fast recalculation method contains the disadvantage that the designer can never be sure to reach an optimal design.In order to overcome this problem the application of inverse design is advantageous.Inverse design is based on the idea to prescribe flow values and to calculate the geometry in such a way that the prescribed flow values are realized.The present paper is restricted to inverse design.
First inverse design methods were developed in the second half of the last century when computer efficiency was increasing.Nowadays a commercial software tool for inverse design is available on the market.The software is based on the work of Zangeneh 2 .One of the first studies concerning inverse design was done by Schwering 3 .He published a work about the design of cascades for incompressible flows with prescribed velocity distribution.The method is based on singularity methods and does not take into account viscous effects.
Goto et al. 4 describe a method whereby a convenient velocity distribution can be calculated along the suction surface by using an inverse boundary layer method.In Goto's work an integral method of Thwaites 5 is applied inversely for the laminar flow region.The shape factor is determined by a parameter m introduced by Thwaites.m is chosen in such a way that the boundary layer does not separate.In the turbulent region a quarter sine wave is prescribed conveniently for the shape factor along the suction surface and by the application of the inverse integral method of Head and Patel 6 the velocity at the outer edge of the boundary layer is calculated.Goto et al. 4 recommend a value of the shape factor H at the trailing edge of H TE 1.7 . . .1.9.
The idea of Goto et al. 4 to prescribe a quarter sine wave of the shape factor along the turbulent boundary layer is advantageous, because the slope dH/ds at the trailing edge is dH/ds TE 0 and therefore boundary layer separation is avoided.The present contribution is based on the approach of Goto et al. 4 .There are some differences in comparison to 4 .These differences are as follows. 1 A finite difference scheme is applied in order to perform the inverse boundary layer calculation.The advantage is that more elaborate models can be used, especially for turbulence modeling.2 The velocity distribution along the pressure side is prescribed by the use of polynomials and done in such a way that a given circulation is realized.3 The inverse potential calculation is done by the application of Schwering's method 3 .
Finally, a result of inverse design is discussed in this contribution.It is shown that there is a good agreement between the prescribed and the realized velocity distribution calculated by the application of potential theory.Both velocity distributions are compared with RANS results and it is pointed out that inverse design is an appropriate method for designing cascades.In order to accomplish the discussion the performance of a NACA65-cascade is compared with the performance of an inverse designed cascade.

Concept of the Method
It is assumed that the flow angles β 1 , β 2 and the spacing ratio t/l are given for the cascade see Figure 1 .The geometry of the cascade is to be calculated by inverse design, that is, the shape of the blades and the stagger angle λ must be determined.
The velocity distribution at the outer edge of the boundary layer is determined by an inverse design method in order to avoid boundary layer separation.A shape factor distribution is prescribed along the suction surface based on an idea of Goto 4 see Figure 2 .
The definition of the shape factor is 7 Following the idea of 4 the suction surface is divided into four regions.In the first region nose region the outer flow is prescribed by and the values of the boundary layer are calculated; that is, for the first region no inverse boundary layer is performed.By formula 2.4 a typical accelerating flow around the nose is described.At s 0 the velocity has a prescribed value C B .The velocity C B is explained in more detail when the velocity distribution on the pressure side is introduced.
In the second region, the laminar flow region, the shape factor H is constant and chosen in such a way that the boundary layer does not separate.The velocity at the outer edge is calculated by an inverse boundary layer method which is described in the next section.Thwaites's results 5 are used to choose conveniently the shape factor in the laminar region.The transition region is characterized by a rapid decrease of the shape factor H s .Again an idea of Goto is applied to describe the shape factor by a half cosine wave as follows: H lam and H turb stand for the shape factor in the laminar and at the beginning of the turbulent region, respectively.They must be chosen conveniently.
In the turbulent region the shape factor H s is prescribed by a quarter sine wave as follows: Index TE stands for the according values at the trailing edge.Goto recommends that H TE should be chosen as H TE 1.7 . . .1.9 4 .At the trailing edge the derivative dH/ds TE is zero which is appropriate in order to avoid flow separation.
The reason why the condition dH/ds TE 0 is appropriate can be explained as follows, at the location where the flow separates the shape factor becomes large see Figure 3 .Some authors state that flow separation occurs, if H exceeds a limit, for example, H 2.5.A more accurate criterion of flow separation is dH/ds TE → ∞ see Figure 3 .So the opposite of dH/ds TE → ∞, that is, of dH/ds TE 0, is an appropriate input condition for avoiding flow separation.
On the basis of the prescribed shape factor the velocity distribution at the outer edge of the boundary layer is calculated for the suction surface by an inverse boundary layer method which is described in the next section.
Along the pressure side of the profiles the velocity distribution is prescribed by polynomials in such a way that the circulation is large enough in order to realize the prescribed turning Δβ β 1 − β 2 see Figure 4 .In the nose region the velocity distribution is prescribed by a can be chosen arbitrarily, for example, a 3, . . ., 5, and is a possibility to have an influence on the shape of the profile nose.The same function has the coordinate x stag .By this coordinate the location of the stagnation point can be prescribed see Figure 4 .The circulation Γ is calculated by The values of the end points of the polynomials are varied until the prescribed circulation is realized see black points in Figure 4 , that is, the following equation must be fulfilled.r stands for a coordinate running in a clockwise direction around the profile.At the very beginning of the design procedure the geometry of the profiles is not known there is no geometry .An arbitrary uncambered profile is taken as an estimation and along the coordinate s the shape factor is prescribed on the suction surface.The inverse design procedure consisting of an inverse boundary layer calculation and inverse application of potential theory is used in order to find a better estimation of the profile.This procedure is repeated until the changes of the geometry fall under a prescribed limit, that is, until the changes of the geometry become smaller as |Δy/y| < , 10 −7 .

Boundary Layer Model
The boundary layer flow is modeled by solving the boundary layer equations:

3.1
μ and μ t stand for the molecular and turbulent viscosity, respectively.μ t is zero in the laminar region.In the transition region the turbulence is modeled by an intermittency factor γ which describes the temporal part of the flow with turbulence 8 .γ is a function of s.At the end of the laminar region γ is equal to zero and at the beginning of the turbulent part γ ≈ 1.The boundary conditions for 3.1 are

3.2
In the case of inverse boundary layer calculation the boundary conditions are formulated as follows:

Laminar Boundary Layer
The boundary layer equations are solved by the application of a finite difference scheme 9, 10 .For region 1 see Figure 2 the finite difference scheme is applied directly, that is, C s is prescribed and the velocity vectors u n , v n T are calculated.An inverse application is performed for the remaining regions.The velocity at the outer edge of the boundary layer C s is estimated and varied until a prescribed shape factor is realized.This process is supported by a variable secant Newton method 10 .
The finite difference scheme is see Figure 5 ρ

3.4
The boundary conditions are u i,j 0 0, v i,j 0 0, for n 0 u i,j max C 3.5 and for inverse application u i,j 0 0, v i,j 0 0, for n 0 H s H prescribed .

3.6
At location i the velocity profile is known and further downstream at location i 1 the velocities can be calculated by the application of the finite difference scheme 3.4 .The application of 3.4 leads to a tridiagonal equation system for u i 1,j j 0, . . ., n max − 1 .By knowing u i 1,j j 0, . . ., n max − 1 the velocity components v i 1,j j 0, . . ., n max − 1 are ISRN Applied Mathematics calculated by the continuity equation by taking into account that v i 1,0 0. The discretisation of the continuity equation is The inverse boundary layer calculation is performed as described below i At location with index i the velocity profile is known.It fulfills the boundary condition u i,j 0 n 0 0, v i,j 0 n 0 0, and H s i H prescribed .
ii C i 1 is estimated.By the application of the finite difference scheme 3.4 the velocities u s i 1 are calculated.The boundary conditions u i 1,j 0 n 0 0 and C C i 1 are taken into account.
iii The velocities v i 1,j are calculated by 3.7 .The boundary condition v i 1,j 0 n 0 0 is taken into account. iv So the procedure starts at the second step again and is repeated until is a small positive value with 10 −6 .The procedure is supported by a variable secant Newton method 10 .
As a rule at each location up to five iterations are necessary.
The inverse laminar boundary layer was performed on an equal spaced grid in n direction.For the transitional and turbulent region a nonuniform grid was applied for the n-direction.500 boundary layer locations were defined along the suction surface s-direction, all equal spaced .In the laminar region the boundary layer profile was resolved by 100 points in n-direction.The turbulent boundary layer profile was calculated by using 300 points in n-direction.In the turbulent part the grid was stretched in n-directions according to see Figure 5 g Δn Δn − , g 1.02, . . ., 1.06.

3.8
The height h of the grid maximum n-coordinate was chosen as h 0.02 • l for the laminar region and as h 0.05 • l for the remaining regions transitional and turbulent boundary layer .The number of grid points used for the inverse boundary layer calculation is large or even very large .But by doing so the numerical mistake is kept very small.The height of the grid is chosen large enough that the outer edge of the boundary layer lies within the grid.Although this large number of grid points is applied the method works quickly on PCs maximum CPU-time is about 30 seconds .
Goto et al. 4 apply an integral method.Integral methods work faster than solution methods based on finite difference schemes.Today PCs work so quickly that the difference in CPU time between both methods is not recognized by the user.
Thwaites's integral method for calculating laminar boundary layers 5 contains a simple correlation between the shape factor H and a parameter m Goto et al. 4 apply Thwaites's integral method inversely : Thwaites's correlation is used in order to prescribe the constant shape factor H in the laminar region.m is a measure of the deceleration of the outer velocity C s .If m is larger than 0.1 the laminar boundary layer separates.

Transitional Boundary Layer
The onset of turbulence is calculated by Walker's criterion: The length of the transitional region Δs trans is assumed to be 15% of the total length of the suction surface.
In the transitional region the turbulent viscosity is modeled by the intermittency function of Dhawan 12 .The intermittency function γ s is 3.12

Turbulence Model
A two-zone turbulence model is used in order to model the turbulent viscosity μ t 7, 13 .It belongs to the group of algebraic turbulence models.μ t is as follows:

ISRN Applied Mathematics
l stands for the mixing length.The boundary layer is considered to exist of an inner and an outer part.For the inner part l is calculated as follows: u τ stands for the friction velocity u τ τ W /ρ. τ W is the shear stress on the wall.In the outer region the turbulent viscosity is calculated by

3.15
The switch from the inner to the outer model is performed at the smallest n-location, where μ t inner ≥ μ t outer .

Potential Flow Model
Schwering's method 3 is applied in order to calculate the geometry of the cascade.Schwering's method works with the application of a resulting stream function consisting of the superposition of the stream function of a translational and circulation flow.The resulting stream function ψ R is see Figure 6 ψ R x, y ψ trans x, y ψ circ x, y , 3.16

3.18
e is a coordinate running around the profile in an anticlockwise direction.It starts at the trailing edge and ends there with the value L. x e and y e are the coordinates of the profile.
x and y stand for the coordinates of an arbitrary point in the flow field.γ is the strength of the vorticity per unit length.Schwering transforms the coordinate e of 3.16 -3.18 on the interval 0 2• π and the points with the coordinates x and y are put on the contour of the profiles.Dimensionless coordinates are introduced.They are

3.19
By using the dimensionless values and the relation tan 3.20 the resulting stream function 3.16 can be written as follows:

3.21
σ and τ start and end at the trailing edge.They run in anti-clockwise direction around the profile.K σ, τ and γ are see 3

3.22
Equation 3.21 is solved iteratively.At the beginning an arbitrary uncambered profile is taken i.e., x τ and y τ are estimated with a stagger angle

3.23
γ τ is known from the inverse calculation of the boundary layer and the calculation of the velocity distribution on the pressure side, respectively.One gets a better approximation of the geometry of the cascade i.e., a better approximation of η when the right-hand side of 3.21 is calculated with the estimated coordinates ξ τ , η τ , and γ τ .The new η-values are taken for a further calculation of the right-hand side of 3.21 in order to improve the approximation of the geometry of the cascade in such a way that the prescribed velocity distribution is realized.This process is repeated until the changes of η fall under a prescribed limit.As a rule five up to ten iterations are needed in order to get a converged solution.
After each iteration the spacing t must be recalculated by

3.24
As a rule the so calculated spacing ratio t/l does not differ significantly from the prescribed value.

Inverse Design
A cascade turning the flow from β 1 40 o to β 2 0 o was designed by the method described previously.The chosen design value are x stag l 0.003.

4.1
The boundary layer values are depicted in Figure 7. c f stands for the friction coefficient.It is

4.2
In the leading edge region the friction coefficient is large and becomes smaller downstream.
Before laminar separation appears the transition from laminar to turbulent onsets and the friction coefficient increases.In the turbulent region the friction coefficient decreases again, but it does not reach zero.In Figure 8 the geometry of the inversely designed cascade is shown.Three according velocity distributions are depicted in Figure 9.

Discussion of the Result
The velocity distribution with the label "Velocity Prescribed" is the result of the inverse boundary layer calculation and the calculation of the velocity distribution for the pressure side, respectively.After the cascade had been designed inversely, the velocity distribution was calculated by the application of Martensen method 14, 15 and a Navier-Stokes solver FLUENT 16 .The prescribed velocity distribution is in a good agreement with the velocity distribution calculated by Martensen method potential theory see Figure 9 .In the trailing edge region no stagnation point is existing for the Martensen solution, because the velocity distribution was modified by a simple trailing edge model described by Gostelow 17 .Gostelow explains a simple "first viscous approximation" by applying a faring-in of the velocity distribution.From the location x/l 0.85 to the trailing edge the velocity is not calculated by potential theory on the suction surface.Instead C is extrapolated from x/l 0.85 to x/l 1.0.For the pressure side C is approximated by a polynomial of second order running from x/l 0.85 to the trailing edge in such a way that C has no gap at x/l 0.85 and x/l 1.0.Additionally, there is no gap in the slope of C at x/l 0.85.
There is a difference between the velocity distribution calculated by potential theory and the Navier-Stokes solver, especially in the trailing edge region.The Navier-Stokes solution indicates that the flow moves around the trailing edge from the suction to the pressure side.This causes an influence on the velocity distribution upstream.Nevertheless, there is a good similarity between the Navier-Stokes solution and the solution based on potential theory.
The grid of the Navier-Stokes-calculation is shown in Figure 10.In the wall region the grid was refined in such a way that the maximum y -value was y 7.1 exception is the region of the stagnation point .The turbulence was modeled by the Spalart-Allmaras turbulence model 18 .On the inflow side not shown in Figure 10, because it is located upstream with the distance t from the left side of the shown grid the velocity components C 1,x and C 1,y were set.The static pressure was calculated by FLUENT.On the outflow side not shown in Figure 10, because it is located downstream with the distance 3•t from the right side of the shown grid a constant static pressure was imposed.Periodic boundary conditions were imposed on the upper and lower side of the calculation domain.A second-order upwind scheme was applied for the numerical calculations.The iterations were stopped after the residuum of all equations continuity, momentum, turbulence had fallen at least five order of magnitudes.
Navier-Stokes calculations were carried out for the inversely designed cascade and for a NACA65-cascade.Both cascades were designed to turn the flow from β 1 40 o to β 2 0 o .The Navier-Stokes calculations were performed for different flow angles β 1 see Figures 11 and  12 .The NACA65-cascade was designed by the data described in 1 .Although the turning Δβ β 1 −β 2 is smaller for the NACA65-cascade than for the inversely designed cascade for all   flow angles β 1 see Figure 12 , the flow losses are larger for the NACA65-cascade than for the inversely designed one see Figure 11 .On the suction side of the inversely designed profiles a concave velocity distribution exists; that is, at the beginning of the turbulent boundary layer the outer velocity is decelerated stronger than in the trailing edge region.Stratford 19 has shown that this kind of deceleration is advantageous concerning the friction losses inside the boundary layer.

Summary
An inverse design method for cascades was described consisting of an inverse boundary layer calculation and an inverse application of potential theory.The inverse boundary layer calculation is restricted to the suction surface of the profiles.On the pressure side the velocity distribution was prescribed in such a way that the flow is turned from the flow angle β 1  to flow angle β 2 .The method was applied for designing a cascade for β 1 40 o to β 2 0 o for a spacing ratio t/l 0.7166.Navier-Stokes simulations show that the performance of the inversely designed cascade possess a better loss performance than an according cascade consisting of NACA65-profiles.
The described method is based on an idea of Goto et al. 4 .It was modified in such a way that the inverse boundary layer method can be applied with different models, especially turbulence models.Secondly, a new procedure for prescribing the velocity distribution on the pressure side was introduced.Navier-Stokes simulations show that the prescribed velocity distribution for the design point can be realized.A concave velocity distribution is calculated on the suction surface in the turbulent region i.e., at the beginning the deceleration is strong and becomes smaller further downstream .This kind of deceleration is advantageous in order to reduce friction losses 19 .

Figure 3 :
Figure 3: Shape factor in the region of separation point.

Figure 4 :
Figure 4: Velocity distribution on the pressure side.

Nomenclature A :
Constant C: Velocity at the outer edge of the boundary layer C m : Maximum outer velocity on the suction surface C 1 : Velocity in front of the cascade C 2 : Velocity behind the cascade H: Shape factor l: Chord length n: Coordinate normal to s-direction Re: Reynolds number p: Static pressure p tot : Total pressure s: Coordinate along suction surface s m : Coordinate along suction surface for C m t: Spacing u: Velocity component of boundary layer in s-direction v : Velocity component of boundary layer normal to s-direction x: x-Coordinate of coordinate system of profile y: y-Coordinate of coordinate system of profile x: x-Coordinate of coordinate system of cascade y: y-Coordinate of coordinate system of cascade α: Constant of turbulence model β: Flow angle Γ: Circulation γ: Strength of vorticity or intermittency function δ: Boundary layer thickness δ * : Displacement thickness η: Dimensionless coordinate of y cascade coordinate system ξ: Dimensionless coordinate of x cascade coordinate system κ: Constant μ: Dynamic viscosity μ t : Turbulent viscosity ν: Kinematic viscosity ρ: Density θ: Momentum thickness ψ: Value of stream function ζ V 1 : Dimensionless pressure loss coefficient Re θ,onset the onset of turbulence occurs.H m and Re θ,m stand for the averaged values of H and Re θ along the laminar boundary layer, respectively.Re θ,inst is the according Reynolds number for the instability point.It is calculated byDrela et al. 11 : .7 − 0.32 • H m • Re θ,m Re θ,inst .3.10 If Re θ C s • θ/ν >