CASCADE-XPERT: SOFTWARE DEVELOPMENT FOR CASCADE ...

7 downloads 445 Views 2MB Size Report
Timisoara, Romania, June 10 - 11, 2005. CASCADE-XPERT: SOFTWARE DEVELOPMENT. FOR CASCADE HYDRODYNAMICS. Romeo F. SUSAN-RESIGA*, ...
Scientific Bulletin of the Politehnica University of Timisoara Transactions on Mechanics Special issue

Workshop on Vortex Dominated Flows – Achievements and Open Problems Timisoara, Romania, June 10 - 11, 2005

CASCADE-XPERT: SOFTWARE DEVELOPMENT FOR CASCADE HYDRODYNAMICS Romeo F. SUSAN-RESIGA*, Professor Department of Hydraulic Machinery “Politehnica” University of Timişoara

Teodora FRUNZĂ, PhD student Department of Hydraulic Machinery “Politehnica” University of Timişoara

Sebastian MUNTEAN, Senior Researcher* Center of Advanced Research in Engineering Sciences Romanian Academy - Timişoara Branch

Sandor BERNAD, Senior Researcher Center of Advanced Research in Engineering Sciences Romanian Academy - Timişoara Branch

Cristian ARMEANA, software engineer Vital Soft Solutions, Timişoara *Corresponding author: Bv Mihai Viteazu 1, 300222, Timisoara, Romania Tel.: (+40) 256 403692, Fax: (+40) 256 403700, Email: [email protected] ABSTRACT The paper presents the current status of the CASCADExpert software developed by the authors for hydrofoil cascade analysis, design and optimization. We first briefly review the problem formulation for the inviscid, incompressible and irrotational flow in 2D linear cascades, and the Finite Element Method for solving the streamfunction boundary-value problem. The boundary conditions for the viscous-inviscid coupling are developed to include the boundary layer displacement thickness. Second, we present the shape parametric representation using a set of special orthonormal function, for which we give the explicit analytical expressions. The actual hydrofoil shape used in computations is obtained via a leastsquares representation outlined in the paper. Third, we present the procedure for averaging the velocity field in order to extract the information for a Q3D turbomachinery computation. The paper is concluded with a review of the present capabilities of our CASCADExpert code. KEYWORDS Cascade hydrodynamics, streamfunction formulation, finite element, viscous-inviscid interaction, modal representation of hydrofoil shape. NOMENCLATURE G V ,V x ,V y ,Vn ,Vτ ,Vtr ,Ve

[-] dimensionless velocity

vector and velocity components, transpiration velocity, boundary layer edge velocity

ψ

[-]

streamfunction



δ [-] boundary layer displacement thickness α in , α out , α s [-] inlet, outlet, and stagger angles p − pref [-] pressure coefficient cp = 2 /2 ρVref Subscripts and Superscripts x, y directions normal and parallel with cascade frontline n,τ normal and tangential directions at a point on the boundary boundary layer edge e cascade inflow section in out cascade outflow section ABBREVIATIONS ref reference section 1. INTRODUCTION Analysis of the flow in hydrofoil cascades is one of the main steps in hydraulic machines design and optimization. Both analytical and numerical methods have been developed for more than a century to accurately compute the cascade flows, [12], and extensive experimental investigations have been performed. However, the modern requirements for hydraulic turbines and pumps design can be met only by using advanced numerical simulation techniques, which are robust and fast to allow iterative geometry optimization. Our current efforts are aimed at developing such a modern numerical tool, therefore

Proceedings of the Workshop on VORTEX DOMINATED FLOWS. ACHIEVEMENTS AND OPEN PROBLEMS, Timisoara, Romania, June 10-11, 2005

60

called CASCADExpert, which incorporates the latest achievements in flow modeling, numerical algorithms, computer implementation, graphics user interface. This paper reviews the current status of our code development, and presents the models we have chosen for implementing the viscous-inviscid interaction. 2. MATHEMATICAL MODEL FOR CASCADE HYDRODYNAMICS The inviscid, incompressible and irrotational flow model employed by our software CASCADExpert is based on the streamfunction formulation [19]. For incompressible flows the divergence of the velocity vanishes everywhere. Vectors satisfying this condition G are called solenoidal (or divergence free), i.e. ∇ ⋅ V = 0 . Any solenoidal vector can be written as the curl of G G another vector, called the vector potential, V = ∇ × β . G G For 2D plane flow, the vector potential is β = ψ e z , G where e z is the unit vector in the direction normal to the plane of flow. The velocity components in the x - and y -directions are ∂ψ ∂ψ , Vy = − . (1) ∂y ∂x Since the gradient of ψ is normal to the velocity G vector, V ⋅ ∇ψ = 0 , along a streamline we have ψ = const . , thus the name streamfunction. The difference in value of streamfunction between any two streamlines is equal to the volumetric flow rate passing between the two streamlines per unit width normal to the plane of motion. A flow that is both irrotational and incompressible G G must satisfy both the condition of zero curl ∇ × V = 0 and zero divergence. The second is automatically satG G isfied by introducing the streamfunction V = ∇ × (ψ e z ) , and enforcing the irrotationality condition we obtain the 2D Laplace equation,

front line; the inlet/outlet boundaries are aligned with the cascade front line. Note that the shape of the periodic boundary segments is arbitrary, the only restriction being the constant width of the strip. • Solid boundary (hydrofoil surface), with the unit length chord-line and stagger angle α s (measured with respect to the cascade front line).

Vx =

∂ 2ψ ∂ 2ψ = 0. + ∂x 2 ∂y 2

(2)

It is useful to recall here that in an irrotational flow of an incompressible fluid in a bounded region, the maximum velocity magnitude and the minimum pressure must occur on the boundary of this region. In our case, the maximum velocity will always occur on the hydrofoil surface. When analyzing the flow in a cascade, one can take into account the geometrical periodicity. As a result, the computational domain will correspond to a periodic strip as shown in Fig. 1. This computational domain in double connected, with the following boundary parts: • Fluid boundary (inlet, outlet, periodic); the periodic boundary segments are obtained by translation with the cascade pitch in the direction of the cascade

Figure 1. Domain of analysis for cascade hydrodynamics. On the cascade inlet section we impose uniform velocity distribution, at angle α in with respect to the cascade front line. The upstream axial velocity component is taken as reference velocity (unit velocity), Vin , x = 1 . As a result, the dimensionless upstream tangential velocity is given by the inflow angle, Vin , y = − cot α in . In terms of the streamfunction, on the inlet we impose a Neumann condition, ∂ψ = − cot α in . (3) ∂n At outlet we consider that the flow downstream the cascade reaches an uniform velocity distribution, with the direction given by the unknown angle α out , also measured with respect to the cascade front line. The downstream axial velocity component is the same as the upstream one, thanks to the continuity equation, Vout , x = 1 . The dimensionless downstream tangential velocity is given by the outflow angle, Vout , y = − cot α out . The corresponding boundary for the streamfunction is ∂ψ − cot α out = 0 , (4) ∂n where cot α out is on the left-hand-side as an additional unknown of the problem.

Proceedings of the Workshop on VORTEX DOMINATED FLOWS. ACHIEVEMENTS AND OPEN PROBLEMS, Timisoara, Romania, June 10-11, 2005

For two points on the upper and lower periodic boundaries, with the same abscissa, velocity vector must be the same, or for the velocity components we must have Vup , n = −Vdo , n and Vup , τ = −Vdo , τ . The periodicity of Vτ leads to periodic Neumann boundary condition  ∂ψ   ∂ψ    +  = 0.  ∂n  up  ∂n  do

(5)

On the other hand, since on inlet/outet sections we have specified only the direction of the flow, we must prescribe somehow the flow rate as well. This can be done via the Dirichlet-type condition for any pair of points on the periodic boundaries, (6) ψ up − ψ do = pitch , which also implies the periodicity of Vn ,  ∂ψ   ∂ψ   = 0.  +   ∂τ  up  ∂τ  do

(6’)

On the hydrofoil surface, the inviscid flow condition corresponds to the impenetrability, i.e. Vn = 0 . In other words, the hydrofoil is a streamline in inviscid flow, and without losing generality we can set (7) ψ = 0 on the hydrofoil. Since we have an additional unknown corresponding to the outflow angle, we need an additional condition to close the problem. This condition is derived from Joukowski’s hypothesis [17, §7.40]: the circulation in the case of a properly designed aerofoil always adjusts itself so that the velocity at the trailing edge is finite. This condition appears to be satisfied with reasonable exactness within the working range of well designed aerofoils. The physical explanation of the origin of the circulation is given by Milne-Thomson [17] as follows. When the motion is just starting, i.e. for low velocity, the flow is ordinary streamline flow, with the stagnation point just ahead of the trailing edge on the upper surface of the aerofoil. As the speed increases, even with small viscosity, the viscous forces increase and the air is no longer able to turn round the sharp edge and a vortex is formed. Since the circulation in any circuit large enough to enclose the aerofoil and the vortex was zero to start with it must remain zero, and so a circulation now exists round the aerofoil equal and opposite to that of the vortex. The vortex gets washed away downstream, and when the steady state is reached the circulation round the aerofoil remains. According to the theorem of Kutta (1902) and Joukowski (1906), an aerofoil at rest in a uniform G G wind of speed V∞ with circulation K round the G G aerofoil, undergoes a lift ρV∞ × K which obviously is perpendicular to the wind. The same theorem holds

61

G for cascades, where V∞ is the vector average of inlet and outlet velocities. Note that according this theorem, in inviscid flow the lift is independent of the form of the profile. The theory gives no drag since no account has been taken of the wake or viscosity. In practice, Joukowski’s condition is implemented by setting the trailing edge to be a stagnation point. Since the normal velocity was already set to zero on the whole hydrofoil surface, at trailing edge the tangential velocity must vanish too, ∂ψ = 0 at trailing edge. (8) ∂n Note that evaluating the normal derivative at an angular point (sharp trailing edge) might be difficult in theory, but this condition can be easily and accurately implemented in the finite element framework where the local flux, rather than the normal derivative, is computed. The problem (2) with conditions (3), (4), (5), (6), (7), and (8) is well defined, and allows the computation of the streamfunction as well as the direction of the outflow. Once the proble is solved, the axial and tangential hydrodynamic force components on the hydrofoil can be readily evaluated as: Vτ = −

(

)

Faxial pitch = cot 2 α out − cot 2 α in ρ 2 chord V x × chord × span 2 Ftang pitch = 2(cot α out − cot α in ) ρ 2 chord V x × chord × span 2 Note that these dimensionless force coefficients depend only on the inflow/outflow angles and the relative pitch. 2.1. Boundary conditions for viscous-inviscid interaction In viscous-inviscid interaction computations, where the irrotational flow is corrected to account for the boundary layer thickness, the displaced boundary of the inviscid region is the edge of the boundary layer. For small boundary layer thicknesses the displaced boundary of the computational region can be modeled by the introduction of the displacement thickness δ ∗ . In this case a mass balance over a boundary layer segment of infinitesimal length dl gives

(

)

d Veδ ∗ , (9) dl where Ve is the velocity at the edge of the boundary layer, and Vtr is the specific volumetric flow rate through the wall, or transpiration velocity [24]. The wall transpiration concept is currently used in wellknown codes such as XFOIL [7] to model the influence Vtr =

Proceedings of the Workshop on VORTEX DOMINATED FLOWS. ACHIEVEMENTS AND OPEN PROBLEMS, Timisoara, Romania, June 10-11, 2005

62

of the viscous layer on the inviscid flow. With the normal unit vector n pointing outside the fluid domain, and the tangent unit vector τ oriented such that the fluid domain remains on the left when walking the boundary, Figure 2, the local velocity components are ∂ψ ∂ψ and Vτ = − (10) ∂n ∂τ The transpiration velocity will always be Vtr = −Vn , but on the upper (suction) side of the hydrofoil we have dl = dτ , while on the lower (pressure) side dl = −dτ . As a result we have Vn =

(

)

(

)

( lo ) sides of the wake, we have the following two conditions: • Since the wake does not support transverse pressure gradients, the velocity on the upper side must equal the velocity on the lower side, for a pair of points Pup and Plo with x Pup = x Plo ∂ψ ∂n

+ up

∂ψ ∂n

=0

(13)

lo

This is the same condition as the one used for periodicity of tangential velocity component. • Since the wake has a displacement thickness, the corresponding transpiration condition is1

∂ ψ + Veδ ∗ = 0 on the upper side ∂τ   ψ −ψ up (11) ∂ψ  or − ∂ψ  = lo Vewake = (14) ∂ ∗ ∗ ψ − Veδ = 0 on the lower side ∂n lo  ∂n up  δ wake ∂τ This condition corresponds to a streamfunction jump By considering ψ = 0 at leading edge stagnation across the wake. Note that in the presence of the point (where Ve = 0 ), the above conditions become boundary layer one does not need any additional condition (like Kutta-Jukowski) to obtain a circulation ψ ± Veδ ∗ = 0 . On the other hand, the boundary layer around the hydrofoil. It is also obvious that at the edge velocity is on the suction side Ve = Vτ = −∂ψ / ∂n , trailing edge, the displacement thickness in the wake while on the pressure side we have Ve = −Vτ = ∂ψ / ∂n . is the sum of displacement thicknesses of the boundary With this observation, the transpiration condition that layers on the suction and pressure sides of the hydrofoil. models the viscous-inviscid interaction in incomA simple iterative calculation of δ ∗ interaction, pressible flow becomes according to Bradshaw et al. [3, §11.4], will proceed ∂ψ ∗ ∂ψ ψ as follows: ψ− δ = 0 or = , (12) ∂n δ ∗ ∂n 1. Compute the pressure distribution with δ ∗ = 0 , i.e. ψ = 0 on the hydrofoil surface; i.e. a homogeneous Cauchy condition on the hydrofoil surface. The boundary layer thickness δ ∗ is obtained by solving the boundary layer equations either in their differential or integral forms.

2. Calculate δ ∗ in boundary layer and wake using this pressure distribution 3. Recalculate pressure distribution on hydrofoil using the mixed boundary condition (12) for the streamfunction, and in the wake using (13) and (14). 4. If not converged, go to 2. However, an analysis of the boundary layer equations reveals the main difficulty of this direct approach, [5]. The steady Von Kármán equation and the entrainment equation can be combined to give dVe b dδ ∗ d + = (15) dl a dl a With coefficients a , b and d given in §2.4 of [5]. It is noted that the coefficient b corresponds to dH 1 / dH , with function H 1 (H ) having a minimum due to which the calculation of H can become impossible when H 1 is given. Using the closure relation for H − H 1 , b is zero at H = H sep , where

Figure 2. Normal and tangential unit vectors definition on the hydrofoil and its wake.

H sep = 4.029 for laminar flow, and H sep = 2.73 for

In the wake, for a pair of points with the same abscissa, but located on the upper ( up ) and lower

1 Of course, for the first iteration with inviscid flow ψ lo = ψ up = 0 since the wake is the streamline ψ

=0.

Proceedings of the Workshop on VORTEX DOMINATED FLOWS. ACHIEVEMENTS AND OPEN PROBLEMS, Timisoara, Romania, June 10-11, 2005

turbulent flow. The skin friction becomes zero at H = H sep which corresponds with the separation point. Furthermore, it can be shown that  > 0 for H < H sep (attached flow) b d > 0. = 0 for H = H sep (separation point) , a a  < 0 for H > H sep (separated flow) When treating equation (15) with the classical direct method, the velocity Ve is known from the external flow calculation and a simple ordinary differential equation for δ ∗ is to be solved: b dδ ∗ d dVe (16) = − a dl a dl At H = H sep the coefficient b / a vanishes and dδ * / dl becomes infinite unless the right-hand side of Eq. (16) vanishes as well,

d dVe − . (17) a dl For a converged solution of the viscous-inviscid problem condition (17) is valid. On the other hand, for an arbitrary prescribed Ve distribution the above relation is usually not satisfied. Moreover, even if (17) is satisfied at H = H sep , in (16) the derivative 0=

dδ * / dl remains undetermined, being of the form 0 / 0 . This problem was identified and explained by Goldstein in 1948, and it was coined in the literature under the name of “Goldstein singularity”. The viscous-inviscid interaction methods have an interesting history, with many subtle developments over the past two decades [23]. The main challenge in solving Prandtl’s boundary-layer equations has been to overcome the singularity at a point of steady flow separation. The current interactive boundarylayer models are based on a coupled problem that can be written as, in principle, two equations with two unknowns:

( ) = B (δ )

External inviscid flow: Ve = E δ ∗ Boundary-layer flow: Ve



(18a) (18b)

Here E denotes the external inviscid-flow operator, whereas B is the boundary-layer operator for prescribed displacement thickness; note that near flow separation the inverse B −1 does not exist. Equation (18b) is practically the boundary layer equation (15) written for a prescribed displacement thickness, ∗

dVe d b dδ (19) = − dl a a dl The term d / a remains positive when H = H sep , which implies that with this inverse scheme the calculations can be continued beyond the separation point.

63

In the classical, or “direct”, method Ve is computed from the inviscid-flow equation (18a), whereas the displacement thickness is determined from the viscous flow (18b), with a break-down of B −1 at separation. In late seventies a number of ideas have been put forward to circumvent the breakdown singularity. The simplest way is to invert the direction of the iterative process in the classical method. One obtains the so-called “inverse” method where the boundary layer is solved with prescribed displacement thickness. For engineering applications, however, the inverse method converges very slowly and it has not being used on a large scale. To speed up the convergence other methods were developed, of which two have survived: the semi-inverse method and the quasisimultaneous method. The semi-inverse method is a mixture of the direct and inverse method: it solves the boundary-layer equations with prescribed displacement thickness, and the inviscid flow in the traditional way (also with prescribed displacement thickness):

( = B (δ

) ) inverse; + ω (V − V ).

VeE = E δ ∗ ( n −1) direct; VeB

∗ ( n −1)

δ ∗ ( n ) = δ ∗ ( n −1)

B e

(20)

E e

In order to obtain convergence, some tuning of the relaxation parameter ω is required, and a fair convergence can be obtained. The quasi-simultaneous method reflects the lack of hierarchy between inviscid and viscous subdomains: in principle, it wants to solve both subdomain problems simultaneously. When the boundary later is modeled by an integral formulation a simultaneous coupling is well feasible, e.g. the XFOIL code [7], or [24]. However, when in both domains a field formulation is chosen, software complexity may prohibit a practical implementation. The idea of the quasi-simultaneous approach was to solve the boundary layer equations simultaneously with a simple but good approximation of the inviscid flow, which was termed the interaction law. The difference between this approximation and the “exact” inviscid flow can be handled iteratively. In this way, the quasi-simultaneous method can be formulated as (21) Ve( n ) − I δ ∗ ( n ) = E δ ∗ ( n −1) − I δ ∗ ( n −1) , Ve( n )

( ) ( − B (δ )= 0 , ∗(n)

) (

)

where the interaction law reads

( )

I δ∗ =

1 π

d δ ∗ dζ ∫Γ dζ τ − ζ ,

(22)

with τ the streamwise boundary-layer coordinate. One can observe that in (21) the interaction law is used in defect formulation, i.e. it does not influence

Proceedings of the Workshop on VORTEX DOMINATED FLOWS. ACHIEVEMENTS AND OPEN PROBLEMS, Timisoara, Romania, June 10-11, 2005

64

the final converged result, but it only enhances the rate of convergence. Let us recall here Veldman’s concluding remarks [23], as the main motivation behind our code CASCAD Expert development. Prandltl’s 1904 boundary-layer theory formed the starting point for the viscous-inviscid interaction methods that have been developed in the last two dacades of the 20th century. They have become very popular, since in comparison with brute-force Navier-Stokes solutions they are about two orders less expensive, whereas for flow conditions with thin shear layers the results are equally useful. Because of their modest computational complexity, they are ideal for use in aero- hydro-dynamic airfoil and airfoil cascades optimization studies. The great challenge has been to understand and resolve the singularity at separation, which occurs when the boundary-layer equations are solved with prescribed pressure. The current development of the CASCADExpert software is to add the viscous-inviscid interaction similar to the algorithm described in [7], [8] and [24]. 3. FINITE ELEMENT METHOD FOR CASCADE HYDRODYNAMICS The boundary value problem for the streamfunction, without or with viscous-inviscid coupling, has to be solved numerically. There are two main approaches: • Solve a boundary integral equation in direct formulation (Boundary Element Method) or indirect formulation (Panel Method) • Solve the boundary value problem in the computational domain from Fig. 1 using a classical Finite Element Method (FEM). In the CASCADExpert code we are using the domain discretization using the FEM, since it provides additional advantages in terms of post-processing domain data. Moreover, the state-of-the-art sparse matrix storage and iterative solvers for large systems of equations [1] have eliminated the gap between the boundary methods (which claim a small number of equations, but have a dense matrix) and the domain methods (which have a much larger number of equations but a very sparse matrix). The starting point in the Finite Element Method is the equivalent variational formulation of the boundary value problem for the streamfunction, 2 1  ∂ψ   ∂ψ F (ψ ) = ∫   + 2 Ω  ∂x   ∂y 

  

2

  dx dy  ,

(22)

− ∫ψ f N dΓ ΓN

where Ω denotes the computational domain, and ΓN is the part of the domain boundary where Neumann condition is imposed, with Neumann data f N . The solution of the original problem will be the one that

realizes the extremum of functional (22), while satisfying the Dirichlet conditions. Note that by using this approach we have only first order partial derivatives instead of second order in the original Laplace equation (2). The computational domain is then discretized with triangular elements, resulting in an unstructured mesh. Then, the streamfunction is approximated using its nodal values ψ j and the finite element shape functions N j (x, y ) , as follows

ψ (x, y ) ≅ ∑ψ j N j (x, y ) ,

(23)

j

Introducing the finite element approximation (23) in (21) transforms the functional in a function depending on the streamfunction nodal values, F (ψ 1 ,ψ 2 ,!) = 2 2 ∂N j   ∂N j   1       ∑ψ j ∂x  +  ∑ψ j ∂y   dx dy 2 Ω∫  j   j   

(24)

  − ∫  ∑ψ j N j  f N dΓ ΓN  j  The extremum of the above function is found by canceling all partial derivatives with respect to the nodal values, leading to the following system of equations:  ∂N ∂N j ∂N i ∂N j  ∂F  dx dy = ∑ψ j ∫  i + ∂ψ i ∂ ∂ ∂ ∂ x x y y j  Ω (25) − ∫ N i f N dΓ = 0 ΓN

The main advantage in using the FEM is found in the implementation of Neumann conditions. For example, at the trailing edge where both impenetrability and the Joukowski conditions must be imposed, we have the following two equations in the system:

ψ TE = 0 and (26)  ∂N TE ∂N j ∂N TE ∂N j   dx dy = 0 + ∂x ∂x ∂y ∂y  j Ω We immediately note that there is no need to actually compute the normal derivative at the sharp trailing edge. In fact, instead of canceling the normal derivative point-wise, we are canceling the local flux. This also can be shown to reduce to the original Kutta-Joukowski condition in the case of cusp trailing edge (as is the case of Joukowski airfoil). The periodic Neumann condition (5), or the pressure continuity across the wake (13), are implemented simply by adding up two rows from the assembled matrix [20]. Once again, there is no need to evaluate the normal derivative, which is an important advantage in terms of accuracy (especially on unstructured meshes) and robustness:

∑ψ j ∫ 

Proceedings of the Workshop on VORTEX DOMINATED FLOWS. ACHIEVEMENTS AND OPEN PROBLEMS, Timisoara, Romania, June 10-11, 2005

 ∂N Pup ∂N j

∂N Pup ∂N j   dx dy +  x x y y ∂ ∂ ∂ ∂ j Ω  (28)  ∂N Plo ∂N j ∂N Plo ∂N j  ∑ψ j ∫  ∂x ∂x + ∂y ∂y  dx dy = 0 j  Ω The actual implementation [11] of our finite element solver is done using the Portable, Extensible Toolkit for Scientific Computation (PETSc), [1].

∑ψ j ∫ 

+

4. HYDROFOIL SHAPE REPRESENTATION An efficient and robust solver like the one described in the previous section must be complemented with an appropriate procedure to automatically build the computational domain, or to modify it according to a design/optimization algorithm. A key ingredient in this pre-processing phase is the parametric description of the hydrofoil. In this section we present the mathematical formulation as well as the implementation of the least-squares modal representation of the hydrofoil, as used in our code CASCADExpert. Miller et al. [15] use a non-uniform rational B-spline (NURBS) to represent a blade design. Their interactive design methodology uses blade sections defined with respect to general surfaces of revolution generally obtained as axisymmetric streamsurfaces in a Q3D turbomachine flow analysis. These algorithms are further developed by Miller in [16], where the geometrical approach for describing blade sections and the process of mapping curves from 2D to 3D space are implemented using B-spline curves and in some cases NURBS curves. The final result is a Bspline surface for the blade. Dahl [6] also uses a single B-spline curve to represent the airfoil shape in a design procedure for wind turbine airfoils. A 35 B-spline control points was considered large enough by Li [14] to accommodate free-form geometry changes in an optimization algorithm. Iollo and Salas [13] use the inverse Theodorsen transform [22] in order to parameterize the airfoil in a procedure of determining airfoils that approximate, in a least squares sense, a given surface pressure distribution. However, in this approach it is not clear how many Fourier coefficients should be taken into account since one does not know how fast the series will converge. 4.1. Modal least-squares representation of hydrofoils The goal of achieving hydrofoil/airfoil shape flexibility with a small number of coefficients, while maintaining the required smoothness was addressed by Chang et al. [4]. The starting point in their approach is a set of functions inspired from the NACA family. Essentially, this is a polynomial basis, with the x added:

65

g1 (x ) = x − x g 2 (x ) = x (1 − x )

g 3 (x ) = x 2 (1 − x )

(29)

g 4 (x ) = x 3 (1 − x )

g 5 (x ) = x 4 (1 − x ) g 6 (x ) = x 5 (1 − x )

The above linearly independent set of functions is orthogonalized using the Gramm-Schmidt process, as follows: h1 (x ) = g1 (x ), h2 (x ) = g 2 (x ) − α 21h1 (x ),

""""""""""

(30)

m −1

hm (x ) = g m (x ) − ∑ α mi hi (x ), i =1

where the α mi coefficients are 1

α mi =

∫ g m (x )hi (x )dx 0

1

.

(31)

∫ (x )dx hi2

0

Each function is further normalized as hi (x ) mi (x ) = . 1

(32)

∫ (x )dx hi2

0

The resulting set of orthonormal functions are: m1 (x ) = 30 g1 (x ) 14  13  10  − g1 (x ) + g 2 (x ) 3  14  27 m3 (x ) = 1190 × 17 73  35   g1 (x ) − g 2 (x ) + g 3 (x ) 81  81  132 442 × m4 (x ) = 13 377 953  123  g1 (x ) + g 2 (x ) − g 3 (x ) − 748 748 748    + g (x )  4   m2 (x ) =

910 10582 × 111 8 5723  671  g1 (x ) − g 2 (x ) + g 3 (x )  5915 35  11830   10123  g 4 (x ) + g 5 (x ) −   5915 

(33a) (33b)

(33c)

(33d)

m5 (x ) =

(33e)

Proceedings of the Workshop on VORTEX DOMINATED FLOWS. ACHIEVEMENTS AND OPEN PROBLEMS, Timisoara, Romania, June 10-11, 2005

66

m 6 (x ) = 108 962 ×

giving the explicit expressions for the modes mi (x ) .

221 551 3467   g1 (x ) + g 2 (x ) − g 3 (x ) (33f) − 5994 5994  11988   10141  13025 g 4 (x ) − g 5 (x ) + g 6 (x ) +  5994  5994 

Since all modes have a term containing

Although the procedure for constructing the set of orthonormal functions is outlined in [4], here we are

slope at origin behaves like 1 / x , thus allowing for the same tangent direction on both suction side and pressure side curves. Fig. 3 shows the six modes given by (33a)…(33f). 3

2

1.5

x , the

1st mode

2nd mode

2

3rd mode

1 1

1 0 0 0.5

−1 −1

0

0

0.2

0.4

0.6

0.8

1

−2

0

0.2

0.4

0.6

0.8

−2

1

2

3

1

2

0

1

0

−1

0

−1

−1

−2

−2

−3

4th mode 0

0.2

0.4

0.6

0.8

1

−2

0

0.2

0.4

0.6

0.8

1

2

5th mode

0

0.2

0.4

1

0.6

0.8

1

−3

6th mode 0

0.2

0.4

0.6

0.8

1

Figure 3. The first six orthonormal modes used for hydrofoil reconstruction This is important in order to avoid angular points at leading edge. However, this requirement does not need to be satisfied at the trailing edge, where an angular point is required for the Kutta-Jukovski condition. 0.1

0.05

suction side modal pressure side modal points suction side points pressure side

0

−0.05

0

0.2

0.4

0.6

0.8

1

Figure 4. Least-squares modal reconstruction of the Clark Y hydrofoil. As an example we present the reconstruction of the Clark Y hydrofoil, Fig. 3. Essentially we are computing the coefficients a i in the expansion below, y ( x ) = ∑ a i mi (x ) i

(34)

separately for suction side and for the pressure side, respectively. In doing so, we are employing the FNLSQ subroutine from the IMSL library, which computes a least-squares approximation of a set of points using user-supplied basis functions. The modal coefficients obtained for the Clark Y hydrofoil are shown in Table 1. One can see that the coefficient values are rapidly decreasing for higher order modes thus motivating the truncation to a rather small number of modes. The overall accuracy is quantified using the sum of errors squared, as provided by the FNLSQ subroutine. Table 1. Modal coefficients for Clark Y hydrofoil Coeff. Suction side Pressure side a1 6.78268E-2 -1.962219E-2 a2 9.65698E-3 3.204551E-3 a3 -1.613487E-3 -1.070799E-3 a4 7.204205E-4 7.302335E-5 a5 -4.714986E-4 -1.120281E-4 a6 Sum of squares of the errors

4.428373E-4

3.308699E-5

6.726589E-6

3.067845E-6

Proceedings of the Workshop on VORTEX DOMINATED FLOWS. ACHIEVEMENTS AND OPEN PROBLEMS, Timisoara, Romania, June 10-11, 2005

NACA 6410, stagger 60, spacing 1, inlet 75

1

0

axial velocity tangential velocity

−1

−2 −0.5

0 0.5 1 dimensionless axial coordinate

1.5

The result of averaging the data from Fig. 5 is presented in Fig. 6. One can see that the corresponding curves are really smooth and correctly represent the desired V x (x ) and V y (x ) functions. The average axial velocity clearly reveals the blade blockage effect, while the average tangential velocity shows the cascade loading from leading edge to trailing edge. These curves are the starting point for computation of the body force which replaces the presence of the actual cascade when computing the axisymmetric turbomachinery swirling flow.

1 dimensionless velocity

250 points on the foil, axial force 0.553, tang. force 1.045 2

Figure 6. Axial and tangential velocity averaged in the direction parallel to the cascade front-line..

250 points on the foil, axial force 0.553, tang. force 1.045 2

0

−1

−2 −0.5

NACA 6410, stagger 60, spacing 1, inlet 75

dimensionless velocity

5. THE Q2D CASCADE MODEL In quasi-three-dimensional (Q3D) turbo-machinery flow model, the flow equations are averaged in circumferential direction, such that only the axial and radial coordinated remain as independent spatial variables. As a result, even if one is computing the flow through cascades on axisymmetric flow surfaces, the axisymmetric flow solver that provides the shape of these streamsurfaces does not account for any circumferential variations. Thus the need to obtain averaged data from cascade computations in order to provide the swirling flow solver with the required law of variation of circulation along the meridional streamlines.

67

0 0.5 1 dimensionless axial coordinate

1.5

Figure 5. Axial and tangential velocity raw data, as computed on the unstructured triangulation fo the computational domain. For example, Fig. 5 shows the axial and tangential velocity components as computed in the centroids of the finite elements, and plotted with respect to the axial coordinate in the abscissa. One can see that the flow distortion in the cascade is seen here a scattered distribution of data. Our goal is to average these data in order to obtain two smooth curves representing the V x (x ) and V y (x ) functions. Note that the velocity components are obtained by computing the partial derivatives of the streamfunction, Eq. 1. When using piecewise linear triangular finite elements, the velocity components will be piecewise constant. When nodal values for velocity components are needed, we have employed the Superconvergent Patch Recovery [10] in order to obtain the best accuracy. The V x (x ) and V y (x ) functions are represented as piecewise polynomials using the BSLSQ subroutine from the IMSL library, which computes a least-squares spline approximation, and returns the B-spline coefficients. The weight vector contains for each point the area of the triangular element where the corresponding velocity components are computed according to Eq. (1).

6. THE CASCADEXPERT CODE In this section we briefly review the present status of our CASCADExpert code, which implements the theoretical developments presented above. The code development [21], [9] and validation [2] led to its first public release. The code has a graphic user interface that allows a quick and easy introduction of input data (hydrofoil shape – including the UIUC extensive database [18], stagger, pitch, inlet angle), and also a set of qualitative/quantitative postprocessing capabilities. Among them, the user can define and plot the pressure coefficient distribution, and can display and examine the streamline pattern. The first major development will be the viscous-inviscid computation, in order to add the capability of evaluating hydraulic losses and optimization from energy viepoint. 7. CONCLUSIONS The paper presents the mathematical formulation in streamfunction for the cascade hydrodynamics. Besides the inviscid, incompressible and irrotational flow we derive the boundary conditions for streamfunction on the hydrofoil and wake in order to implement a viscous-inviscid interaction algorithm. The FEM used for solving the boundary value problem for streamfunction, is revisited with emphasis on the implementation of boundary conditions involving normal derivatives. For the pre-processing part we

68

Proceedings of the Workshop on VORTEX DOMINATED FLOWS. ACHIEVEMENTS AND OPEN PROBLEMS, Timisoara, Romania, June 10-11, 2005

have implemented a robust and flexible approach for least-squares modal reconstruction of the hydrofoil shape. Thanks to its accuracy and ability to produce correct streamlined smooth shapes, this method is well suited for hydrofoil design. Moreover, since a small number of parameters are involved in the hydrofoil parametrization, this representation is very useful within an optimization procedure. For the post-processing part we are presenting an averaging procedure for the velocity components in order to obtain the information required in a Q3D turbomachinery computation. Our method preserves high accuracy by using a weighted piecewise B-splines least-squares approximation of velocity data from the unstructured mesh. Finally, we review the current status in the development of our CASCADExpert software, and underline the futures developments. ACKNOWLEDGMENTS The present work has been supported from the National University Research Council (CNCSIS) Consortium Grant 33/2005 “Vortex Hydrodynamics and Applications”. The CASCADExpert software is developed at the National Center for Engineering of Systems with Complex Fluids, Politehnica University of Timişoara, in association with Vital Soft Solutions, Timişoara. The first author would like to thank EPFLLMH where he completed part of this work as Visiting Professor in April-May 2005. REFERENCES 1. Balay S., Buschelman K., Eijkhout V., Gropp W., Kaushik D., Knepley M., Curfman McInnes L., Smith B., Zhang H. (2004) PETSc Users Manual, Argonne National Laboratory, ANL-95/11 Revision 2.2.1. 2. Balint D. and Susan-Resiga R. (2000) Numerical Investigation of a NACA 65 Series Low Speed Compressor, Buletinul Ştiinţific al Univ. “Politehnica” Timişoara, Seria Mecanică, Tom 45(59), pp. 7-14. 3. Bradshaw P., Cebeci T., Whitelaw J. H. (1981) Engineering Calculation Methods for Turbulent Flow, Academic Press. 4. Chang I.C., Torres F.J., Tung C. (1995) Geometric Analysis of Wing Sections, NASA Technical Memorandum 110346. 5. Coenen E. G. M. (2001) Viscous-Inviscid Interaction with the Quasi-Simultaneous Method for 2D and 3D Aerodynamic Flow, PhD thesis, Rijksuniversiteit Groningen. 6. Dahl K. S., Fuglsang P. (1998) Design of the Wind Turbine Airfoil Family RISO-A-XX, Riso National Laboratory, Roskilde, Denmark, Riso-R-1024(EN). 7. Drela M. (1989) XFOIL: An Analysis and Design System for Low Reynolds Number Airfoils, in Mueller T. J. (ed.), Low Reynolds Number Aerodynamics, Lecture Notes in Engineering #54, Springer Verlag, pp.1-12. 8. Drela M. and Giles M. B. (1987) Viscous-Inviscid Analysis of Transonic and Low Reynolds Number Airfoils, AIAA Journal 25(10), pp. 1347-1355.

9. Frunză T., Susan-Resiga R. (2002) Finite Element Analysis of Cascade Flows, in Proceedings of the “Classics and Fashion in Fluid Machinery” Conference, Beograd, pp. 179-186. 10. Frunză T., Susan-Resiga R. (2003) Superconvergent Patch Recovery for Velocity Computation on Unstructured Mesh and Cascade Flow Application, Proceedings of the Workshop on Numerical Methods in Fluid Mechanics and FLUENT Applications, Ed. Orizonturi Universitare, Timişoara, pp. 64-75. 11. Frunză T., Susan-Resiga R. (2001) Software Development for Cascade Flow Simulation, in Proceedings of the Workshop on Numerical Simulation for Fluid Mechanics and Magnetic Liquids, Ed. Orizonturi Universitare, Timişoara, pp. 57-68. 12. Gostelow J. P. (1984) Cascade Aerodynamics, Pergamon Press. 13. Iollo A., Salas M. D. (1999) Optimum transonic airfoils based on the Euler equations, Computers & Fluids, 28, 653-674 14. Li, W. (2003) Profile Optimization Method for Robust Airfoil Shape Optimization in Viscous Flow, NASA TM-2003-212408. 15. Miller P. L., Oliver J. H., Miller D. P., Tweedt D. L. (1996) BladeCAD: An Interactive Geometric Design Tool for Turbomachinery Blades, NASA Technical Memorandum 107262. 16. Miller P. L. (2000), Blade geometry description using B-splines and general surfaces of revolution, PhD thesis, Iowa State University. 17. Milne-Thomson, L. M. (1968) Theoretical Hydrodynamics, The MacMillan Press, 5th edition. 18. Selig M. (2002) UIUC Airfoil2 Data Site, http://www.aae.uiuc.edu/m-selig/ads.html, University of Illinois at Urbana-Champain, Department of Aerospace Engineering. 1. Susan-Resiga R. (2003) Numerical Fluid Mechanics, Orizonturi Universitare Publishing House, Timişoara, ISBN 973-638-014-9 (in Romanian). 20. Susan-Resiga R. and Muntean S. (1999) Periodic Boundary Conditions Implementation for the Finite Element Analysis of the Cascade, Buletinul Ştiinţific al Univ. “Politehnica” Timişoara, Seria Mecanică, Tom 44(58), pp. 151-160. 21. Susan-Resiga R., Muntean S., Anton I. (2000) Numerical Analysis of Cascade Flow. Part I: Finite Element Analysis of the Inviscid Flow, Buletinul Ştiinţific al Univ. “Politehnica” Timişoara, Seria Mecanică, Tom 45(59), pp. 159-166. 22. Theodorsen T. (1933) Theory of Wing Sections of Arbitrary Shape, NACA Report 411. 23. Veldman A. E. P. (2001) Matched asymptotic expansions and the numerical treatment of viscous-inviscid interaction, Journal of Engineering Mathematics, 39, pp. 189-206. 24. Yiu K.F.C. and Giles M. B. (1995) Simultaneous Viscous-Inviscid Coupling via Transpiration, Journal of Computational Physics 120, pp. 157-170. 2 In the UIUC database coordinates are in x,y format, starting from trailing edge, along the upper surface to the leading edge and back around the lower surface to trailing edge, i.e. counterclockwise. The CASCADExpert convention is that points on hydrofoil are ordered clockwise (positive τ direction), starting from and ending at trailing edge.