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

Authors: Martin Böhle
Journal: ISRN Applied Mathematics
 (2012)

---

## Abstract

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.


---

## Body

## 1. 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 HTE=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.


## 2. 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.Figure 1
Nomenclature of the cascade.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](2.1)δ*(s)=∫0∞(1-u(n)C(s))⋅dn,(2.2)Θ(s)=∫0∞(1-u(n)C(s))⋅u(n)C(s)⋅dn,(2.3)H(s)=δ*(s)Θ(s).
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(2.4)C(s)=(Cm-CB)⋅1-(1-ssm)5+CB,0≤s≤sm,
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 CB. The velocity CB is explained in more detail when the velocity distribution on the pressure side is introduced.Figure 2
Nomenclature of the cascade.In the second region, the laminar flow region, the shape factorH 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 factorH(s). Again an idea of Goto is applied to describe the shape factor by a half cosine wave as follows:(2.5)H(s)=Hlam-Hlam-Hturb2⋅(1+cos⁡[π⋅(1+s-sonsetsturb-sonset)]),sonset<s≤sturb.Hlam and Hturb 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 factorH(s) is prescribed by a quarter sine wave as follows:(2.6)H(s)=Hturb+(HTE-Hturb)⋅sin⁡(π2⋅s-sturbsTE-sturb),sturb<s≤sTE.
Index TE stands for the according values at the trailing edge. Goto recommends that HTE should be chosen as HTE=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.Figure 3
Shape factor in the region of separation point.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(2.7)C(s)=CB⋅(1-1-(1-xxstag)a),0≤x≤xstag.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 xstag. By this coordinate the location of the stagnation point can be prescribed (see Figure 4). The circulation Γ is calculated by(2.8)Γ=C1⋅t⋅(tanβ1-tanβ2)⋅cos⁡β1.
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(2.9)Γ=∮rC(r)⋅dr
must be fulfilled. r stands for a coordinate running in a clockwise direction around the profile.Figure 4
Velocity distribution on the pressure side.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 coordinates 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.


## 3. Models

### 3.1. Boundary Layer Model

The boundary layer flow is modeled by solving the boundary layer equations:(3.1)∂u∂s+∂v∂n=0,ρ⋅(u⋅∂u∂s+v⋅∂v∂n)=ρ⋅C⋅dCds+(μ+γ(s)⋅μt)⋅∂2u∂s2.μ 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)u=0,v=0,forn=0u=C,forn⟶∞.In the case of inverse boundary layer calculation the boundary conditions are formulated as follows:(3.3)u=0,v=0,forn=0u=C,H(s)=Hprescribed.


#### 3.1.1. 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 Figure5)(3.4)ρ⋅(ui,j⋅ui+1,j-ui,jΔs+vi,j⋅ui+1,j+1-ui+1,j-1Δn++Δn-)=ρ⋅Ci⋅Ci+1-CiΔs+2Δn++Δn-⋅((μ+μt)i,j+1/2⋅ui+1,j+1-ui+1,jΔn+)-2Δn++Δn-⋅((μ+μt)i,j-1/2⋅ui+1,j-ui+1,j-1Δn-).
The boundary conditions are(3.5)ui,j=0=0,vi,j=0=0,forn=0ui,j=max⁡=CFigure 5
Nomenclature of finite difference.and for inverse application(3.6)ui,j=0=0,vi,j=0=0,forn=0H(s)=Hprescribed.At locationi 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 ui+1,j(j=0,…,nmax-1). By knowing ui+1,j(j=0,…,nmax-1) the velocity components vi+1,j(j=0,…,nmax-1) are calculated by the continuity equation by taking into account that vi+1,0=0. The discretisation of the continuity equation is(3.7)ui+1,j+1+ui+1,j-ui,j+1-ui,j2⋅Δs+vi+1,j+1-vi+1,jΔn+=0.
The inverse boundary layer calculation is performed as described below(i)
At location with indexi the velocity profile is known. It fulfills the boundary condition ui,j=0(n=0)=0, vi,j=0(n=0)=0, and H(si)=Hprescribed.(ii)
Ci+1 is estimated. By the application of the finite difference scheme (3.4) the velocities u(si+1) are calculated. The boundary conditions ui+1,j=0(n=0)=0 and C=Ci+1 are taken into account.(iii)
The velocitiesvi+1,j are calculated by (3.7). The boundary condition vi+1,j=0(n=0)=0 is taken into account.(iv)
Hi+1 is calculated (see (2.3)).(v)
If|(Hi+1-Hprescribed)/Hprescribed|>ϵ the estimated value Ci+1 is corrected in order to reduce the difference Hi+1-Hprescribed. So the procedure starts at the second step again and is repeated until |(Hi+1-Hprescribed)/Hprescribed|<ϵ. ϵ 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 inn 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)(3.8)g=Δn+Δn-,g=1.02,…,1.06.The heighth 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):(3.9)H(m)=0.07310.14-m+2.088,m=-Θ2ν⋅dCds,0≤m≤0.1.
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.


#### 3.1.2. Transitional Boundary Layer

The onset of turbulence is calculated by Walker’s criterion:(3.10)ReΘ,onset=(1.7-0.32⋅Hm)⋅Reθ,m+Reθ,inst.
If Reθ=C(s)·θ/ν>Reθ,onset the onset of turbulence occurs. Hm 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 by Drela et al. [11]:(3.11)log⁡10(Reθ,inst)=(1.415H-1-0.489)⋅tanh(20H-1-12.9)+3.295H-1+0.44.
The length of the transitional region Δstrans 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)γ(s)=1-e-0.412⋅(s-sonset/λ)2,sonset≤s≤sturb,λ2=-0.412ln⁡(0.01)⋅Δstrans⁡2.


#### 3.1.3. 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:(3.13)μt=ρ⋅l2|∂u∂n|.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:(3.14)l=κ⋅n(1-e-n+/A),κ=0.4.,A=26,n+=n⋅uτν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)μt=α⋅ρ⋅C(s)⋅δ*α=0.0168⋅1.551+0.55⋅(1-e(-0.234⋅z-0.298⋅z))z=ReΘ425-1.
The switch from the inner to the outer model is performed at the smallest n-location, where (μt)inner≥(μt)outer.


### 3.2. 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)(3.16)ψR(x¯,y¯)=ψtrans⁡(x¯,y¯)+ψcirc(x¯,y¯),(3.17)ψtrans⁡=C∞,x⋅y¯-C∞,y⋅x¯,(3.18)ψcirc=-12⋅π⋅∫e′=0e′=Lγ(e′)⋅ln⁡(Π)⋅de′,Π=sinh⁡2[πt⋅(x¯-x¯(e′))]+sin⁡2[πt⋅(y¯-y¯(e′))].
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.Figure 6
Nomenclature of Schwering’s method [3].Schwering transforms the coordinatee′ of (3.16)–(3.18) on the interval [02·π] 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)ξ=2⋅π⋅x¯t,η=2⋅π⋅y¯t,γ*=γC∞,x,ρ¯=(dξdτ)2+(dηdτ)2,ψ¯R=ψRC∞,x⋅t,γ¯=γ*⋅ρ¯,0≤σ≤2⋅π,0≤τ≤2⋅π.
By using the dimensionless values and the relation(3.20)tanβ∞=tanβ1+14⋅π∫τ=02⋅πγ¯(τ)⋅dτ,the resulting stream function(3.16) can be written as follows:(3.21)η(σ)=ψ¯R+tanβ1+12⋅π∫τ=02⋅πγ¯(τ)⋅Π⋅dτΠ=K(σ,τ)+12⋅ξ(σ).σ and τ start and end at the trailing edge. They run in anti-clockwise direction around the profile. K(σ,τ) and γ¯ are (see [3])(3.22)K(σ,τ)=12⋅ln⁡(cosh⁡[ξ(σ)-ξ(τ)]-cos⁡[η(σ)-η(τ)]2),γ¯=γ(τ)C∞,x⋅(dξdτ)2+(dηdτ)2.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)λ=arctan(12⋅(tanβ1+tanβ2)).γ(τ) 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 spacingt must be recalculated by(3.24)t=∮rC(r)⋅drC1⋅cos⁡β1⋅(tanβ2-tanβ1).As a rule the so calculated spacing ratiot/l does not differ significantly from the prescribed value.


## 3.1. Boundary Layer Model

The boundary layer flow is modeled by solving the boundary layer equations:(3.1)∂u∂s+∂v∂n=0,ρ⋅(u⋅∂u∂s+v⋅∂v∂n)=ρ⋅C⋅dCds+(μ+γ(s)⋅μt)⋅∂2u∂s2.μ 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)u=0,v=0,forn=0u=C,forn⟶∞.In the case of inverse boundary layer calculation the boundary conditions are formulated as follows:(3.3)u=0,v=0,forn=0u=C,H(s)=Hprescribed.


### 3.1.1. 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 Figure5)(3.4)ρ⋅(ui,j⋅ui+1,j-ui,jΔs+vi,j⋅ui+1,j+1-ui+1,j-1Δn++Δn-)=ρ⋅Ci⋅Ci+1-CiΔs+2Δn++Δn-⋅((μ+μt)i,j+1/2⋅ui+1,j+1-ui+1,jΔn+)-2Δn++Δn-⋅((μ+μt)i,j-1/2⋅ui+1,j-ui+1,j-1Δn-).
The boundary conditions are(3.5)ui,j=0=0,vi,j=0=0,forn=0ui,j=max⁡=CFigure 5
Nomenclature of finite difference.and for inverse application(3.6)ui,j=0=0,vi,j=0=0,forn=0H(s)=Hprescribed.At locationi 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 ui+1,j(j=0,…,nmax-1). By knowing ui+1,j(j=0,…,nmax-1) the velocity components vi+1,j(j=0,…,nmax-1) are calculated by the continuity equation by taking into account that vi+1,0=0. The discretisation of the continuity equation is(3.7)ui+1,j+1+ui+1,j-ui,j+1-ui,j2⋅Δs+vi+1,j+1-vi+1,jΔn+=0.
The inverse boundary layer calculation is performed as described below(i)
At location with indexi the velocity profile is known. It fulfills the boundary condition ui,j=0(n=0)=0, vi,j=0(n=0)=0, and H(si)=Hprescribed.(ii)
Ci+1 is estimated. By the application of the finite difference scheme (3.4) the velocities u(si+1) are calculated. The boundary conditions ui+1,j=0(n=0)=0 and C=Ci+1 are taken into account.(iii)
The velocitiesvi+1,j are calculated by (3.7). The boundary condition vi+1,j=0(n=0)=0 is taken into account.(iv)
Hi+1 is calculated (see (2.3)).(v)
If|(Hi+1-Hprescribed)/Hprescribed|>ϵ the estimated value Ci+1 is corrected in order to reduce the difference Hi+1-Hprescribed. So the procedure starts at the second step again and is repeated until |(Hi+1-Hprescribed)/Hprescribed|<ϵ. ϵ 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 inn 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)(3.8)g=Δn+Δn-,g=1.02,…,1.06.The heighth 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):(3.9)H(m)=0.07310.14-m+2.088,m=-Θ2ν⋅dCds,0≤m≤0.1.
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.


### 3.1.2. Transitional Boundary Layer

The onset of turbulence is calculated by Walker’s criterion:(3.10)ReΘ,onset=(1.7-0.32⋅Hm)⋅Reθ,m+Reθ,inst.
If Reθ=C(s)·θ/ν>Reθ,onset the onset of turbulence occurs. Hm 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 by Drela et al. [11]:(3.11)log⁡10(Reθ,inst)=(1.415H-1-0.489)⋅tanh(20H-1-12.9)+3.295H-1+0.44.
The length of the transitional region Δstrans 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)γ(s)=1-e-0.412⋅(s-sonset/λ)2,sonset≤s≤sturb,λ2=-0.412ln⁡(0.01)⋅Δstrans⁡2.


### 3.1.3. 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:(3.13)μt=ρ⋅l2|∂u∂n|.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:(3.14)l=κ⋅n(1-e-n+/A),κ=0.4.,A=26,n+=n⋅uτν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)μt=α⋅ρ⋅C(s)⋅δ*α=0.0168⋅1.551+0.55⋅(1-e(-0.234⋅z-0.298⋅z))z=ReΘ425-1.
The switch from the inner to the outer model is performed at the smallest n-location, where (μt)inner≥(μt)outer.


## 3.1.1. 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 Figure5)(3.4)ρ⋅(ui,j⋅ui+1,j-ui,jΔs+vi,j⋅ui+1,j+1-ui+1,j-1Δn++Δn-)=ρ⋅Ci⋅Ci+1-CiΔs+2Δn++Δn-⋅((μ+μt)i,j+1/2⋅ui+1,j+1-ui+1,jΔn+)-2Δn++Δn-⋅((μ+μt)i,j-1/2⋅ui+1,j-ui+1,j-1Δn-).
The boundary conditions are(3.5)ui,j=0=0,vi,j=0=0,forn=0ui,j=max⁡=CFigure 5
Nomenclature of finite difference.and for inverse application(3.6)ui,j=0=0,vi,j=0=0,forn=0H(s)=Hprescribed.At locationi 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 ui+1,j(j=0,…,nmax-1). By knowing ui+1,j(j=0,…,nmax-1) the velocity components vi+1,j(j=0,…,nmax-1) are calculated by the continuity equation by taking into account that vi+1,0=0. The discretisation of the continuity equation is(3.7)ui+1,j+1+ui+1,j-ui,j+1-ui,j2⋅Δs+vi+1,j+1-vi+1,jΔn+=0.
The inverse boundary layer calculation is performed as described below(i)
At location with indexi the velocity profile is known. It fulfills the boundary condition ui,j=0(n=0)=0, vi,j=0(n=0)=0, and H(si)=Hprescribed.(ii)
Ci+1 is estimated. By the application of the finite difference scheme (3.4) the velocities u(si+1) are calculated. The boundary conditions ui+1,j=0(n=0)=0 and C=Ci+1 are taken into account.(iii)
The velocitiesvi+1,j are calculated by (3.7). The boundary condition vi+1,j=0(n=0)=0 is taken into account.(iv)
Hi+1 is calculated (see (2.3)).(v)
If|(Hi+1-Hprescribed)/Hprescribed|>ϵ the estimated value Ci+1 is corrected in order to reduce the difference Hi+1-Hprescribed. So the procedure starts at the second step again and is repeated until |(Hi+1-Hprescribed)/Hprescribed|<ϵ. ϵ 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 inn 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)(3.8)g=Δn+Δn-,g=1.02,…,1.06.The heighth 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):(3.9)H(m)=0.07310.14-m+2.088,m=-Θ2ν⋅dCds,0≤m≤0.1.
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.


## 3.1.2. Transitional Boundary Layer

The onset of turbulence is calculated by Walker’s criterion:(3.10)ReΘ,onset=(1.7-0.32⋅Hm)⋅Reθ,m+Reθ,inst.
If Reθ=C(s)·θ/ν>Reθ,onset the onset of turbulence occurs. Hm 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 by Drela et al. [11]:(3.11)log⁡10(Reθ,inst)=(1.415H-1-0.489)⋅tanh(20H-1-12.9)+3.295H-1+0.44.
The length of the transitional region Δstrans 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)γ(s)=1-e-0.412⋅(s-sonset/λ)2,sonset≤s≤sturb,λ2=-0.412ln⁡(0.01)⋅Δstrans⁡2.


## 3.1.3. 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:(3.13)μt=ρ⋅l2|∂u∂n|.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:(3.14)l=κ⋅n(1-e-n+/A),κ=0.4.,A=26,n+=n⋅uτν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)μt=α⋅ρ⋅C(s)⋅δ*α=0.0168⋅1.551+0.55⋅(1-e(-0.234⋅z-0.298⋅z))z=ReΘ425-1.
The switch from the inner to the outer model is performed at the smallest n-location, where (μt)inner≥(μt)outer.


## 3.2. 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)(3.16)ψR(x¯,y¯)=ψtrans⁡(x¯,y¯)+ψcirc(x¯,y¯),(3.17)ψtrans⁡=C∞,x⋅y¯-C∞,y⋅x¯,(3.18)ψcirc=-12⋅π⋅∫e′=0e′=Lγ(e′)⋅ln⁡(Π)⋅de′,Π=sinh⁡2[πt⋅(x¯-x¯(e′))]+sin⁡2[πt⋅(y¯-y¯(e′))].
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.Figure 6
Nomenclature of Schwering’s method [3].Schwering transforms the coordinatee′ of (3.16)–(3.18) on the interval [02·π] 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)ξ=2⋅π⋅x¯t,η=2⋅π⋅y¯t,γ*=γC∞,x,ρ¯=(dξdτ)2+(dηdτ)2,ψ¯R=ψRC∞,x⋅t,γ¯=γ*⋅ρ¯,0≤σ≤2⋅π,0≤τ≤2⋅π.
By using the dimensionless values and the relation(3.20)tanβ∞=tanβ1+14⋅π∫τ=02⋅πγ¯(τ)⋅dτ,the resulting stream function(3.16) can be written as follows:(3.21)η(σ)=ψ¯R+tanβ1+12⋅π∫τ=02⋅πγ¯(τ)⋅Π⋅dτΠ=K(σ,τ)+12⋅ξ(σ).σ and τ start and end at the trailing edge. They run in anti-clockwise direction around the profile. K(σ,τ) and γ¯ are (see [3])(3.22)K(σ,τ)=12⋅ln⁡(cosh⁡[ξ(σ)-ξ(τ)]-cos⁡[η(σ)-η(τ)]2),γ¯=γ(τ)C∞,x⋅(dξdτ)2+(dηdτ)2.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)λ=arctan(12⋅(tanβ1+tanβ2)).γ(τ) 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 spacingt must be recalculated by(3.24)t=∮rC(r)⋅drC1⋅cos⁡β1⋅(tanβ2-tanβ1).As a rule the so calculated spacing ratiot/l does not differ significantly from the prescribed value.


## 4. Results

### 4.1. Inverse Design

A cascade turning the flow fromβ1=40o to β2=0o was designed by the method described previously. The chosen design value are(4.1)tl=0.8,Re=ρ⋅C1⋅lμ=6.5⋅105,sMl=0.07,CMC1=1.4,m=0.023,Hturb=1.4,HTE=1.7,x1,pressl=0.1,CBC1=0.82,x2,pressl=0.4,xstagl=0.003.The boundary layer values are depicted in Figure7. cf stands for the friction coefficient. It is(4.2)cf=τw(ρ/2)⋅C12.Figure 7
Values of inverse boundary layer calculation.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 Figure8 the geometry of the inversely designed cascade is shown. Three according velocity distributions are depicted in Figure 9.Figure 8
Geometry of the inverse designed cascade.Figure 9
Velocity distribution of inversely designed cascade.


### 4.2. 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 Figure10. 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 C1,x and C1,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.Figure 10
Numerical grid for Navier-Stokes calculation.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=40o to β2=0o.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.Figure 11
Dimensionless loss coefficient in dependence ofβ1.Figure 12
Δβ in dependence of β1.


## 4.1. Inverse Design

A cascade turning the flow fromβ1=40o to β2=0o was designed by the method described previously. The chosen design value are(4.1)tl=0.8,Re=ρ⋅C1⋅lμ=6.5⋅105,sMl=0.07,CMC1=1.4,m=0.023,Hturb=1.4,HTE=1.7,x1,pressl=0.1,CBC1=0.82,x2,pressl=0.4,xstagl=0.003.The boundary layer values are depicted in Figure7. cf stands for the friction coefficient. It is(4.2)cf=τw(ρ/2)⋅C12.Figure 7
Values of inverse boundary layer calculation.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 Figure8 the geometry of the inversely designed cascade is shown. Three according velocity distributions are depicted in Figure 9.Figure 8
Geometry of the inverse designed cascade.Figure 9
Velocity distribution of inversely designed cascade.


## 4.2. 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 Figure10. 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 C1,x and C1,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.Figure 10
Numerical grid for Navier-Stokes calculation.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=40o to β2=0o.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.Figure 11
Dimensionless loss coefficient in dependence ofβ1.Figure 12
Δβ in dependence of β1.


## 5. 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=40o to β2=0o 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].



---

*Source: 148607-2012-06-05.xml*