VIEWS: 427 PAGES: 10 POSTED ON: 7/10/2011 Public Domain
Scientific Bulletin of the 3rd Workshop on “Politehnica” University of Timisoara Vortex Dominated Flows Transactions on Mechanics Timisoara, Romania Tom 52(66), Fascicola 3, 2007 June 1 - 2, 2007 COMSOL MULTIPHYSICS VERSUS FLUENT: 2D NUMERICAL SIMULATION OF THE STATIONARY FLOW AROUND A BLADE OF THE ACHARD TURBINE Andrei-Mugur GEORGESCU, Assoc. Prof.* Sanda-Carmen GEORGESCU, Assoc. Prof. Hydraulics and Environmental Protection Department Hydraulics and Hydraulic Machinery Department Technical Civil Engineering University Bucharest University “Politehnica” of Bucharest Sandor BERNAD, Senior Researcher Costin Ioan COŞOIU, Assist. Prof. Center of Advanced Research in Engineering Sciences Hydraulics and Environmental Protection Department Romanian Academy - Timisoara Branch Technical Civil Engineering University Bucharest *Corresponding author: Bd Lacul Tei 24, Sector 2, 020396, Bucharest, Romania Tel.: (+40) 212433660, Fax: (+40) 212433660, Email: andreig@mail.utcb.ro ABSTRACT x, y, z [m] Cartesian coordinates α [grd] angle of attack Two-dimensional numerical modelling of the sta- tionary flow around a blade of the Achard turbine, a λ [−] tip speed ratio new water-current turbine concept, is performed both ν [m2/s] water cinematic viscosity with COMSOL Multiphysics 3.3 and with Fluent 6.0.1 ω [rad/s] rotational velocity software, in order to compare the results and the sof- ρ [kg/m3] water density tware capabilities. The k − ε turbulence model has θ [grd] azimuth angle of the blade been selected and same geometry and boundary con- ditions were considered within computations. Subscripts and Superscripts D drag (along flow direction) KEYWORDS L lift (cross flow direction) Achard turbine, cross-flow marine turbine, NACA * dimensionless variable airfoil 1. INTRODUCTION NOMENCLATURE In 2001, the Geophysical and Industrial Fluid Flows c [m] airfoil chord length Laboratory (LEGI) of Grenoble, France, launched the cD [-] drag coefficient HARVEST Project (Hydroliennes à Axe de Rotation cL [-] lift coefficient VErtical STabilisé), to develop a suitable technology for marine and river hydro-power farms using cross- c0 [m] airfoil mean camber line length flow current energy converters piled up in towers [1, d [m] airfoil maximum thickness 2 & 8]. The hydro-dynamic of these systems is studied F [N] hydrodynamic force at LEGI with the support of the R&D Division of the FD [N] drag force EDF Group, while other laboratories of the Rhône- FL [N] lift force Alpes Region are charged with mechanical aspects g [m/s2] gravity (3S-INP of Grenoble and LDMS-INSA of Lyon), as H [m] turbine height well as with electrical aspects (LEG-INP Grenoble). k [m2/s2] turbulent kinetic energy In 2006, the Technical Civil Engineering Univer- p [Pa] pressure sity Bucharest, in collaboration with the University R [m] turbine radius “Politehnica” of Bucharest and the Romanian Academy - Timisoara Branch, started the THARVEST Project, Rec [−] chord based Reynolds number within the CEEX Program sustained by the Roma- u = ωR [m/s] transport velocity nian Ministry of Education and Research [5]. The V0 [m/s] upstream velocity THARVEST Project aims to study experimentally and w [m/s] relative velocity numerically the hydro-dynamics of a new concept 14 Proceedings of the 3 rd Workshop on Vortex Dominated Flows. Achievements and Open Problems, Timisoara, Romania, June 1 - 2, 2007 of water-current turbine, called Achard turbine, in radial supports are shaped with NACA 4518 airfoils, collaboration with the LEGI partners involved in the while the circular rims are shaped with lens type airfoil. HARVEST Project. The Achard turbine, a cross-flow The turbine radius is R = 0.5 m, and the turbine marine or river turbine with vertical axis and delta height is H = 1 m. In Figure 2 we present the Achard blades is studied in France mainly with regard to marine turbine geometry generated in MATLAB (the upper applications, to extract energy from tidal currents in and lower rims are not represented here). Along each costal locations. But the Achard turbines are also delta blade, the airfoil mean camber line length c0 suitable to be placed in big rivers, as the Danube, and varies from 0.18 m at z = 0 , to 0.12 m at the ex- to produce the desired power by summing elementary tremities, where z = ±0.5 m. Between the leading power provided by small turbine modules [6 & 8]. edge of the blade’s extremity and the leading edge of the blade at mid-height of the turbine, there is a 30º azimuth angle. Figure 1. The Achard turbine [LEGI courtesy] Figure 2. The Achard turbine geometry In Figure 1 we present the Achard turbine that is studied now in Romania. The turbine description and At z = ± H 4 = ±0.25 m, the horizontal cross- the generation of its geometry are explained within section of the runner gives airfoils with c0 = 0.15 m the section 2 of this paper. (the mean value of the mean camber line length In this paper we focus on the 2D numerical mod- along the delta wing), and with the chord length elling of the stationary flow around a blade of the c = 0.1494 m (see section 3). Within this paper, the Achard turbine. The blades are shaped with NACA 2D computations correspond to the cross-plane airfoils. The 2D modelling corresponds to the flow placed at z = 0.25 m level (see Figure 3), where the around an airfoil in a horizontal cross-section of the runner, situated at one quarter of the turbine height, 0.5 starting from the top. Different azimuth angles of the o 0.4 θ=0 blade are analysed. Two types of airfoils are considered 0.3 within the section 3: a curved one (NACA 4518) and a straight one (NACA 0018). The simulations are 0.2 performed both with COMSOL Multiphysics 3.3 0.1 software and with Fluent 6.01 software. The results y [m] 0 and the software capabilities are compared within the section 4 with experimental data. The paper con- −0.1 clusions are summarized in section 5. −0.2 o θ = 120 o θ = 240 2. ACHARD TURBINE DESCRIPTION −0.3 The vertical axis Achard turbine from Figure 1 −0.4 consists of a runner with three vertical delta blades, −0.5 sustained by radial supports at the mid-height of the −0.5 −0.4 −0.3 −0.2 −0.1 0 x [m] 0.1 0.2 0.3 0.4 0.5 turbine, and stiffened with circular rims at the upper and lower part of the turbine. The blades and their Figure 3. Computational runner cross-section Proceedings of the 3 rd Workshop on Vortex Dominated Flows. Achievements and Open Problems, Timisoara, Romania, June 1 - 2, 2007 15 three blade profiles have a mean camber line length Its dimensionless value is m ∗ = m c = 0.03758 . As of c0 = 0.15 m. The values of the azimuth angle of percentage of the chord, m = 3.758% ≅ 4% , so the { the blades in Figure 3 are θ = 0 o ; 120 o ; 240 o , in } first digit of the NACA airfoil is 4. counter clockwise direction. (a) NACA 4518 0.02 3. AIRFOIL SELECTION 0.01 y [m] 0 The runner blades are shaped with NACA airfoils −0.01 of four-digit series [13], where the first digit is the −0.02 maximum upper camber m (as percentage of the chord), 0 0.05 0.1 0.15 x [m] the second digit is the distance p of the maximum upper camber from the airfoil leading edge (in tens of 0.02 (b) NACA 0018 percents of the chord), and the last two digits describe 0.01 the maximum thickness of the airfoil, d, as percent y [m] 0 of the chord length. For the Achard turbine blades, −0.01 we consider d = 18 %. In a xOy plane, the airfoil coordinates {x, y} are −0.02 0 0.05 0.1 0.15 x [m] nondimensionalised with respect to the chord length Figure 4. Blade profiles for c0 = 0.15 m: (a) curved c, with x ∗ = 0 at the leading edge and with x ∗ = 1 airfoil NACA 4518; (b) straight airfoil NACA 0018 at the trailing edge (the dimensionless variables are denoted with an asterisk). Thus, the airfoil type corresponding to the Achard ∗ { The dimensionless coordinates x ∗ , y s of the airfoil } turbine blades is NACA 4518, an airfoil with the mean camber line along the runner circumference mean camber line are defined as: (Figure 4a). The computations from section 4, are performed for the profile NACA 4518 and also for the ∗ ys = m∗ ∗2 (2 p ∗ ) 2 − 1 x ∗ for 0 ≤ x ∗ < m ∗ , and straight profile NACA 0018, which can be generated p ∗ from (2) for y s = 0 , since m * = 0 in (1). For the ∗ ys = m∗ ( ) ⎡ 1 − 2 p ∗ + 2 p ∗ x ∗ - x ∗ 2 ⎤ for (1) NACA 0018 airfoil, the chord length is c = c0 = 0.15 m (1 − p ) ∗ 2 ⎢ ⎣ ⎥ ⎦ (Figure 4b). That last choice is due to the fact that experimental and numerical data are available for such m∗ ≤ x∗ ≤ 1 , a straight profile, corresponding to a Darrieus marine where p ∗ = 0.5 , since the mean camber line of the turbine, a vertical axis cross-flow turbine with two straight blades of 0.15 m chord length [3, 4, 7 & 9]. airfoil is along the circle of radius R = 0.5 m, as in Figure 3. So, the second digit of the NACA airfoil is 4. NUMERICAL RESULTS 5, because the maximum upper camber is placed at a half-distance between the leading edge and the The numerical simulations are performed with trailing edge. COMSOL Multiphysics 3.3 [11] (a Finite Element { } The coordinates x ∗ , y ∗ of the upper and lower Method based software), and with Fluent 6.01 [12] (a Finite Volume Method based software), using the surfaces of the NACA airfoil are defined by: k − ε turbulence model. The stationary 2D flow d ⎛ around a blade of the Achard turbine is analysed both ∗ y∗ = ys ± ⎜ 0.29690 x ∗ − 0.12600 x ∗ − 20 ⎝ for the profile NACA 4518 (Figure 4a), for which (2) the mean camber line is along the circle of radius − 0.35160x ∗ + 0.28430x ∗ − 0.10150x ∗ ⎞ . 2 3 4 ⎟ R = 0.5 m, and for the straight profile NACA 0018 ⎠ (Figure 4b), for which other results are available in The chord length c can be expressed upon the literature. runner radius and the mean camber line length c0 : 4.1. Geometry of the models c = 2 R sin (c0 2 R ) . (3) The 3D blades geometry is realised in MATLAB as in Figure 2. At a certain z-level, each blade profile For R = 0.5 m and c0 = 0.15 m, we obtain the mean is generated in 101 nodes on the upper airfoil surface, chord length c = 0.1494 m at z = 0.25 m. and 101 nodes on the lower airfoil surface, using (2). The maximum upper camber is defined as: The dimensional (x, y) coordinates of those 202 nodes m = R (1 − cos(c0 2 R )) . (4) are then imported within the two software used for 16 Proceedings of the 3 rd Workshop on Vortex Dominated Flows. Achievements and Open Problems, Timisoara, Romania, June 1 - 2, 2007 the flow simulation. The 2D flow domain extent in the on a PC with 2 GB of memory and 3 GHz Intel proc- xOy plane is scaled with respect to the NACA profiles essor, while in Fluent it is less than 25 minutes, on a mean camber line length that is c0 = 0.15 m. The PC with 4 GB RAM and 3.2 GHz dual-core Intel origin of the system is placed on the profile chord, at Xeon processor. one quarter from the leading edge. The rectangular Within both software, the following boundary con- flow domain extents ranges from − 15c0 to + 25c0 on ditions are considered: At the left side of the domain, on the water inflow boundary, a constant upstream x-direction, and from − 3.5c0 to + 3.5c0 on y-direction. velocity V0 = 4.71 m/s and a turbulent intensity of Both COMSOL Multiphysics software and Fluent 2% are imposed. At the right side of the domain, on software use the same geometry of the flow domain. the water outflow boundary, a zero relative pressure The discretization has the same number of nodes on is considered. On the upper wall, as well as on the the profiles, but the total number of cells is slightly lower wall of the domain, a slip symmetry condition different from one software to the other (see Figure 5). is selected. On the profile surface, a logarithmic wall The meshing consists of triangular cells in COMSOL function is selected. In Fluent, a discretization of sec- and quadrilateral cells in Fluent. For about 33000 cells, ond order is used both for the velocity and pressure. the resulting run-time in COMSOL exceeds 10 hours, Figure 5. Zoom of the domain discretization at α = 0 o , in COMSOL Multiphysics (upper image) and in Fluent (lower image) Proceedings of the 3 rd Workshop on Vortex Dominated Flows. Achievements and Open Problems, Timisoara, Romania, June 1 - 2, 2007 17 The chord based Reynolds number is defined as within the range V0 ≤ w ≤ 3V0 . The velocity ratio Rec = w c ν , where w is the relative velocity. The tip w V0 upon the azimuth angle, that is w V0 = f (θ ) , speed ratio λ = ωR V0 has the imposed value λ = 2 . as well as the angle of attack variation α = α (θ ) are The relative velocity w on the blade at the leading presented in Figure 6. edge is obtained by composing the upstream velocity In this paper, most of the results are computed V0 and the transport velocity u = ωR = λ V0 : for discrete values of the pair of angles {α ;θ } . Due to the symmetry, for the straight profile NACA 0018, w = V0 1 + 2λ cosθ + λ2 = V0 1 + 4(1 + cosθ ) , (5) the flow behaviour obtained at a certain positive α where the azimuth angle θ defines the position of the value, is identical to the one obtained at the corre- blade around the circle, in counter clockwise direction. sponding negative α value. For our simulations, the Reynolds numbers exceed 4.2. COMSOL Multiphysics versus Fluent 7 ⋅ 10 5 for the whole range of the azimuth angle, We present the results, namely the pressure field [ ] θ ∈ 0 o ; 360 o , placing the phenomenon behaviour and/or the velocity field around the blades, obtained within the self modelling region with respect to the with both software, for the two types of profiles, Reynolds number. In this paper, 2D numerical com- NACA 4518 and NACA 0018, at different azimuth putations are performed on a fixed blade, the value of angles θ (see Figures 7÷9). the upstream velocity being taken so that the Reynolds The results show that the flow behaviour descrip- number on the fixed blade exceeds 10 5 , thus the flow tion is similar in both COMSOL Multi-physics and may be assumed to have the same characteristics as Fluent software. Unfortunately, it is difficult to fit the in the real rotating case. same colour scheme on both software post-processors. The angle of attack α , considered between the Fluent seems to give more accurate results. For chord and the w-direction at the leading edge, is also MATLAB users, the advantage of using COMSOL related to the azimuth angle θ through the following Multiphysics is obvious in pre-processing step, where relations: the geometry can be generated directly with MATLAB sources. α = arccos (λ 2 ) − 1 V02 + w 2 − arcsin c , for 4.3. Influence of the number of elements 2λ V0 w 2R Within this paragraph we will discuss the influence [ θ ∈ 0 o ; 180 o , and ) (6) of the number of computational cells. The simulations α = arccos (λ 2 ) − 1 V02+w 2 + arcsin c , for realised in COMSOL Multiphysics need less than 10 2λ V0 w 2R minutes when the grid consists of about 17000 cells, [ θ ∈ 180 o ; 360 o . ) while the run-time is 60 times greater for twice as many cells (e.g. about 10 h for near 33000 cells). In Figure 10 we present the results obtained using 4 40 COMSOL for the velocity field around the blade velocity ratio profiles NACA 4518 and NACA 0018, at α = 30 o ( ) 3.5 30 angle of attack 3 20 θ = 120 o , for about 16500 cells, and 32800 cells respectively. It seems that the flow behaviour (e.g. the 2.5 10 vortices development) doesn’t change substantially when multiplying by two the number of elements. w/V0 2 0 α 1.5 -10 So, to get a general view of the flow, computations with about 17000 cells can be acceptable due to the 1 -20 significant run-time economy. 0.5 -30 Of course, as long as pressure coefficients are concerned, the values obtained for 33000 cells are 0 -40 closer to the experimental existing results. 0 90 180 270 360 θ Highly accurate computations are reported by Ervin et al [4] for the flow within a Darrieus marine Figure 6. Velocity ratio w V0 [−] and angle of turbine with two-blades of NACA 0018 profile. They attack α [grd] versus the azimuth angle θ [grd] used the Turb’Flow software, which is less dissipative than COMSOL and Fluent. For more than 140000 cells, So, the resulting angle of attack varies within the the computations are really time consuming (about ( ) range α ∈ − 30o ; 30o , more precisely, from − 29.85 o 50 days) on a workstation with 8 GB of memory and to + 29.85 o . For λ = 2 , the relative velocity (5) varies 3.2 GHz dual-core Intel Xeon processor. 18 Proceedings of the 3 rd Workshop on Vortex Dominated Flows. Achievements and Open Problems, Timisoara, Romania, June 1 - 2, 2007 Figure 7. NACA 4518 profile (left images) and NACA 0018 profile (right images), at θ = 0 o , in COMSOL Multiphysics (1st & 3rd row) and in Fluent (2nd & 4th row): pressure (first 2 rows) and velocity (last 2 rows) Proceedings of the 3 rd Workshop on Vortex Dominated Flows. Achievements and Open Problems, Timisoara, Romania, June 1 - 2, 2007 19 Figure 8. Velocity field for NACA 4518 profile (left images) and NACA 0018 profile (right images), at θ = 60 o , in COMSOL Multiphysics (upper row) and in Fluent (lower row) Figure 9. Pressure field (left images) and velocity field (right images) for NACA 4518, at θ = 120 o , in COMSOL Multiphysics (upper row) and in Fluent (lower row) 20 Proceedings of the 3 rd Workshop on Vortex Dominated Flows. Achievements and Open Problems, Timisoara, Romania, June 1 - 2, 2007 Figure 10. Velocity field around the NACA 4518 (left images) and NACA 0018 (right images), at α = 30 o ( ) θ = 120 o , in COMSOL Multiphysics, for about 16500 cells (upper row) and 32800 cells (lower row) Figure 11. Pressure field (left image) and velocity field (right image) around the NACA 4518, at α = −30 o ( ) θ = 240 o , in COMSOL Multiphysics (about 33000 cells) In Figure 11 we present the results obtained using the flow around the NACA 4518 and NACA 0018 COMSOL for the pressure field and the velocity field, profiles, at different values of the angle of attack α . around the blade profile NACA 4518, at the angle of The drag coefficient c D and the lift coefficient attack α = −30 o and azimuth angle θ = 240 o , for c L are computed for each profile, based on the ele- about 33000 cells. The flow behaviour can be com- mentary drag force (dFD ) distribution, and on the pared with the one obtained for the blade placed in mirror position (with respect to the flow direction) in elementary lift force (dFL ) distribution around the the Figure 10 (see left image, lower row). profile boundary, using the following relations: 4.4. Results at different angle of attacks 2 2 Within this paragraph we present some results cD = ρ c w2 ∫ dF D and c L = ρ c w2 ∫ dF L . (7) obtained with COMSOL Multiphysics software for Proceedings of the 3 rd Workshop on Vortex Dominated Flows. Achievements and Open Problems, Timisoara, Romania, June 1 - 2, 2007 21 We mention that the above curvilinear integrals, In Figure 13, our computed results for the NACA which can be computed within the software pre- 4518 profile are plotted together with our results for processor, give the resulting drag force and lift force NACA 0018, in order to show that there are no on the profile, as: significant differences between the drag and lift coefficient values, due to the slight curvature of the ∫ FD = dFD and FL = dFL . ∫ (8) 4518 airfoil with respect to the straight-one. The numerical values of the drag coefficients and 5. CONCLUSIONS lift coefficients are compared in Figure 12 with the experimental data obtained at Sandia Laboratories In this paper, 2D numerical computations are [10] for the NACA 0018 airfoil, both in c D = c D (α ) performed both with COMSOL Multiphysics 3.3 software and Fluent 6.0.1 software, in order to depict and c L = c L (α ) representations, but also in polar the stationary flow behaviour around a cross-section diagram c D = c D (c L , α ) . We notice that there are of a fixed blade of the Achard turbine. The value of some differences between our computed values and the upstream velocity is taken so that the Reynolds the experimental ones, especially at low angles of number on the fixed profile exceeds 10 5 , thus the attack (see the polar diagram). flow may be assumed to have the same characteris- NACA 0018 tics as in the real rotating case. 2 1.5 Sandia Labs [10] As mentioned above, both software produce simi- 1.5 COMSOL results lar results. Fluent results seem to be more accurate, c D 1 1 but this is probably due to our greater experience in using Fluent. The run-time also seems to favour the 0.5 use of Fluent for more complicated problems (i.e. 0 0.5 greater amount of grid cells), but as the two computers 0 50 100 150 200 α [grd] c L used for the numerical simulations are far from having 1.5 the same computational power, this statement must be observed with reserves. 0 1 On the other hand, for profiles that can be generated cL 0.5 −0.5 from mathematical known equations (like the NACA 0 profiles), the interoperability between COMSOL and −0.5 MATLAB software gives a plus on the geometry −1 0 50 100 150 200 −1 0 0.5 1 1.5 2 generation aspect for COMSOL. α [grd] cD Although global results with respect to the pressure Figure 12. Drag and lift coefficients versus the coefficients on the airfoils agree well with experi- angle of attack (left images), and polar diagram mental static results, it is obvious that such compu- (right image) for the NACA 0018 airfoil: numerical tational approaches cannot predict the dynamic stall results in COMSOL and experimental data [10] phenomenon that is reported in the case of vertical COMSOL Multiphysics results axis cross-flow turbines. Moreover, 2D simulations 0.8 can be acceptable for blades with constant cross- 0.6 NACA 4518 NACA 0018 sectional profiles along the z-axis, but do not permit c D an accurate description of the flow for vertical axis 0.4 cross-flow turbines with varying blade cross-section 0.2 along the z-axis like the delta blade of the Achard 0 turbine. −30 −20 −10 0 α [grd] 10 20 30 Those last two aspects are to be considered in further work performed by our team. 1.5 1 0.5 ACKNOWLEDGMENTS c L 0 Authors gratefully acknowledge the CEEX Pro- −0.5 gramme from the Romanian Ministry of Education −1 and Research, for its financial support under the −1.5 −30 −20 −10 0 α [grd] 10 20 30 THARVEST Project no. 192/2006. Special thanks are addressed to Dr Jean-Luc Achard, CNRS Research Figure 13. Variations c D = c D (α ) and c L = c L (α ) Director, and to PhD student Ervin Amet from LEGI for the NACA 4518 and 0018 airfoils Grenoble, France, for consultancy and documentation on the Achard turbine. 22 Proceedings of the 3 rd Workshop on Vortex Dominated Flows. Achievements and Open Problems, Timisoara, Romania, June 1 - 2, 2007 7. Maître T., Achard J-L., Guittet L., Ploeşteanu C. (2005) REFERENCES Marine turbine development: Numerical and experi- mental investigations. Sci. Bull. of the Politehnica 1. Achard J.-L., Imbault D., Maître T. (2005) Dispositif University of Timisoara, Transactions on Mechanics, de maintien d’une turbomachine hydraulique. Brevet vol 50 (64), pp 59-66 déposé, Code FR 05 50420, Titulaire: Institut National 8. Maître T., Antheaume S., Buvat C., Corre C., Achard Polytechnique de Grenoble, France J.-L. (2007) An innovative modeling approach to 2. Achard J.-L., Maître T. (2004) Turbomachine hydrau- optimize the design configurations of marine (river) lique. Brevet déposé, Code FR 04 50209, Titulaire: cross-flow current energy converters farm. In: Proc. Institut National Polytechnique de Grenoble, France of the 7th European Wave and Tidal Energy Conference, 3. Amet E., Pellone C., Maître T. (2006) A numerical Porto, Portugal, 7p. (to appear) approach for estimating the aerodynamic characteristics 9. Ploeşteanu C. (2004) Étude hydro-dynamique d’une of a two bladed vertical Darrieus wind turbine, Scientific type d’hydraulienne à axe vertical pour les courants Bulletin of the Politehnica University of Timisoara, marins. Thèse de Doctorat Institut National Poly- Transactions on Mechanics, vol 51 (65), pp 95-102 technique de Grenoble, France 4. Amet E., Pellone C., Maître T., Achard J.-L. (2006) 10. Sheldahl R. E., Klimas P. C. (1981) Aero-dynamic Décrochage tourbillonnaire à l’arrière des pales d’une characteristics of seven airfoil sections through 180- turbine Darrieus. In: 18ème Congrès Français de Mé- degrees angle of attack for use in aerodynamic analysis canique, Grenoble, France, 9p. (to appear) of vertical axis wind turbines, SAND80-2114, Sandia 5. Georgescu A.-M., Georgescu S.-C., Bernad S. et al. National Laboratories, USA (2006-2008) Interinfluence of the vertical axis, stabilised, 11. *** (2006) COMSOL Multiphysics 3.3. User’s Guide, Achard type hydraulic turbines (THARVEST), CEEX COMSOL AB., Stockholm, Sweden Programme, no. 192/20.07.2006, AMCSIT Politehnica 12. *** (2001) FLUENT 6. User’s Guide, Fluent Inc., Bucharest, Romania, Available at: http://hidraulica. Lebanon, USA utcb.ro/tharvest/ 13. *** (2006) UIUC Airfoil Coordinates Database. Ver- 6. Georgescu S.-C., Georgescu A.-M., Damian R. M., sion 2.0 (over 1550 airfoils), UIUC Applied Aerody- Achard J.-L. (2007) Past and future of water turbines namics Group, Department of Aerospace Engineering, in Romania. In: Proceedings of the 5th International University of Illinois, Urbana-Champaign, USA, At: Water History Association Conference on Pasts and http://www.ae.uiuc.edu/m-selig/ ads/coord_database. Futures of Water, Tampere, Finland, 8p. (to appear) html