fa Univerza v Mariboru Fakulteta za energetiko Journal of ERGY TECHNOLOGY FEBRUARY 2013 Projekt ja sofinanciran v okviru Fatrolovaga programa za izboljšanje energetske učinkovitosti. * j| Energija za življenje Projekt je sofinanciran v okviru Petrolovega programa za izboljšanje energetske učinkovitosti. JOURNAL OF ENERGY TECHNOLOGY jet) VOLUME 6 / Issue 1 Revija Journal of Energy Technology (JET) je indeksirana v naslednjih bazah: INSPEC©, Cambridge Scientific Abstracts: Abstracts in New Technologies and Engineering (CSA ANTE), ProQuest's Technology Research Database. The Journal of Energy Technology (JET) is indexed and abstracted in the following databases: INSPEC©, Cambridge Scientific Abstracts: Abstracts in New Technologies and Engineering (CSA ANTE), ProQuest's Technology Research Database. JOURNAL OF ENERGY TECHNOLOGY Ustanovitelji / FOUNDERS Fakulteta za energetiko, UNIVERZA V MARIBORU / FACULTY OF ENERGY TECHNOLOGY, UNIVERSITY OF MARIBOR Izdajatelj / PUBLISHER Fakulteta za energetiko, UNIVERZA V MARIBORU / FACULTY OF ENERGY TECHNOLOGY, UNIVERSITY OF MARIBOR Izdajateljski svet / PUBLISHING COUNCIL Zasl. prof. dr. Dali DONLAGIČ, Univerza v Mariboru, Slovenija, predsednik / University of Maribor, Slovenia, President Prof. dr. Bruno CVIKL, Univerza v Mariboru, Slovenija / University of Maribor, Slovenia Prof. ddr. Denis DONLAGIČ, Univerza v Mariboru, Slovenija / University of Maribor, Slovenia Prof. dr. Danilo FERETIČ, Sveučilište u Zagrebu, Hrvaška / University in Zagreb, Croatia Prof. dr. Roman KLASINC, Technische Universität Graz, Avstrija / Graz University Of Technology, Austria Prof. dr. Alfred LEIPERTZ, Universität Erlangen, Nemčija / University of Erlangen, Germany Prof. dr. Milan MARČIČ, Univerza v Mariboru, Slovenija / University of Maribor, Slovenia Prof. dr. Branimir MATIJAŠEVIČ, Sveučilište u Zagrebu, Hrvaška / University in Zagreb, Croatia Prof. dr. Borut MAVKO, Inštitut Jožef Stefan, Slovenija / Jozef Stefan Institute, Slovenia Prof. dr. Greg NATERER, University of Ontario, Kanada / University of Ontario, Canada Prof. dr. Enrico NOBILE, Universita degli Studi di Trieste, Italia / University of Trieste, Italy Prof. dr. Iztok POTRČ, Univerza v Mariboru, Slovenija / University of Maribor, Slovenia Prof. dr. Andrej PREDIN, Univerza v Mariboru, Slovenija / University of Maribor, Slovenia Prof. dr. Jože VORŠIČ, Univerza v Mariboru, Slovenija / University of Maribor, Slovenia Prof. dr. Koichi WATANABE, KEIO University, Japonska / KEIO University, Japan Odgovorni urednik / EDITOR-IN-CHIEF Andrej PREDIN Uredniki / CO-EDITORS Jurij AVSEC Miralem HADŽISELIMOVIČ Gorazd HREN Milan MARČIČ Jože PIHLER Iztok POTRČ Janez USENIK Peter VIRTIČ Jože VORŠIČ Uredniški odbor / EDITORIAL BOARD Prof. dr. Jurij AVSEC, Univerza v Mariboru, Slovenija / University of Maribor, Slovenia Prof. ddr. Denis OONLAGIČ, Univerza v Mariboru, Slovenija / University of Maribor, Slovenia Doc. dr. Željko HEDERIČ, Sveučilište Josipa Jurja Strossmayera u Osijeku, Hrvatska / Josip Juraj Strossmayer University Osijek, Croatia Prof. dr. Miralem HADŽISELIMOVIČ, Univerza v Mariboru, Slovenija / University of Maribor, Slovenia Doc. dr. Gorazd HREN, Univerza v Mariboru, Slovenija / University of Maribor, Slovenia Prof. dr. Roman KLASINC, Technische Universität Graz, Avstrija / Graz University Of Technology, Austria dr. Ivan Aleksander KODELI, Institut Jožef Stefan, Slovenija / Jožef Stefan Institute, Slovenia Prof. dr. Jurij KROPE, Univerza v Mariboru, Slovenija / University of Maribor, Slovenia Prof. dr. Alfred LEIPERTZ, Universität Erlangen, Nemčija / University of Erlangen, Germany Prof. dr. Branimir MATIJAŠEVIČ, Sveučilište u Zagrebu, Hrvaška / University of Zagreb, Croatia Prof. dr. Matej MENCINGER, Univerza v Mariboru, Slovenija / University of Maribor, Slovenia Prof. dr. Greg NATERER, University of Ontario, Kanada / University of Ontario, Canada Prof. dr. Enrico NOBILE, Universita degli Studi di Trieste, Italia / University of Trieste, Italy Prof. dr. Iztok POTRČ, Univerza v Mariboru, Slovenija / University of Maribor, Slovenia Prof. dr. Andrej PREDIN, Univerza v Mariboru, Slovenija / University of Maribor, Slovenia Prof. dr. Aleksandar SALJNIKOV, Univerza Beograd, Srbija / University of Beograd, Serbia Prof. dr. Brane ŠIROK, Univerza v Ljubljani, Slovenija / University of Ljubljana, Slovenia Doc. dr. Andrej TRKOV, Institut Jožef Stefan, Slovenija / Jožef Stefan Institute, Slovenia Prof. ddr. Janez USENIK, Univerza v Mariboru, Slovenija / University of Maribor, Slovenia Doc. dr. Peter VIRTIČ, Univerza v Mariboru, Slovenija / University of Maribor, Slovenia Prof. dr. Jože VORŠIČ, Univerza v Mariboru, Slovenija / University of Maribor, Slovenia Prof. dr. Koichi WATANABE, KEIO University, Japonska / KEIO University, Japan Prof. dr. Mykhailo ZAGIRNYAK, Kremenchuk Mykhailo Ostrohradskyi National University, Ukrajina / Kremenchuk Mykhailo Ostrohradskyi National University, Ukraine, Doc. dr. Tomaž ŽAGAR, Univerza v Mariboru, Slovenija / University of Maribor, Slovenia Doc. dr. Franc ŽERDIN, Univerza v Mariboru, Slovenija / University of Maribor, Slovenia Tehniška podpora / TECHNICAL SUPPORT Tamara BREČKO BOGOVČIČ, Sonja NOVAK, Janko OMERZU; Izhajanje revije / PUBLISHING Revija izhaja štirikrat letno v nakladi 150 izvodov. Članki so dostopni na spletni strani revije -www.fe.um.si/si/jet.html / The journal is published four times a year. Articles are available at the journal's home page - www.fe.um.si/en/jet.html. Cena posameznega izvoda revije (brez DDV) / price per issue (VAT not included in price): 50,00 EUR Informacije o naročninah / subscription information: http://www.fe.um.si/en/jet/subscriptions.html Lektoriranje / LANGUAGE EDITING Terry T. JACKSON Oblikovanje in tisk / DESIGN AND PRINT Vizualne komunikacije comTEC® d.o.o. Oblikovanje revije in znaka revije / JOURNAL AND LOGO DESIGN Andrej PREDIN Na naslovni sliki je prikazan model vertikalne vetrnice v Laboratoriju za Aero- in Hidroenergetske tehnologije (LAHET) na Fakulteti za energetiko. Cover picture shows a model of vertical wind turbine in the Laboratory for Aero-and Hydro-energy technologies (LAHET) at the Faculty of Energy. Revija JET je sofinancirana s strani Javne agencije za knjigo Republike Slovenije ter v okviru Petrolovega Programa velikih zavezancev za zagotavljanje prihrankov energije pri končnih odjemalcih / The Journal of Energy Technology is co-financed by the Slovenian Book Agency and Petrol's Programme for providing energy savings of end consumers. Mikro elektrarna - pot k zadostitvi električnih potreb naših domov Cene energentov, tudi električne energije, nezadržno naraščajo iz leta v leto. Bližajo se časi, ko bodo mikro elektrarne tudi ekonomsko zanimive. S kombinacijo npr. sončne in vetrne energije, bi lahko v veliki meri zadostili potrebe po električni energiji v naših domovih. S kompenzacijo lastnega odjema elektrike v naših domovih bi si lahko bistveno znižali strošek električne energije. Navkljub dejstvu, da se podpore RS proizvedeni energiji iz obnovljivih virov iz leta v leto nižajo, se investicija v mikro elektrarno izplača pred desetimi leti obratovanja. Podpora se torej še naprej niža, strošek oz. prispevek za rabo obnovljivih energetskih virov pa se na naših položnicah veča. Pred koncem lanskega leta, 28. 11. 2012, je vlada RS sprejela spremembe Uredbe o podporah električni energiji, proizvedeni iz obnovljivih virov energije. Za sončne elektrarne, ki bodo vključene v elektro omrežje RS po 1. 12. 2012, uredba določa 2 % mesečno zniževanje podpore. S 1. 1. 2013 pa na novo uvaja 5% bonus za tiste sončne elektrarne, katerih skupna moč ne presega 5 kW. Zahtevan je tudi t.i. priklop »za lastnim« števcem stavbe, ki omogoča lasten odjem pred oddajo proizvedene električne energije v omrežje. S tem ukrepom uredba vzpodbuja vključitev mikro sončnih elektrarn v interno omrežje končnega odjemalca - uporabnika, torej lastnika stavbe, kar je bilo do sedaj le izjemoma dovoljeno. S kombinacijo sonca in npr. majhne vetrnice bi lahko na ugodnih lokacijah razširili proizvodnjo tudi v nočni čas in tako pokrili lastne potrebe po električni energiji doma tudi v nočnem času. S takšno kombinacijo, na primer 3,5 kWp sončne in 1,5 kWp vetrne elektrarne, bi ostali znotraj 5 kW skupne moči. Investicijska ocena takšne kombinacije je okrog 11.300 €. V mesecu januarju je trenutno veljavna podpora zagotovljenega odkupa (ZO) 0,147 €/kWh in 5% bonus za moči do 5kW. Upoštevamo življenjsko dobo solarnih modulov - 25 let. Upoštevamo še 0,5% propad izkoristka sončnih modulov letno. Upoštevana življenjska doba vetrnice pa je 25 let brez padca izkoristka. Na koncu upoštevamo še prodajo električne energije po 15-tih letih proizvodnje po ceni 0,065 €/kWh in letno stopnjo rasti cene 1 %. Izračuni, upoštevajoč navedene podatke, kažejo, da je donosnost takšne naložbe bistveno višja (več kot dvakrat), kot na primer varčevanje sredstev v bankah. Doba vračanja take naložbe je med 7 do 11 let, odvisno od tega, kako smo se sposobni prilagoditi lastni porabi in s tem kompenzaciji lastnega odjema. To pomeni, da svojo porabo v največji možni meri prilagajamo v čas lastne proizvodnje. Tako v Evropski uniji, kot tudi v Sloveniji, prihajajo časi, ko bo kombinacija sončne in vetrne energije postala najcenejši in najbolj ekonomični vir električne energije v naših domovih. Zato ne bo več optimalno izkoriščati ves razpoložljiv prostor na lastni strehi, ampak le toliko, da zadostimo lastnim potrebam - do 5 kW instalirane moči, kar zadošča porabi povprečne 4 - članske družine. S takšno investicijo bi vsekakor razbremenili tudi okolje in hkrati vzpodbudili rodnost v Sloveniji, torej da bi spet prišli do povprečne 4 - članske družine. Bralce obveščamo, da se nam je pri tisku naše revije JET, št. 3/12 zgodila neljuba napaka, saj je bil pomotoma objavljen tekst članka z naslovom: Dež kot zeleni energetski vir?. Pravilni tekst objavljamo v tej številki - članek z naslovom: Generični model za kreiranje lopatic vetrnih turbin, avtorjev G. Hren, I. Žagar in A. Predin. Za neljubo napako se vam opravičujemo. Krško, februar 2013 Andrej PREDIN Micro-power plants - meeting the electrical needs of our homes Energy prices, including electricity, are rising steadily from year to year. The time when micro-power plants will become economically appealing is approaching. Through a combination of solar and wind power, they could essentially satisfy the demand for electricity in our homes. By offsetting the consumption of electricity in our homes, we could significantly reduce the cost of electricity. Despite the fact that the support of the Republic of Slovenia for renewable energy has been on a downward trend, an investment in micropower can be returned over ten years of operation. Support therefore continues to decrease, although the cost of contribution for the use of renewable energy resources has increased. Late last year, the Republic of Slovenia adopted amendments to its regulation on support for electricity produced from renewable energy sources. For solar power plants connected to the electrical grid, support has decreased by 2%. Starting in 2013, there is a new 5% bonus for solar power plants with a total power up to 5 kW of installed power. Furthermore, a so-called connection "after own" building meter, which enables the electrical energy produced to be consumed before being transmitted to electrical grid, is required. With this measure, connections of the micro-solar power plant systems to user's internal electricity grid are now encouraged, whereas they had previously rarely been allowed. The combination of solar power and small wind turbines (on favorable locations) could expand the household production of electricity throughout the entire day. For example, a 3.5 kWp solar and 1.5 kWp wind turbine could yield 5 kW of total power for an initial investment evaluation of around €11,300. In January, the current support guaranteed purchase price is € 0.147/kWh, and a 5% bonus for solar power plants up to 5kW of installed power. The estimated lifespan of solar modules is up to 25 years, with a 0.5% efficiency decrease per year. In the calculations, the estimated lifespan wind turbine is 25 years without falling efficiency. Finally, we take into account the sale of electricity after 15 years of the production up to €0.065/kWh and annual growth rates of 1%. From calculations tak ing into account the above data, it is shown that the return on such investment is substantially higher (more than twice) than the interest offered on deposits in banks. The return on investment period is between 7 and 11 years, depending on how we are able to best modify our consumption of electricity according to our production. A time is coming when a combination of solar and wind energy will become the cheapest and most economical source of electricity in our homes, in the European Union, as well as in Slovenia. Therefore, the optimal utilization is not all the available space on the our roof but only as much to meet our own needs - up to 5 kW of installed power capacity is enough to cover the needs of the average four-member family. Such an investment would certainly provide relief to the environment. Announcement to Readers: In our print magazine JET, no. 5-3, there was an error, showing an article entitled Rain as a green energy source? The correct text is, in fact, published in this issue, now entitled Generic model for wind turbine blades design by G. Hren, I. Žagar, and A. Predin. We apologize for this unfortunate error and any trouble it may have caused. Krško, February, 2013 Andrej PREDIN Table of Contents / Kazalo Design of electrical machines by using conformal mapping / Konstruiranje električnih strojev z uporabo konformnih preslikav Jan Šlamberger, Peter Virtič........................................................................................................13 Impact of the characteristics of overhead ground wires on the current reduction factor, and their effect on the change of grounding system potential / Vpliv lastnosti nadzemnega žičnega ozemljila na koeficient znižanja toka in posledično na spremembo potenciala ozemljitvenega sistema Mario Havranek, Tomislav Strinic, Goran Kneževic....................................................................19 The unsteady static-stall aerodynamic characteristics of an S809 airfoil at low Reynolds numbers / Aerodinamične karakteristike nestacionarnega zastoja na profilu S809 pri nizkih Reynoldsovih številih Matej Fike, Gorazd Bombek, Matjaž Hriberšek, Aleš Hribernik.................................................33 Fracture evaluation of energy components with local brittle zones / Vrednotenje loma v energetskih komponentah z lokalnimi krhkimi področji Zdravko Praunseis........................................................................................................................51 Generic model of wind turbine blades / Generični model lopatic vetrne turbine Gorazd Hren, Andrej Predin, Ivan Žagar......................................................................................61 Instructions for authors...............................................................................................................69 JET Volume 6 (2013), p.p. 13 - 18 Issue 1, February 2013 http://www.fe.um.si/en/jet/e-jet.html DESIGN OF ELECTRICAL MACHINES BY USING CONFORMAL MAPPING KONSTRUIRANJE ELEKTRIČNIH STROJEV Z UPORABO KONFORMNIH PRESLIKAV Jan Šlambergerw, Peter Virtič Keywords: Electrical machine, Conformal mapping, Slot opening, Schwarz-Christoffel transformation Abstract The design of electrical machines requires good working knowledge of magnetic fields in air gaps, which is very difficult or analytically unsolvable in most cases. With the development of computers, numerical methods came to the fore, enabling very good approximations of real values to be calculated. One of the major disadvantages of numerical methods is the duration of the calculations, particularly with the construction of prototypes, in which the structure changes, thus requiring more calculations. With the goal of bringing about quicker calculations, analytical methods were used, and combinations of analytical methods with numerical methods were re-started. This article will present an analytical calculation of magnetic fields using conformal mappings. Povzetek Načrtovanje električnih strojev med drugim zahteva dobro poznavanje megnetnega polja v zračni reži, kar pa je v večini primerov zelo zahtevno oziroma analitično nerešiljivo. Z razvojem računalnikov so prišle v ospredje numerične metode, s katerimi lahko izračunamo zelo dobre približke realnim vrednostim. Ena večjih slabosti numeričnih metod so dolgotrajni izračuni, še posebaj pri konstruiranju prototipov, kjer se konstrukcija spreminja in s tem potrebujemo več izračunov. Z željo, po čim hitrejših izračunih so se ponovno začele uporabljati analitične metode m Corresponding author: Asst. Jan Šlamberger B.Sc.E.E., Mailing address: University of Maribor, Faculty of Energy Technology, Hočevarjev trg 1, SI-8270 Krško Tel.: +386 3 777 0407, Fax: +386 3 777 0413, E-mail address: jan.slamberger@uni-mb.si jet Journal of Energy Technology Jan Šlamberger, Peter Virtič JET Vol. 6 (2013) Issue 1 ter kombinacije analitičnih metod z numeričnimi metodami. V tem članku bo predstavljen analitičen izračun magnetnega polja s pomočjo konformnih preslikav. One of the most important necessities in the design of electrical machines is knowledge of magnetic fields in air gaps. There are several types of approaches to solving the issues of magnetic fields in air gaps; one such approach is analytical methods, a combination of analytical and numerical methods. In this paper, analytical conformal mapping is presented, which connects the symmetrical slotted air gap to the slot-less air gap. The connection between the two gaps was given by Zarko, Ban and Lipo, [1], and Gibbs, [2]. The slot opening of a symmetrical slotted air gap in the Z and W planes was given by Markovic, Jufer and Perriardin, [3], and Zhu and Howe, [4]. There are four conformal transformations necessary to transform the slotted air gap into a slot-less air gap. A single slot of the original geometry is shown in Fig.1. This geometric shape needs to be transformed into its linear model in the Z plane, shown in Fig. 2, using a logarithmic conformal transformation defined as where s = m + jn = re]e, z = x + jy. The link between the coordinates in the S and Z plane is 1 INTRODUCTION 2 SLOT OPENING z = ln(s), (2.1) x = 0 (2.2) and (2.3) 5 4 e e 2 ^ Figure 1: Slot opening in the S plane w=-1 4 w=-a y W=o> w=a w=1 3 Figure 2: Slot opening in the Z plane The coefficients b0 and g are defined as b0 - 2-e * - ln|R (2.4) The second transformation is to transform the geometric structure in the Z plane into the upper half of the W plane, using a Schwarz-Christoffel (SC) transformation, shown in Fig. 3. In the symmetrical slotted air gap, the SC transformation will have the form -b fOpEfc^E dw ( w-1)( w +1) (2.5) The unknown coefficient a, which represents the values of w at the corner points, is defined as [4] a - „ (1 + V b0 J (2.6) 6 4,5 2,3 1 w=-a w=-1 w=1 w=a Figure 3: Slot opening in the Wplane The next transformation is required from the T plane where the field is regular to the W plane. The slot opening in the T plane represents two parallel plates extending an infinite distance in both directions, as shown in Fig. 4. The transformation from the 7" plane into the W plane is given by b 5 2 g x v u b r t =— f jn (w-i)(w +1) dw (2.7) w=-1 q 2 T g W=x T W=1 i- -i,-► Figure 4: Slot opening in the T plane The transformation of linear geometry in the T plane into curved geometry in the K plane (Fig. 5) requires an exponential function in the form k=e. (2.8) 2 Figure 5: Slot opening in the K plane 3 FIELD SOLUTION IN THE SLOTTED AIR GAP The field solution in the K plane, which represents a slot-less air gap, can now be mapped back to the S plane. The connection between magnetic field in the K and the S plane is, [1], = Bk if I , os (3.1) where (dk/ds) conjugate value of (dk/ds) is. The partial derivate (dk/ds) can be expressed as 8k Ok 0t dw Oz p Os 8t dw 8z 8s 4 3 P The partial derivatives in (3.2) are defined by conformal transformations between the corresponding complex planes 3k t In k , — = e = e = k dt dt b 1 dw jn w2 -1 dw jn w2 -1 2 , (3.3) dz b ddL-1 ds s Considering (3.2) and (3.3) in (3.1) yields Bs = Bk fk 1 ^ V s Vw2 - a2 J (3.4) The variables k and s can both be expressed as a function of w. Combining (2.1) and (2.5) yields (3.5) b ( . -1 ( w^ va2 -1, va2 -w2 +w\$a2 -1 ^ I sin 1 — +--ln nl Va J 2 ^ä^-wyfa2-. s = e Combining (2.7) and (2.8) yields jb, |w+1| ln- k = e2n |w-1. 4 CONCLUSION (3.6) This paper presents the analytical conformal mapping that connects a slotted air gap with a slot-less air gap. With this mapping, the magnetic field in slotted air gap can be transformed into a slot-less air gap, be solved and then be mapped back to a slotted air gap. This method is very useful in the early design stages of electrical machines, because it is much faster than numerical methods. References [1] D. Žarko, D. Ban, T.A. Lipo: Analytical Calculation on Magnetic Field Distribution in the Slotted Air Gap of a Surface Permanent-Magnet Motor Using Complex Relative Air-Gap Permeance, IEEE Trans. Magn., vol. 42, no. 7, pp. 1828-1837, 2006 [2] W.J. Gibbs: Conformal Transformations in Electrical Engineering, Chapman & Hall, 1958 [3] M. Markovic, M. Jufer, Y. Perriard: Determination of Tooth Cogging Force in a Hard-Disk Brushless DC Motor, IEEE Trans. Magn., vol. 41, no. 12, pp. 4421-4426, 2005 Jan Šlamberger, Peter Virtič JET Vol. 6 (2013) Issue 1 [4] Z. Q. Zhu, D. Howe: Instantaneous Magnetic Field Distribution in Brushless Permanent Magnet dc Motors, Part III: Effect of Stator Slotting, IEEE Trans. Magn., vol. 29, no.1, 1993 JET Volume 6 (2013), p.p. 19 - 32 Issue 1, February 2013 http://www.fe.um.si/en/jet.html IMPACT OF THE CHARACTERISTICS OF OVERHEAD GROUND WIRES ON THE CURRENT REDUCTION FACTOR, AND THEIR EFFECT ON THE CHANGE OF GROUNDING SYSTEM POTENTIAL VPLIV LASTNOSTI NADZEMNEGA ŽIČNEGA OZEMLJILA NA KOEFICIENT ZNIŽANJA TOKA IN POSLEDIČNO NA SPREMEMBO POTENCIALA OZEMLJITVENEGA SISTEMA Mario Havranek, Tomislav Strinic, Goran KneževicM Keywords: surface potential, current reduction factor, grounding, touch voltage, step voltage Abstract This paper shows the impact of the characteristics of overhead ground wires on the change of grounding system potential in a power system substation during a single line to ground fault. The touch and step voltages can be affected by changing the overhead ground wire radius. The calculation of the current reduction factor depending on the size of the overhead ground wire is also presented. Furthermore, simulations of the grounding system potential, the surface potential as well as the touch and step voltages of a study case have been made using CYMGRD simulation software, which has implemented the IEEE 80-2000 standard (IEEE Guide for Safety in AC Substation Grounding). The power system substation observed in the simulation includes two incoming 110 kV transmission lines with two 110 kV bus bars as well as 35 kV and 10 kV bus bars. The entire grounding system of the substation is comprised of a conducting net aligned horizontally, and mutually interconnected by vertically placed conducting rods. The simulations M Corresponding author: Goran Kneževic, MEng, University of Osijek, Faculty of Electrical engineering, Tel.: +385 91 224 6126, Fax: +385 31 224605, Mailing address: Kneza Trpimira 2b, 31000 Osijek, Croatia, E-mail address: goran.knezevic@etfos.hr Journal of Energy Technology performed show the change of the observed potential values for a given size of the overhead ground wire. In order to determine the maximum fault current that can occur in the observed substation (which is needed for calculation of the substation ground potential), a single-phase short circuit is simulated on the 110 kV bus bars, using the simulation software DIgSILENT Power Factory. Povzetek V tem članku je predstavljen vpliv lastnosti nadzemnega žičnega ozemljila na spremembo potenciala ozemljitvenega sistema v razdelilni transformatorski postaji v primeru enofaznega kratkega stika z zemljo. Sprememba polmera nadzemnega žičnega ozemljila lahko povzroči spremembo napetosti dotika in napetosti koraka. Prav tako je predstavljen tudi izračun koeficienta znižanja toka, ki je odvisen od dolžine nadzemnega žičnega ozemljila. Predstavljene pa so tudi simulacije potenciala ozemljitvenega sistema, površinskega potenciala, napetosti dotika in napetosti koraka dotičnega primera obravnave, ob uporabi CYMGRD simulacijskega programskega orodja, zasnovanega na standardu IEEE 80-2000 (IEEE Navodila za varnost pri ozemljevanju izmeničnih razdelilnih postaj). Obravnavana razdelilna transformatorska postaja vključuje dva 110 kV dovodna daljnovoda z dvema 110 kV zbiralkama ter 35 kV in 10 kV zbiralkama. Pri tem je celotni ozemljitveni sistem razdelilne postaje sestavljen iz horizontalno položene prevodne ozemljitvene mreže, ki je medsebojno povezana z navpičnimi prevodnimi ozemljitvenimi palicami. Izvedene simulacije so pokazale spremembo opazovanih vrednosti potenciala pri določeni dolžini nadzemnega žičnega ozemljila. Z namenom določitve maksimalne vrednosti kratkostičnega toka, ki lahko steče v obravnavani razdelilni postaji in je potreben za izračun ozemljitvenega potenciala, je z uporabo simulacijskega programskega orodja DIgSILENT Power Factory izvedena simulacija enofaznega kratkega stika na 110 kV zbiralkah. 1 INTRODUCTION When a ground fault occurs in a substation, the flow of the fault current to the earth produces potential gradients within and around the grounding system. In such cases, the 'step voltage' is defined as the difference between the surface voltage experienced by a person whose stride is 1 m long, [1]. The surface voltage where individuals stand while their hands are in contact with a grounded structure is called a 'touch voltage', [1]. Good design of a substation grounding system should have low grounding system resistance with considerable touch and step voltages that are tolerable to a person, [2]. For the grounding system resistance of large power substations, in the majority of cases, their resistance is about 1 Q or even lower, while the grounding system resistance of smaller power substations is usually between 1 and 5 Q, [3]. It is challenging to estimate the resistance of a grounding system based on the length of its conductors, especially if conducting rods are used. The effort for calculation is greatly eased by using specialized simulation software tools, mostly available commercially, based on numerical methods for the calculation of the desired variables. CYMGRD is simulation software specialized for grounding systems that have implemented the IEEE 80-2000 standard, [4]. Based on a plotted grounding system, CYMGRD calculates the allowed touch and step voltages for a given single-phase short circuit current. In this paper, the impact of the overhead ground wire (OHGW) characteristics on the change of a power system substation ground potential is presented. Simulation cases are based on a real power system substation of 110/35/10 kV with its corresponding grounding system. Determination of single- phase short circuit current, needed for grounding system analysis, is performed using DIgSILENT Power Factory software, [5]. 2 DETERMINATION OF THE CURRENT RESPONSIBLE FOR THE INCREASE OF GROUNDING SYSTEM POTENTIAL The total amount of the single-phase short circuit current If can be presented as the sum of the transformer neutral currents ITR and the sum of the short circuit currents of each transmission line 3I0 , [6]. The equivalent circuit diagram of a single-phase short circuit current distribution in case of the fault in power system substation with two incoming transmission lines is shown in Figure 1. B B f/g- — Sf 1 3I01 + Sf 2 31o Figure 1: The short circuit current distribution of a power system substation with two incoming transmission lines (A - phase, B - neutral) In order to determine the current responsible for the increase of the grounding system potential, the amount of the current that flows through the grounding system needs to be known. In case of a short circuit inside the power system substation, the fault current is divided between the local ground grid and several metallic return paths associated with the incoming and outgoing lines, i.e. overhead ground wires, cable sheaths and feeder neutrals, [7]. The IEEE 80-2000 standard's Guide for Safety in AC Substation Grounding describes a way of calculating the reduction factor based on the information of incoming transmission lines and power system substation bus bars, [8]. It describes the procedure for the calculation of the current that goes through the grounding system into the ground, and gives the data needed for the calculation of the root mean square value of the asymmetrical failure current for the grid value of XjR in respect to the transient state time. The quotient of the current which goes through the grounding system Ig and the short circuit current 3I0 gives the reduction factor Sf given in the form of Equation (2.1): Sf = A. (2.1) f 3Io where: Sf is the current reduction factor; Ig is the value of the current which flows through the grounding system; 3I0 is the short circuit current through the neutral conductor. The current that flows through the grounding system needs to be corrected by the factor Df as shown by Equation (2.2): Ig = D fig (2.2) where: Ig is the maximum current which can appear in the grounding system; Df is the transient state factor. Expressing the symmetrical value of the grounding network current Ig by rearranging the Equation (2.1); substituting it into Equation (2.2) yields Equation (2.3) for the maximum grounding system current: Ig = DfSf 3Io (2.3) The Df factor takes into account the network value of XR during the transient state. In case of a failure time larger than 30 periods for the nominal current harmonic of 50 Hz, the value of the factor Df is approximately equal to 1, as presented in Figure 2. 1JS0 1 40 I l.ä) s j IjM \ 0.80 I o.eo D. 40 0.20 0.00 0.03 OOS 0.0? 0.CÖ 0.J Oti 0.3 0.4 O.S S-Vsult duration [n] Figure 2: Decrement factor in relation to transient state time The current reduction factor Sf depends on the overhead protection wire resistance; it can be calculated using Equation (2.1), giving Equation (2.4), [9]: / 7 Sf = = 1 - (2.4) f 3/o Ze where: 7ie is the mutual impedance of the loop overhead ground wire and the line conductor with common earth return; ZE is the loop impedance of the overhead ground wire and the earth return. In order to determine previously mention variables following equations will be used, [9]: (2.5) Ze=R + m»+j -E 8 J mo» fln f S1 + m Zie = + j S = 2n I I r I 4 M0» Jn f S 2n ^ dj, 1.85 (2.6) (2.7) im» , Pe where: 7ie is the mutual impedance per unit length of the loop overhead ground wire and the line conductor with common earth return; ZE is the loop impedance per unit length of the overhead ground wire and the earth return; dIE is the distance between the conductors and the overhead protection wire; S is the depth of the earth path; pE is the soil resistivity; R is the resistance per unit length of the overhead protection wire. 3 INFLUENCE OF THE FAULT PROTECTION SYSTEM SPEED ON THE ALLOWED CURRENT VALUE THROUGH A HUMAN BODY The current value that a person can feel when in contact with overvoltage is relatively low: about 1 mA. The amount of current from 1 mA to 6 mA is called the 'releasing current'; by definition, it represents the current value that does not cause muscle contraction. For women, its maximum value is taken to be 6 mA, while for men its maximum value is taken to be 9 mA. The values between 60 mA to 100 mA represent a danger of ventricular fibrillation and the halting of heartbeat. Non-fibrillation current of amplitude /B with the duration time between 0.03 s and 3 s is related to absorbed energy that the body receives; it is described by Equation (3.1): where: SB = Ilts (3.1) IB is the root mean square value of the current which flows through the body; ts is duration time of failure; SB is the empirical constant related to allowed absorbed energy of the body. The amplitude and the duration time of the current that flows through the human body needs to be lower than the one that causes ventricular fibrillation. The duration of the current flow through the human body that can be tolerated in the majority of people is directly related to the current amplitude, [8]. It is assumed that 99.5% of people weighing about 50 kg can withstand the current with amplitude IS and duration time ts represented with the equation (3.2): Ib = ■ 0.116 For people with weight from about 70 kg upwards, Equation (3.3) is usually used: 0.157 Ib =■ VtT (3.2) (3.3) Equations (2.6) and (2.7) are both valid if the duration time of failure is between 0.03 s and 3 s; it is not suitable for very short or a very long duration failure time. In Figure 3, the allowed current values in respect to the shock duration time with values from 0.03 s to 3 s are presented. Figure 3: The value of the allowed current through the body Impact of the characteristics of overhead ground wires on the current reduction factor, and their effect on the change of grounding system potential 4 GROUNDING SYSTEM SIMULATION The power system substation observed in the simulation is the Valpovo transformer station, which includes two incoming 110 kV transmission lines with two 110 kV bus bars as well as 35 kV and 10 kV bus bars. The grounding grid is placed 0.8 m underground and made up of 12 transverse and 9 longitudinal conductors mutually interconnected creating the grid of 110 m x 80 m with a total conductor length of 1950 m. In order to lower the grounding system resistance, 25 grounding 3 m long conducting rods are connected to the grounding grid. The grounding system is presented in Figure 4. Figure 4: The Grounding system layout 4.1 Soil Analysis For simulation of the grounding system, it is necessary to have the soil resistivity data, which in this case has been obtained by measuring using Wenner method. The measured values were imported into the CYMGRD simulation software. The two-layer soil model simulation has been performed. The results of the soil analysis are presented in Table 1. Table 1: Results of the soil analysis Upper Layer Thickness (m) 9.02 Upper Layer Resistivity (Qm) 72.34 Lower Layer Resistivity (Qm) 11.96 4.2 Safety Parameter Analysis The allowed touch and step voltages depend on different parameters, such as shock duration, body weight, surface layer thickness and material. With a surface layer thickness of 0.1 m and resistivity of 500 Qm for a body weight of 70 kg, the allowed touch and step voltages are calculated depending on the shock duration. The simulation results obtained using CYMGRD software according to IEEE 80-2000 standard are presented in Figure 5. 0.1 0.2 0.3 0.4 Shock duration (s) 0.5 0.6 0.7 Maximum Permissible Touch Maximum Permissible Step Figure 5: Maximum permissible values of touch and step voltages in respect to shock duration time 4.3 Single-phase short circuit current determination In order to determine the increase in the grounding system potential, it is necessary to determine the single-phase short circuit current; in this case, it has been done using DIgSILENT Power Factory simulation software. The simulation model covers the entire Croatian transmission system. The results are shown in Figure 6. Impact of the characteristics of overhead ground wires on the current reduction factor, and their effect on the change of grounding system potential $ DOr" 6.00-- 4jOO -...... 2 00 O.OO- -2M- -0.1000 0.03B4 - M It a Oil r Q.1767 0.3151 0.4534 [t] 0.581 Figure 6: Single-phase failure duration time Disconnection time of the short circuit protection system is set to 0.5 s. The value of the singlephase short circuit current equals 7480 A. 4.4 Grounding System Analysis The current reduction factor is calculated using Equations (4-7) considering two incoming transmission lines for different radii of OHGW, which are some of the standard sizes of zinc-coated steel OWHG that could be used in the observed transmission lines. The geometric mean distance between the line conductors and OWHG is 5.76 m, while the soil resistivity is assumed to be 72.34 Qm. The results are presented in Fig. 7. 0.9 0.85 0.8 o 0.75 ■M £ 0.7 C o 0.65 Č 0.6 ■c tU t£. 0.55 0.5 0.45 0.4 m 35 50 70 95 Overhead ground wire size (mm2) 120 Figure 7: Current reduction factor in respect to some standard sizes of overhead ground wires For the several standard sizes of OWHG, the grounding system potential, the maximum surface potential, the minimum surface potential and the maximum touch voltage are calculated using CYMGRD. The fault duration is set to 0.5 s. The transient state factor Df is assumed to be 1. The value of the single-phase short circuit current is 7480 A. The results are presented in Table 2. Table 2: Results of the grounding analysis for different sizes of the overhead protection wire OHGW size (mm2) Current reduction factor Grounding system potential (V) Max. surface potential (V) Min. surface potential (V) Max. touch voltage (V) 35 0.84 876.4 784.6 515.6 360.8 50 0.78 813.7 728.5 478.7 335 70 0.72 751.1 672.5 441.9 309.2 95 0.64 667.7 597.8 392.8 274.9 120 0.56 584.2 523.1 343.7 240.5 As it can be seen in Table 2, the grounding system potential and the surface potential decrease with the increase of OHGW size. Due to the results of the safety parameter analysis presented in Chapter 4.2, for the fault duration of 0.5 s and the body weight of 70 kg, the maximum permissible touch voltage is 344.34 V. In the case in which 35 mm2 OHGW is used, the maximum touch voltage that occurred exceeds the maximum permissible touch voltage. The distributions of the touch voltage and the surface potential around the grounding grid for the case of 35 mm2 OHGW are presented in Figure 8 and Figure 9, respectively. Figure 8: The distribution of the touch voltage around the grounding grid for the case with OHGW size of 35 mm2 Figure 9: The distribution of the surface potential around the grounding grid for the case with OHGWsize of 35 mm2 A profile plot along the main diagonal of the grid using a stride of 1 m is presented in Figure 10. The maximum step voltage that occurred equals 49.1 V, which is less than the maximum permissible step voltage. Figure 10: Profile plot along the main diagonal of the grid using a step of 1 m for the case with OHGW size of 35 mm2 The size of OHGW that is actually used on 110 kV transmission lines that are connected to the observed transformer station is 70 mm2. In this case, the maximum permissible touch and step voltages are not exceeded. Distributions of the touch voltage and the surface potential around the grounding grid for the case with 70 mm2 OHGW are presented in Fig. 11. and Fig. 12. A profile plot along the main diagonal of the grid using a stride of 1 m for the case with OHGW size of 70 mm2 is presented in Fig. 13. Maximum occurred step voltage equals 39.5 V. «0 isa ng 4 ! 1« 111 1H IIS 1« Langlh (mstar*) Surface Pudendals -Step PotefVlials Touch TotenlMis Ground Paential Rise ■ 751 1SS votts Meimiiutti Pt-TPinssiblc Silt* - 71Mawmum Pern»«**? Touch - >H.3$vott$ Figure 13: Profile plot along the main diagonal of the grid using a step of 1 m for the case with OHGW size of 70 mm2 4 CONCLUSION This paper shows the impact of the characteristics of overhead ground wires on the change of the grounding system potential in the observed power system substation during a single line-to-ground fault. Calculation of the current reduction factor depending on the size of the overhead ground wire is presented. Simulations of the grounding system potential, the surface potential, and the touch and step voltages have been made using CYMGRD simulation software for a particular study case. The simulation results show that the touch and step voltages, as well as the surface potential decrease with the increase of the overhead protection wire radius. Permissible touch and step voltages are exceeded if 35 mm2 overhead ground wire is used in the observed substation. The simulation results for the case of the mm2 overhead ground wire, which is actually used on the transmission lines connected in the observed transformer station, show that permissible touch and step voltages are not exceeded. References [1] L.H. Chena, J.F. Chena, T.J. Liang, W.I. Wang: Calculation of ground resistance and step voltage for buried ground rod with insulation lead, Electric Power System Research 2008, Vol.78., pp. 995-1007. [2] D. Prasad, H.C. Sharma: Significance of Step and Touch Voltages, International Journal of Soft Computing and Engineering 2011, Vol.1(5), pp. 193-197. [3] J. Ma, F. P. Dawalibi: Modern Computational Methods for the Design and Analysis of Power System Grounding, Proceedings of International Conference on Power System Technology, p.p. 122-126, 18-21 August 1998. [4] http:// www.cyme.com [5] http:// www.digsilent.de [6] Franjo Majdandžic: Groundings and grounding systems, Graphis, pp. 133-143, 2004. [7] P.L. Buccheri, S. Mangione: Analysis of ground fault current distribution along nonuniform multi-section lines, Electric Power System Research 2008, Vol. 78(9), pp. 1610-1618. [8] IEEE: IEEE Guide for Safety in AC Substation Grounding, The Institute of Electrical and Electronics Engineers, p.p. 1-192, 2000. [9] J. Schlabbach: Short circuit current, The Institution of Engineering and Technology, pp. 172-179, 2005. Nomenclature (Symbols) (Symbol meaning) Df transient state factor Sf current reduction factor Ig the value of the current which flows in grounding system 3I0 short circuit current through the neutral conductor IG maximum current that can appear in the grounding system IB root mean square value of current that flows through the body ts fault duration time SB empirical constant related to allowed absorbed energy of the body mutual impedance of the loop overhead ground wire and line conductor with common earth return mutual impedance per unit length of the loop overhead ground wire and line conductor with common earth return ZLE ZE loop impedance of overhead ground wire and earth return zle ZE loop impedance per unit length of overhead ground wire and earth return dLE distance between conductors and overhead protection wire S depth of the earth path pE soil resistivity R resistance per unit length of the overhead protection wire 'J Technology ^ Journal of Energy JET Volume 6 (2013), p.p. 33 - 50 Issue 1, February 2013 http://www.fe.um.si/en/jet/e-jet.html THE UNSTEADY STATIC-STALL AERODYNAMIC CHARACTERISTICS OF AN S809 AIRFOIL AT LOW REYNOLDS NUMBERS AERODINAMIČNE KARAKTERISTIKE NESTACIONARNEGA ZASTOJA NA PROFILU S809 PRI NIZKIH REYNOLDSOVIH An investigation was conducted to study the unsteady static-stall characteristics of an S809 airfoil whose aerodynamic characteristics are representative of a horizontal axis wind-turbine. It is very difficult to experimentally investigate this phenomenon, especially when attempting to study the development of flow-separation when the blades are rotating. The application of wind-tunnel tests or CFD simulation is obviously more appropriate. In order to investigate unsteady static-stall regarding the aerodynamic characteristics of an S809 airfoil, a comparative study was performed to obtain numerical and experimental results for the flow-separation position, airfoil pressure distribution, and velocity profiling that included velocity oscillations. The experimental results were acquired using PIV (Particle Image Velocimetry) and the study was performed at various angles of attack (AOA). No separation was observed at low AOA, but at 9.6° AOA the separation vortex comprised 50% of the airfoil's chord length, whilst a complete stalling of the airfoil occurred at 20° AOA. The observed separation zone was not steady but was found to oscillate around its mean-position at an interval of ± 10% of the chord's length. Neither of the applied turbulence models k-e nor SST used in 2D unsteady m Corresponding author: Matej Fike, Faculty of Mechanical Engineering, University of Maribor, E-mail address: matej.fike@gmail.com ŠTEVILIH Matej Fikew, Gorazd Bombek, Matjaž Hriberšek, Aleš Hribernik Keywords: S809 airfoil, low Reynolds number application, unsteady static-stall Abstract simulation predicted this oscillation, although the numerical results agreed fairly well with those experimentally obtained, especially the averaged velocity and vorticity fields around the suction side of the airfoil when using the SST model. Povzetek Opravljena je bila raziskava, kjer smo preučevali karakteristike nestacionanega statičnega zastoja na profilu S809, katerega aerodinamične karakteristike so reprezentativne za vetrne turbine z horizontalno osjo. Eksperimentalno je zelo težko preučevati ta pojav, še posebej kadar želimo preučevati odcepljanje toka na rotirajočih lopaticah. Izvedba meritev v zračnem tunelu oziroma simulacij računalniške dinamike tekočin je primernejše. Z namenom raziskovanja nestacionarnega statičnega zastoja glede aerodinamičnih karakteristik S809 profila je bila opravljena primerjalna študija, tako so bili primerjani eksperimentalni in numerični rezultati položaja odcepitve toka, porazdelitve tlaka po profilu in hitrostih profilov vključno z oscilacijami hitrosti. Eksperimentalni rezultati so bili pridobljeni z uporabo PIV (meritev hitrosti z odslikavo delcev) in pri različnih napadnih kotih. Pri nizkih napadnih kotih nismo zaznali odcepitve toka, pri napadnem kotu 9.6° je odcepljen vrtinec zavzemal 50% dolžine tetive profila, pri napadnem kotu 20° zastojni vrtinec obsega ves profil. Opazovana področje odcepitve ni bilo stacionarno ampak smo ugotovili, da oscilira okrog njene srednje pozicije v intervalu ±10% dolžine tetive. Uporabljena turbulentna modela, k-£ in SST, v 2D časovno odvisnih simulacijah nista napovedala teh oscilacij, čeprav so se numerični rezultati dokaj dobro ujemali z eksperimentalno pridobljenimi, še posebej polja povprečne hitrosti in polja vrtinčnosti na sesalni strani profila pri uporabi SST turbulentnega modela. 1 INTRODUCTION Wind-turbines operate within a broad span of wind-speed by starting their operations at 4 m/s and operating up to 24 m/s, this being the usual shut-down wind speed. A turbine's rotational speed is usually constant and the blade-pitch angle fairly rigid. However, if the relative windangle increases, this causes an increase in the angle of attack that may lead to flow-separation at the suction sides of the turbine's blades, and eventually to blade-stalling. Although flow-separation should be avoided in aviation unless firm-braking is intended, it is commonly used for wind-turbines. These are designed with a maximum power level that can be reached dozens of times per year. Predicting peak rotor power and post-peak power is important when designing constant-speed and variable-speed stall-regulated rotors. Thus, information about the state of the boundary layer over the suction surface of the airfoil is necessary when designing wind-turbine blades, in order to predict and consider the flow characteristics within the transitional regime. Furthermore, predicting the location of the transition-point is important for those wind turbine applications that operate at low Reynolds numbers. This is because the laminar/turbulent properties of the flow-field have an important influence on the skin-friction and separation that significantly affect both the lift and drag characteristics of the blade. Since the introduction of the boundary-layer concept, scientists and engineers have been facing a constant challenge to minimize its adverse effects, and then use it to their advantage. The present research was conducted to study the unsteady static-stall characteristics of an S809 airfoil, the aerodynamic characteristics of which are representative of a horizontal axis windturbine. This was in order to provide a greater understanding of those unsteady flow-separation processes that are effectively used for the stall-regulation of a wind-turbine's rotor. Low Reynolds number aerodynamics are characteristic for wind turbine applications in which low inflow velocities result in operational regimes with low-to-moderate wing-chord Reynolds numbers (i.e. chord Reynolds numbers ranging from 10,000 to 500,000). Whilst conventional aeronautic design principles usually either neglect the viscosity effect or restrict its influence to a thin region near the airfoil's body at high Reynolds numbers, a predominance of the fluid -viscosity effect for the low Reynolds number applications would result in boundary layers growing rapidly and separating from the surfaces of the airfoils easily, which should never be neglected. The behaviour of the laminar boundary layer on low-Reynolds number airfoils would significantly affect the aerodynamic performance of the airfoils. Since laminar boundary layers are unable to withstand any significant adverse pressure gradient, laminar flow separation is usually found on low-Reynolds-number airfoils. It was suggested by Lissaman, [1], and Mueller, [2], that the separated laminar boundary layers around low-Reynolds-number airfoils would behave more like free shear layers, which are highly unstable and, therefore, the rolling-up of Kelvin-Helmholtz vortex structures and transition to turbulent flows would be readily realised. When the adverse pressure gradient over the airfoil surface becomes large enough, the transition of the separated laminar boundary layer to turbulent flow could be conducted rapidly, and the increased entrainment of the turbulent flow could cause the turbulent flows to reattach to the airfoil's surface as a turbulent boundary layer. This would form what is termed a 'laminar separation bubble'. The reattached turbulent boundary could stay firmly attached to the airfoil surface up to the airfoil trailing edge. As the adverse pressure gradient increased with any increasing angle of attack, a second separation (separation of turbulent boundary layer) could occur. The second separation point would firstly move gradually along the suction surface with an increasing angle of attack, and then, starting at a certain AOA when the adverse pressure gradient had become more severe, it would suddenly (almost instantaneously) jump to the leading edge. This would be because of the sudden bursting of the laminar separation bubble, which would then result in airfoil stall. It has been reported that the position of the separation point in the case of a turbulent boundary layer is not static, [3]. Furthermore, Nishimura and Taniike, [4], showed a correlation between the frequencies of the separation point-position fluctuations and von Karman instability, which subsequently causes fluctuations of the lift and drag coefficient. However, common practice is to use the average values for lift and drag coefficient obtained by steady-state simulations, [5-6], or experimentally measured, [7-8]. The purpose of the presented research was to evaluate the oscillation zone of the turbulent layer's separation point, using a combination of flow visualization by the PIV method and numerical simulation done with CFD code, and to provide a greater understanding of those mechanisms responsible for flow-separation and its influence on lift and drag regarding the S809 airfoil. 2 MATERIALS AND METHODS 2.1 Airfoil section For this study, an airfoil was chosen whose aerodynamic characteristics were representative of horizontal axis wind-turbine (HAWT) airfoils: the S809. It is a 21% thick, laminar-flow airfoil designed specifically for HAWT applications. A sketch of the airfoil is shown in Figure 1. As reported by Somers, [8], at positive angles of attack below approximately 5°, the flow remains laminar over the forward half of the airfoil. It then undergoes laminar separation followed by a turbulent reattachment. As the angle of attack is increased further, the upper-surface transition point moves forward, and the airfoil begins to experience small amounts of turbulent trailing edge separation. At approximately 9°, the last 5 % to 10% of the upper surface is separated. The upper-surface transition point has moved forward to approximately the leading edge. As the angle of attack increases to 15°, the separated region moves forward to about the mid-chord. With further increases in the angle of attack, the separation moves rapidly forward to the vicinity of the leading edge, so that at about 20°, most of the upper surface is stalled. According to this, the presented investigations were focused on four different angles of attack: 5.1°, 9.6°, 15.1°, and 20.0°, respectively. 2.2 Wind tunnel The experiments were performed within a small wind-tunnel (cross-section 200 x 200 mm). The maximum air velocity in the tunnel was 36.3 m/s, and the turbulence intensity, measured by hot wire anemometry, was 1.5%. The size of the airfoil was designed according to the tunnel's dimensions and the desired pressure taps for pressure measurement. The airfoil was produced using rapid prototyping, and the chord's length (c) was selected as 100 mm, resulting in Reynolds number Re = 2 x 105 at maximum air-speed. A transparent extension was added to the wind-tunnel to enable PIV measurements. This extension and the airfoil fixture enabled a setting of the airfoil's angle of attack. The airfoil and the locations of the pressure measurement points are presented in Fig 1. _«o_ _ 40 SO Ji <0 :co Figure 1: Geometry of the airfoil and the pressure-tap's locations Sixteen positions for pressure measurements were applied, and 50 mm-long pneumatic tubes were used for connecting with the sensors. The sinusoidal response of the plastic tubing and sensors was, therefore, satisfactory for frequencies up to 300 Hz, since the gain factor was less than 1.05, [9]. Therefore, the frequency response of the system was satisfactory because the phenomena investigated, in particular the Von Karman vortex shedding, had a characteristic frequency below 300 Hz. GMSD 2.5 MR and GMSD 25 MR pressure sensors were used, with different measurement uncertainties of 0.5 Pa and 5 Pa, respectively. Pressure signals within durations of 10 seconds were acquired with 1 kHz sampling frequency. 2.3 PIV system The Dantec PIV (Particle Image Velocimetry) system was used to capture the flow-field around the airfoil (Fig. 2). The camera and laser were placed on a lightweight traverse system originally used for LDA measurements. A measurement plane was placed in the middle of the airfoil to reduce the influence of the wall. A two-cavity Nd: YAG laser was used, operating at high power with 50 mJ pulse energy. The frequency of bursts was 4 Hz. A fog-generator was used for seeding in order to produce seeding particles with average diameters of 1 ^m. The laser was placed downstream at approx. 1 m from the airfoil. A CCD camera with 1280 x 1024 pixels resolution was used, and the area covered was approx. 120 x 100 mm. The time-interval between the laser pulses was 20 ^s, and 32 x 32 pixel-size interrogation areas were used for velocity calculation. Cross-correlation and adaptive correlation were used, both with 25% overlap. Transparent solid sheet Figure 2: The PIV measurement system Several parameters have to be considered when estimating the uncertainties of PIV velocity measurements, [10]. Systematic errors occur due to uncertainties in determining the geometrical parameters and the fabrication tolerances of the camera's devices and lenses. Non-systematic errors are mainly due to the uncertainty when determining the average particle displacement within the interrogated region. These depend on the size of the interrogated region, the time separation between the laser pulses, the magnification of the recording, the out-of-plane velocity component, the turbulence, the length-scale of the flow etc. As the flow within the test-section was quasi-two-dimensional, the out-of-plane component of the vectors only caused negligible errors. Lehr and Boelcs, [11], showed that for these conditions the standard measurement uncertainty for the mean velocity field is less than 0.04 u^, and in the regions of strong velocity gradients it is smaller than 0.05 2.4 Numerical simulations 2D numerical simulations were performed for airfoil flow analysis. The simulations were conducted at the same geometry and angles of attack as used during the experiment. Steady-state and transient simulations were made, and the turbulence model selections were limited to k-e and SST (Shear Stress Transport) models. The Ansys 12.0 CFX code, [12], and the 2nd-order high-resolution scheme were used. An RMS residual below 10-6 or a stationary value of lift and drag coefficient were used as the convergence criterion. Computational meshes were created using ICEM-CFD. Precise analysis of the laminar transition boundary layer was possible only if the spatial resolution of the grid near the wall satisfied the condition y+<1, whilst the value y+>30 was required in the case of the k-e model. Since a scalable wall-function was used, the same computation mesh satisfying the condition y+<1 was used when applying the k-e model. A block-structured mesh type C with 119,000 nodes was used for 2D simulations. The simulation area was extended to two chord-lengths in front of the airfoil and to five chord lengths behind it. The tunnel's actual height (two chord lengths) was used. Fig. 3 shows an excerpt of the computational mesh close to the airfoil. No slip-boundary condition was used on the walls of the wind-tunnel. A constant velocity with a measured turbulence level of 1.5% and a static pressure boundary condition were applied at both the inlet and outlet, respectively. The GCI (Grid Convergence Index method) was used, [13], in order to calculate numerical uncertainty. The numerical uncertainty of the separation point location was 0.5%. Figure 3: Detail of the computational mesh close to the airfoil 3 RESULTS As mentioned, the presented investigations were focused on the study of turbulent layer separation from the suction side of the airfoil. This separation may happen at moderate AOA after the reattached turbulent boundary layer separates again from the airfoil surface. The turbulent boundary layer separation is greatly influenced by von Karman instability, which causes the separation point to oscillate. Simpson et al., [3], proposed a set of quantitative definitions on the detachment state near the wall using definitions based on the fraction of time the flow moves downstream. Four characteristic points were defined: incipient detachment (ID) occurring with an instantaneous backflow 1% of the time, transitory detachment (TD) occurring with an instantaneous backflow 50% of the time, and detachment (D) occurring where the time-averaged wall shear stress was equal to zero. Available data indicated that the TD and D points were at the same location, [3]. Finally, the characteristic point Dcp determines a critical position downstream of which the flow is detached at any time. In the following, the experimental results based on PIV images and pressure measurements are presented, followed by the numerical results and a comparison. 3.1 PIV images analysis The PIV system operated at frequencies of up to 4 Hz, and disallowed direct time-dependent analysis of the flow-field around the airfoil. Therefore, at least 60 consecutive images were acquired and their corresponding velocity fields established. The averaged velocity fields around the airfoil measured at four different AOA are presented in Fig. 4. No turbulent boundary layer separation occurred at 5° AOA. At 9.6° and 15.1° AOA, backflow vortex in the backward half of the airfoil indicated turbulent layer separation and a partially-stalled airfoil, whilst the burst of laminar separation bubble near the leading edge appeared at 20° AOA, and most of the suction surface was stalled, which agreed with the findings of Somers, [8]. Velocity [mfc]: 0 5 10 iS 2D 25 30 354Ö45 50 a) AOA = 5.1° b) AOA = 9.6° c) AOA = 15.1° d) AOA = 20.0° Figure 4: Velocity fields around the airfoil acquired by PIV, Re = 2 x 10 The fluctuating nature of the flow-separation was evident from successive PIV images. Two characteristic images that proved the oscillation of the flow-separation point at 15.1° AOA are presented in Fig. 5. Whilst in the right image flow separation occurred at x/c = 0.4, it shifted towards the trailing-edge of the airfoil by x/c = 0.6 in the left image. a) Separation point at x/c = 0.6 b) Separation point at x/c = 0.4 Figure 5: Two characteristic PIV images proving the oscillation of the flow separation point AOA = 15.1°, Re = 2 x 105 By careful analysis of a series of PIV images taken at the same AOA over a longer period of time, the separation zone could be further studied in order to predict not only both extreme points ID and Dcp, but also to estimate the percentage of back-flow at any point within the separation zone. A program was developed in LabVIEW to use the PIV results and extract all the values for each vector (velocity components, position and status). It was possible to analyse the data and calculate the medium value and standard deviation for the velocity and its direction (angle). It was also possible to count those cases in which the angle differed from the chord angle by more than a certain value. If the value was set at 90°, the result represented backflow occurrence. This value could be divided by the number of images, resulting in a backflow occurrence ratio. The results can be presented as intensity graphs, as shown in Fig. 6. a) AOA=15.1° b) AOA=20.0° Figure 6: Intensity graph of the backflow occurrence ratio; Re = 2 x 105 As can be seen from Fig. 6a (AOA = 15.1°), the alternations between the attached and separated flow (ID point) occurred at x/c = 0.4 (area with 0 to 10% of backflow). At x/c = 0.47 the flow was reversed 50% of the time; therefore, this point corresponded to the transitory detachment (TD) point, and at x/c = 0.62 the flow was detached at any time and corresponded to the Dcp point. The separation zone stretched between x/c = 0.4 and x/c = 0.62 at 15.1° AOA, but at 20° AOA (Fig. 6b) it comprised only a short region near the leading edge, actually representing the laminar separation bubble (LSB) which had moved forward and burst due to severe adverse pressure gradient. Similarly, the series of images at both the other two AOAs were analysed, and the results are presented in Table 1. The separation point oscillation-zone moved towards the trailing-edge when AOA decreased. In the case of 20° AOA, permanent flow-separation occurred at x/c = 0.11 with the bursting of the LSB, which was spread between 0.02 < x/c < 0.11 and the airfoil suction side was totally stalled. Partial-stall was observed at 15.1° AOA and 9.6° AOA. The turbulent boundary layer was first reattached and then the second separation took place with the oscillation zones between 0.40 < x/c < 0.62 for AOA = 15.1° and 0.48 < x/c < 0.66 for AOA = 9.6°, respectively. No separation of turbulent boundary layer was observed in the case of 5.1° AOA. Table 1: Comparison between the experimentally obtained positions (x/c) of the characteristic points (ID, TD, Dcp) of the flow separation zone and the predicted results Angle of attack (°) 5.1 9.6 15.1 20 ID TD Dcp ID TD Dcp LSB burst PIV images analysis NS* 0.48 0.52 0.66 0.40 0.47 0.62 0.02-0.11 Pressure signal analysis NS 0.5 0.6 0.3 0.6 0.0-0.2 Simulation, k-e turb. model NS NS 0.63 0.37** Simulation, SST turb. model 0.92 0.53 0.42 0.01** *NS - no separation occurred **LSB burst was not predicted by the k-£ and SST turbulence models; therefore, the data corresponds predicted flow separation position 3.2 Pressure tap signal analysis Eight pressure taps on each side of the airfoil (Fig. 1) were used for instantaneous pressure signal measurements. When averaged, these signals could be used to obtain pressure profile around the airfoil. Pressure is usually presented in the form of a pressure coefficient: r ■ kpuL (3.1) where p^ is the free-stream static pressure and u^ the free-stream velocity far upstream from the airfoil. In order to validate the obtained pressure profiles, a comparison was made with the results obtained by Sommers et al., [8]. Fig. 7 shows a comparison between pressure coefficient profiles at AOA = 9.6°. The agreement was good and ensured that the presented measurements are correct. & 0,2 0,4 0,6 0,8 x/c Figure 7: Pressure coefficient profile for AOA = 9.6° 3 2 1 0 1 2 0 1 Besides obtaining the tap position's averaged pressure, the acquired pressure signals may also be used to calculate the signal standard deviation, which represents the measure of pressure fluctuation intensity at a particular tap position. Along with the averaged pressure coefficient, its standard deviation for the airfoil suction side is also presented in Fig. 7. As can be seen, it varied alongside the airfoil with a maximum value at approximately x/c = 0.5. According to Sicot et al. [14], the oscillation zone of the flow-separation point can be obtained using standard deviation of the signals from the pressure taps along the suction-side of the airfoil. The curves connecting these values for AOA = 5.1°, where (according to PIV images) no stall was presented and AOA = 9.6°, where the airfoil was partially stalled, were compared and are shown in Fig. 8. The difference was quite evident. The curve for AOA = 5.1° has no local extremes, whilst these are clearly evident in the curve for AOA = 9.6° which had two characteristic points. The first one corresponds to the incipient detachment point (ID). This ID point represents the position of the separation point where the alterations occurred between the attached and separated flows. Therefore, the pressure fluctuations at the ID point were the greatest, thus coinciding with the local maximum value for normalized standard deviation. Upstream of ID, the flow is attached at any moment. The characteristic point Dcp downstream from which the flow is detached at any time corresponds, according to Sicot et al. [14], to the local minimum of standard deviation. The interval between the ID and Dcp points locates the oscillation zone of the separation point. The error made using this method comes from the discontinuous distribution of the pressure taps. The resolution was quite low in the presented case, since only 8 pressure taps were applied on the airfoil suction's surface. Similarly, the series of pressure signals at all other AOA were analysed, and the results are presented in Table 1 for comparison with the results obtained from PIV images, and the numerical predictions. The presented pressure signal based method adopted from Sicot et al. [14] predicts only both extreme points (ID and Dcp) within the separation point-oscillation zone, and gives no quantitative insight into it unlike the PIV method. However, it is very robust and may be applied for the first approximation of the position and length of the separation point oscillation zone. In the presented case (Table 1), the predicted separation point oscillating zone agreed fairly well with the results of the PIV based method. The agreement was very good especially at 9.6° AOA. The results at 15.1° AOA were less accurate, whilst the total airfoil suction surface stall at 20° AOA was again correctly predicted. Figure 8: Variations in the normalised standard deviations of pressure signals from successive pressure taps on the suction side of the airfoil for AOA = 5.1° and AOA = 9.6°, Re = 2 x 10 3.3 Numerical simulation The last part of the presented study was devoted to numerical simulations. The idea was to obtain better insight into the separation area, in order to understand any possible mechanisms responsible for turbulent layer separation, and oscillation of the separation point. However, the presented work was limited to the application of 2D U-RANS simulation, due to limited computational resources otherwise necessary to run 3D LES or DNS. According to, [15], the 3D U-RANS simulation results appear to be close to a spanwise repetition of the 2D field. The reason for this is due to the role of the 3D spanwise random perturbations that affect the vortex shedding within the real flow. These effects are not described by U-RANS computation, which is intrinsically deterministic, [16]. The 3D U-RANS simulations were therefore omitted. Two-dimensional steady state simulations were performed during the first step in order to optimise the computational mesh and to obtain the initial conditions for unsteady simulations. Two turbulence models, k-e and SST model, were applied. The k-e turbulence model is the most-widely used turbulence model and enables acceptable predictions when separation is absent (small angles of attack), [17]. However, this model does not consider the transport of shear stress under conditions where separation is caused by an adverse pressure gradient, thus resulting in an overestimation of eddy viscosity and, consequently, estimating both the separation point and separation area inaccurately, resulting in a stalled delay phenomenon. Menter et al., [18], developed the SST turbulence model by improving some weaknesses of the k-e model. Compared to the k-e model, SST provides more accurate predictions of the separation point's position, and the separation area caused by an adverse pressure gradient. In order to check for the accuracy of simulations, a comparison was made between the measured and computed lift and drag coefficients. The results are presented in Fig. 9. At instances of lower AOA of up to 9.6°, the lift and drag coefficients were well predicted, whilst both the k-e and SST models, respectively, overestimated the lift coefficient and underestimated the drag coefficient at higher AOA, where the airfoil was partially or fully stalled. As was discovered during the experiments performed by the PIV, as well as by instantaneous pressure measurements, the separation point oscillated significantly, especially when the airfoil was partially stalled. It was, therefore, suspected that the poor accuracies of the simulations at high AOA were caused by the steady-state numerical computation. Thus, transient 2D simulations using a SST turbulence model were then carried out; however, the results did not differ significantly from those obtained using steady-state simulations, except for when the AOA = 20°. An oscillating Von Karman vortex street, formed behind the airfoil after the simulation was started, slowly died away, and the final result was a totally-static flow field around the airfoil with no oscillations presented, and the separation point remained static. It is assumed that the reason for this is the overly large turbulent viscosity of RANS turbulent models, which are, therefore, highly dissipative and are unlikely to be triggered into unsteady mode unless the flow instabilities are large, [19]. The latter was the case at AOA = 20°. At this very high AOA, the reattachment of a turbulent boundary layer failed, and total stall took place with a fully turbulent wake behind the leading edge. The predictions were thus in good qualitative agreement with the experimental results (Fig. 4d); however, the damping of oscillations was still too high. The predicted Strouhal number, which should have been around St = 0.2 [14], was less than 0.1 and the mean predicted lift and drag coefficients remained overestimated. 1,2 lift (experiment) drag (experiment) lift (k-£) • •••!»••• drag (k-£) — — lift (SST) - •• - drag (SST) 0 5 10 15 20 25 Angle of attack [°] Figure 9: Comparison between experimentally and numerically obtained lift and drag coefficients; Re = 2 x 10 Since the 2D numerical simulations of partially stalled airfoils at 9.6° and 15.1° AOA could not show any oscillation of the flow-separation point location and predicted only the permanently-detached flow vortex (Fig. 11), it was impossible to determine all the characteristic separation points proposed by Simpson et al., [3], except for the Dcp point location. However, when considering what the steady-state result of an oscillating phenomenon is, this point can be regarded as point TD, where the flow is separated 50% of the time. Table 1 presents both the experimentally and numerically-predicted flow-separation points. There were substantial differences between the numerical results from different turbulence models. Both the experimental results and the k-e turbulence model's prediction showed no separation on the suction side of the airfoil, whilst the SST turbulence model predicted a small vortex on the trailing edge at 5.1° AOA. The difference increased with any increase in the angle of attack. At 9.6° AOA, no vortex was predicted by the k-e model, although some reduction in velocity was observed, whilst the SST model predicted the flow-separation at x/c = 0.53. At 15.1° AOA, a difference in vortex location was also observed. The flow-separation point predicted by the k-e model was located at x/c = 0.63, whilst the SST model predicted this location at x/c = 0.42. In the case of 20° AOA, a simulation using the k-e model predicted separation at x/c = 0.37, whilst the SST model predicted no reattachment of the turbulent layer, and the flow separated completely at x/c = 0.01 (Fig. 10). Fig. 10 shows a comparison between flow separations predicted using the SST model at 15.1° AOA and 20.0° AOA. At AOA = 15.1° (Fig. 10a), a small vortex was formed at the suction side of the airfoil approximately 0.04 c from the leading edge. The flow separated due to sudden directional change when hitting the airfoil and then reattached to the airfoil surface using the high turbulent energy of the main flow. It remained attached to the airfoil until another separation occurred at x/c = 0.42 due to the adverse pressure gradient. No leading edge vortex was formed at AOA = 20.0° (Fig. 10b). The main flow turbulent energy was insufficient to reattach the flow again, and a total stalling of the airfoil occurred. It has to be pointed out that the vortex predicted at AOA = 15.1° should not be The unsteady static-stall aerodynamic characteristics of an S809 airfoil at low Reynolds numbers confused with the laminar separation bubble, since the simulation was fully turbulent and could not predict it. a) AOA = 15.1° b) AOA = 20.0° Figure 10: Flow separation predicted using the SSTmodel; Re = 2 x 105 Although both models, k-£ and SST, failed to predict the time-dependent behaviour of the turbulent layer separation phenomenon that was established by the experiment at 9.6° and 15.1° AOA, their results could be used as a time-averaged state and (as will be seen later) do agree fairly well with the averaged experimental results, especially for the results of the SST model. Fig. 11 shows a comparison between the predicted and experimentally-obtained time-averaged flow-fields. The latter were obtained by averaging all the instantaneous PIV images corresponding to a particular AOA. A much better agreement with the experimental data was achieved by the SST model, which more accurately predicted the separation vortex at the suction side of the airfoil, in comparison with the k-£ model. In contrast, both turbulence models predicted similar velocity fields within the leading edge area and along the whole pressure side of the airfoil. No comparison between velocities within the boundary layer was possible, due to the intense laser light reflection near the airfoil's surface, which prevented any velocity measurements within this area. i«» o* 0 4 02 ■«M|H«| t i »H»»»»«»«» ? 1J 01 0« 04 02 * a) PIV b) k-£ model c) SST model Figure 11: Comparison between the velocity fields of the numerical and averaged experimental results AOA = 15.1°; Re = 2 x 105 The experimentally-obtained average vorticity field (Fig. 12a) showed the formation of two shear layers: one from the leading edge and one from the trailing edge. The formation length was greater from the leading edge (/2) than from the trailing edge (/1). A sudden growth of the average shear layer, combined with a roll-up of the shear layer, was observed downstream of the airfoil. Similar results were obtained by Sicot et al., [14], and Yang et al., [20], when they both investigated the flow-field around an isolated airfoil at moderate Re numbers. The numerical simulation using the k-e model predicted no flow separation at 9.6° AOA (see Table 1); thus, no shear-layer was observed from the leading edge (Fig. 12b). The SST model predicted both shear layers quite accurately, although their structures were more compact, and both were longer. a) PIV b) k-e model c) SST model Figure 12: Comparison between the vorticity fie/ds of the numerica/ (2D) and averaged experimenta/ resu/ts AOA = 9.6°; Re = 2 x 105 4 CONCLUSIONS The unsteady static stall aerodynamic characteristics of an S809 airfoil at low Re numbers were analysed, using both experimental and numerical approaches. Several angles of attack were studied. The PIV system was used to acquire a series of captures in order to obtain instantaneous velocity fields around the suction-side of the airfoil. At low AOA, no flow separation was observed, whilst the already moderate AOA caused the reattached turbulent boundary layer to separate again near the trailing edge of the airfoil, which was therefore partially stalled. With any increasing of AOA, the separation point moved increasingly towards the leading edge until reattachment failed and the airfoil became totally stalled. Instantaneous PIV captures revealed the oscillating nature of the separation point. A simple numerical procedure was applied to obtain a characteristic interval for the separation point's positional oscillation spreading between the permanent separation point and the starting separation point. The experimental part also included pressure profile measurements along the airfoil, and time-dependent measurements of the pressure-oscillation within the specified points around it. The data were analysed and specific points were located for defining the permanent separation points and starting separation points. The results from both experimental methods agree fairly well. An attempt was also made to simulate the separation point's positional oscillations with CFD simulations. A transient 2D approach using k-e and an SST turbulence model was used; however, both models failed to predict the separation point's positional oscillations. The possible reason for this is due to the role of the 3D spanwise random perturbations that affect the vortex shedding in the real flow, [21]. Although the simulations were unsuccessful in predicting the oscillatory nature of the separation point, their results showed reliable agreement with the averaged velocity and vorticity fields around the suction side of the airfoil obtained by PIV. The SST model showed better agreement with the experimental data. Although the 3D effects were omitted, the predicted shear layer formations at the leading and trailing-edges agreed well with the averaged experimental results. References [1] Lissaman, P.B.S., Low-Reynolds-number airfoils, Annual Review of Fluid Mechanics, Vol. 15, 1983, pp. 223-239. [2] Mueller, T.J. (ed.), Fixed and flapping wing aerodynamics for micro air vehicle applications, Progress in Astronautics and Aeronautics, ISBN 1-56347-517-0, AIAA, 2001. [3] Simpson, R.L., Chew, Y.T., Shivaprasad, B.G., The structure of a separating turbulent boundary layer: part I, mean flow and Reynolds stresses; part II, higher order turbulence results. The Journal of Fluid Mechanics 113, 1981, pp. 23-73. [4] Nishimura, H., Taniike, Y., Aerodynamics characteristics of fluctuating forces on a circular cylinder. Journal of Wind Engineering and Industrial Aerodynamics, 89, 2001, pp. 713723. [5] Kim, B.S.; A Study on the Optimum Blade Design and the Aerodynamic Performance Analysis for the Horizontal Axis Wind Turbines, Ph.D. Thesis, Korea Maritime University, 2005. [6] Wolfe, W.P., Ochs., S.S., CFD calculation of S809 aerodynamic characteristics. AIAA- 970973, 1997. [7] Devinant, P., Laverne, T., Hureau, J., Experimental study of wind-turbine airfoil aerodynamics in high turbulence, Journal of Wind Engineering and Industrial Aerodynamics, 90, 2002, pp. 689-707 [8] Sommers, D.M., Design and Experimental Results for the S809 Airfoil, NREL/SR-440-6918, 1997. [9] Chapin, W. G., Dynamic-pressure measurement using an electronically scanned pressure module, NACA TM-84650, 1983. [10] Raffel, M., Willert, C., Kompenhans, J., Particle Image Velocimetry: A Practical Guide, ISBN 3-540-63683-8, Springer-Verlag., 1998. [11] Lehr, A., Bölcs, A., Application of a particle image velocimetry system to the investigation of unsteady transonic flows in turbomachinery, 9th International Symposium on Unsteady Aerodynamics, Aeroacoustics and Aeroelasticity of Turbomachines, Lyon, 2000, pp. 4-8. [12] Ansys. CFX-SolverTheory Guide, Release 12.1, (2009). [13] Celik, I. B., Ghia, U., Roache, P. J., Procedure for estimation and reporting of uncertainty due to discretization in CFD applications. Journal of Fluids Engineering, Vol. 130, No. 7., 2008. [14] Sicot, C., Aubrun, S., Loyer, S., Devinan, P., Unsteady characteristics of the static stall of an airfoil subjected to freestream turbulence level up to 16%. Experiments in Fluids, 41, 2006, pp. 641-648. [15] Jacob, M.C., Boudet, J., Casalino, D., A rod-airfoil experiment as benchmark for broadband noise modeling, Theoret. Comput. Fluid Dynamics, 19, 2005, pp. 171-196. [16] Casalino, D., Jacob, M.C.: Prediction of aerodynamic sound from circular rods via spanwise statistical modeling. Journal of Sound and Vibrations, 262, 2003, pp. 815-844. [17] Jošt, D., Lipej, A., Numerical prediction of the vortex rope in the draft tube. Proceedings of the 3rd IAHR International Meeting of the Workgroup on Cavitation and Dynamic Problems in Hydraulic Machinery and Systems, Brno: University of Technology, Brno, 2009, pp. 14-16. [18] Menter, F.R., Two-equation eddy-viscosity turbulence models for engineering applications. AIAA Journal, 32(8), 1994, pp. 1598-1605. [19] Davidson, L., Evaluation of the SST-SAS model: channel flow, asymmetric diffuser and axi-symmetric hill, European Conference on Computational Fluid Dynamics ECCPMAS CFD 2006, Delft, the Netherlands, 2006. [20] Yang, Z., Haan, F. L., Hui, H., Ma, H., An Experimental Investigation on the Flow Separation on a Low-Reynolds-Number Airfoil, 45th AIAA Aerospace Sciences Meeting and Exhibition, Reno Nevada, USA, AIAA-2007-0275, 2007. [21] Lammers, P., Jovanovic, J., Durst, F., Numerical experiments on wall turbulence at low Reynolds numbers, Thermal Science; Vol. 10, No 2, 2006, pp. 33-62 JET Volume 6 (2013), 51 - 60 Issue 1, February 2013 http://www.fe.um.si/en/jet/e-jet.html FRACTURE EVALUATION OF ENERGY COMPONENTS WITH LOCAL BRITTLE ZONES VREDNOTENJE LOMA V ENERGETSKIH KOMPONENTAH Z LOKALNIMI KRHKIMI PODROČJI Zdravko Praunseism Keywords: fracture, energy components, sectioning technique brittle zones Abstract The fracture toughness in the testing of multi-pass welds and heat-affected zones is remarkably sensitive to the microstructures in the vicinity of the crack tip of test specimen. Therefore, the introduction of the pre-crack to the weld and heat-affected zone specimen should be done most carefully. However, since there is an uncertainty of the crack tip position in fatigue pre-cracking, it becomes common to section near the fatigue pre-crack front after testing is complete, and to examine the cross section in order to identify the position in weld and heat-affected zone microstructures, the so-called local brittle zones. Concerning this sectioning technique, the precise experimental procedure is specified in this article. Povzetek Lomna žilavost večvarkovnih zvarov in toplotno vplivanih področjih je odvisna od vrste mikrostrukture na konici razpoke v preizkušancu. Zato je nujna natančna postavitev utrujenostne razpoke v zvar in toplotno vplivano področje preizkušanca. Postavitev utrujenostne razpoke ni natančno določena v nobeni proceduri, je pa postala praksa, da se m Corresponding author: Zdravko Praunseis, PhD, Faculty of Energy Technology, University of Maribor, Tel.: +386 31 743 753, Fax: +386 7 620 2222, Mailing address: Hočevarjev trg 1, SI-8270 Krško, Slovenia, E-mail address: zdravko.praunseis@uni-mb.si jet Journal of Energy Technology razrez preizkušancev izvrši v neposredni bližini utrujenostne razpoke in sicer tako, da prerez omogoča natančno določitev lokacije utrujenostne razpoke v zvaru in toplotno vplivanem področju ter pripadajoče mikrostrukture, tako imenovanih lokalno krhkih področij. V članku je natančno opisana eksperimentalna procedura za razrez in ocenitev preizkušancev. 1 INTRODUCTION The crack tip opening displacement (CTOD) test has been become a common method of measuring the fracture toughness of steel weldments. Nevertheless, the commonly used fracture mechanics testing standards, including the CTOD testing standards, such as BS 5762, [1], assume the use of metals with high degrees of homogeneity, although not explicitly emphasised. As already mentioned, welded joints have typical macroscopic heterogeneity and residual stresses as a result of welding. In order to clarify the applicability of the common testing methods, a basis of knowledge taking the above heterogeneity into account must been established outside the standards. Recently, some activities have been conducted for establishing the CTOD testing procedure of steel welds and some recommended practices/guidelines for CTOD tests of weldments have been published. It is widely understood that the fracture toughness is considerably affected by the shape of the crack front of the fracture toughness specimen [2-4]. Therefore, in the common fracture toughness specimen, attention is carefully paid to realise a straight crack front perpendicular to the plate surface. However, in the welded joint it is sometimes very difficult to obtain a straight crack front of the fatigue pre-crack due to the existence of weld residual stresses. In order to avoid the confusion due to the irregularity of the crack front and to realize reproducibility, in the current standards it is required that as straight a crack front as possible be achieved, [5]. 2 NOTCHING AND SECTIONING PROCEDURE TECHNIQUE In order to achieve a uniform fatigue crack shape that meets the standard requirements, some treatments, i.e. residual stress-relieving treatment, have to be applied for notched specimens of welded joints. A different method for relieving residual stresses is to impose a local plastic strain on the region suffering from residual stresses; the following techniques [2-5] are currently in use: • local compression, • reverse bending, • the use of a high R-ratio in the cycle, and step wise high R-ratio method, • both sides hole method. Table 1 gives a summary of the relative merits of the three methods of them. In the Recommended Procedure proposed by The Welding Institute, [6, 7], the mechanical relief of residual stresses by local compression, in which a plastic strain of 1% of the specimen thickness, is recommended. According to the institute, "the use of reverse bending prior to fatigue pre-cracking as a means of redistributing welding stresses is not recommended." Moreover, "the effect of a high R- ratio on the fracture toughness is not well understood and so until more work has been completed on this technique its use is not generally recommended." However, for very thick section welds, "the use of high R-ratios during fatigue pre-cracking has been found to be successful in obtaining acceptable crack front profiles." Table 1: Characteristics of materials used in stress relieving Advantages Disadvantages LOCAL COMPRESSION. -method well published -method been in use since 1975 -uses normal fatigue pre-cracking procedures -requires extra operation -requires high capacity compression rig and tools -toughness may be conservative for some materials -specimen must be flat REVERSE BENDING. -special equipment not needed -conservative toughness measurements expected -uses conventional fatigue pre-cracking procedures -requires extra operation -toughness may be significantly lower -little information published USE OF HIGH R-RATIO. -no extra operation needed -no extra equipment needed -required loads and R-ratios in conflict with limits of current standards -little information published -non-conservative assessments of toughness are expected In the common fracture toughness test, the use of a notch sharpened by a pre-crack produced by fatigue loading of the test piece is generally required in order to simulate sharp macroscopic defects in the structure and to provide a conservative assessment of toughness. In order to avoid confusion and to realize reproducibility, the condition of the fatigue pre-cracking loading must also be kept within limits. After the CTOD test is conducted, both halves (or the half containing the weld metal) of the broken specimen are sectioned and metallurgically examined. The cut into the fracture face is taken just behind, but within 2.5 mm, of the fatigue-crack front. The cross section may contain a portion of the fracture surface near one or both surfaces due to the fatigue-crack front curvature. Each such portion is not wider than 10% of the specimen thickness. For CTOD specimens that are notched to sample the coarse grain (CG) regions, quantification is as shown in Fig. 1, where the linear fraction of the CGHAZ region sampled by the fatigue crack is calculated. A similar procedure is used for the inter-critical coarse grain (IC) and subcritical coarse grain (SC) HAZ areas. Fatigue-crack sampling calculations are made by examining enlarged photographs (3-6 times magnification) of the CTOD cross sections. Figure 1: Sectioning both halves of a HAZ CTOD specimen to calculate CGHAZ percentage By using both halves of the broken CTOD specimen and the enlarged photographs, fatigue-crack sampling calculations can be made with reasonable accuracy without microscopic examination. Each HAZ specimen should be sectioned to determine the regions of microstructure sampled by the fatigue crack. Figure 2: Example of sections taken from a HAZ, through-thickness-notched CTOD specimen to identify microstructures sampled by fatigue crack and at the fracture initiation point. In the case of through-thickness notched specimens, this is best achieved by sectioning at a small distance behind the fatigue crack tip, so as to include as much of the fatigue crack front as possible (Fig. 2). With surface notched specimens, a similar approach could be used. However, where the region being sampled is small and/or the fatigue crack front is bowed, misleading results may be obtained. For this situation, a better approach is to section as shown in Fig. 3(b), and if necessary, take a series of sections. Figure 3: Example of sectioning techniques for: a) through-thickness notched, and b) surface notched specimens It is recommended that similar sectioning procedures are applied to all tests (HAZs and weld metals) carried out to measure the fracture toughness associated with known cracks, [5]. It is generally agreed that there is an additional requirement to establish the microstructure at the fracture initiation point; detailed fractography is necessary to determine the microstructure at that point and hence locate the position from which the section has to be taken. A practical example, [5], of the sectioning procedure is shown in Figure 4 and Figure 5. 3 DISCUSSION OF RESULTS The cause for differences, [2-4], of HAZ fracture toughness (Fig. 4) was the effect of the root weld metal strength mismatch (overmatched root of homogeneous weld and under-matched soft root layer of heterogeneous weld) on the deviating nature of stable crack growth and on the initiation of the specimens' final brittle fractures. Specifically, in CTOD specimens with homogeneous and heterogeneous welds, the whole fatigue crack tip front was located in CG HAZ (Fig. 4 - Cross-section G-G, Detail A and B), as also indicated with the high values of micro-hardnesses (« 340 HV1). All specimens with homogeneous welds reached CTOD values Su after initial stable crack growth (Fig.4 - Fracture). Owing to the root overmatched weld metal and its shielding effect, the stable crack growth path deviated towards the base material (Fig.4 - Cross-section P-P), so that the crack tip front first crossed the brittle part of CG HAZ (appearance of a small pop-in); then, the tougher fine grain (FG) HAZ followed by the final specimen fracture in IC HAZ (Fig.4 - Detail A). The influence of root overmatched weld metal (shielding effect) of the under-matched homogeneous weld on a deviation of fracture path was (in this example) as expressed as for the specimens with so-called composite fatigue crack front (a/W=0.5) in HAZ, where the crack growth path also turned toward the softer base material. Soon after the start of the loading, local brittle zones (LBZ) in the form of pop-ins appeared in the HAZ-notched CTOD specimens taken from heterogeneous weld, because the crack tip reached the region of low toughness - IC CG HAZ (Fig. 4). The presence of expressive local strength mismatches at narrow CG HAZ and the global mismatching effect of the softer under-matched root layer has caused deviations of stable crack growth paths to the regions of low toughness, fusion lines and under-matched weld metals (Fig. 4 - Cross-section G-G and Detail A). The specimen brittle fracture appeared in the bainitic microstructures (Fig.4 - Detail C) of weld filler metal. The brittle fracture initiation point was indicated on the specimen fracture surface (Fig. 4 - Fracture) by EDX analysis as a carbide Fe3C, [5] The CTOD - S5 values, [5], for through-thickness notched specimens (Fig. 5) in the heterogeneous weld indicate a reduction of the fracture toughness with an increase of soft root volume, i.e. for soft root passes thickness h > 8mm. In the aforementioned CTOD specimens, the brittle fracture occurred in the soft root layer centre (Fig. 5 - Fracture surface) without previous origin, at the position where the mismatch factor was the lowest, as well as hardness, 161 HV1 (Fig. 5 - Detail A), being a consequence of very low toughness of soft root layer. The main cause for the low soft root layer toughness was a change of the microstructure of the all-weld metal obtained by the welding wire, which was exposed to different alloying mechanisms during welding in the root region. Specifically, the all-weld metal microstructure of soft root weld metal was bainitic, having high toughness, but because of the aforementioned reasons its microstructure was changed to inconvenient ferritic microstructures (Fig. 5 - Section L-L and Detail A), causing low toughness of the soft root layer. fatique crack stabile crack growth fracture Detail B Section G-G fatique crack crack growth path deviation brittle fracture brittle fracture brittle fracture initiation point Figure 4: Position of surface (a/W = 0.254) fatigue crack in CG HAZ and crack path deviation in specimen with heterogeneous weld joint Fracture Section K-K Figure 5: Fracture and microstructure at the vicinity of brittle fracture initiation point of through-thickness notched specimen in heterogeneous weld 4 CONCLUSIONS The crack tip opening displacement (CTOD) test has been become a common method of measuring the fracture toughness of steel weldments. The fracture toughness in the multi-pass welds and heat-affected zones testing is remarkably sensitive to the microstructures in the vicinity of the crack tip of the test specimen. In the common fracture toughness test, the use of a notch sharpened by a pre-crack produced by fatigue loading of the test piece is generally required in order to simulate sharp macroscopic defects in the structure and to provide a conservative assessment of toughness. After CTOD test is conducted, both halves (or the half-containing the weld metal) of the broken specimen are sectioned and metallurgically examined. Detailed fractography is necessary to determine the microstructure at that fracture initiation point and hence locate the position from which the section has to be taken. Acknowledgements The authors wish to acknowledge the financial support of the Slovenian Foundation of Science and Technology and the Japanese Promotion of Science References [1] BS 5762, Methods for crack opening displacement (COD) testing, The British Standards Institution, London 1979. [2] Praunseis, Z., Toyoda, M., Sundararajan, T.: Fracture behaviours of fracture toughness testing specimens with metallurgical heterogeneity along crack front. Steel res., Sep. 2000, 71, no 9, [3] Praunseis, Z., Sundararajan, T., Toyoda, M., Ohata, M.: The influence of soft root on fracture behaviors of high-strength, low-alloyed (HSLA) steel weldments. Mater. manuf. process., 2001, vol. 16, 2 [4] Toyoda, M.: Fracture toughness evaluation of steel welds, Book of Mechanics, Osaka University, 1989. [5] Praunseis, Z.: The influence of Strength Under-matched Weld Metal containing Heterogeneous Regions on Fracture Properties of HSLA Steel Weld Joint (Dissertation in English). Faculty of Mechanical Engineering, University of Maribor, Slovenia, 1998 [6] ASTM E 1152-87, Standard test method for determining J-R curves, Annual Book of ASTM Standards, Vol. 03.01, American Society for Testing and Materials, Philadelphia, 1990. [7] ASTM E 1290-91, Standard test method for crack-tip opening displacement (CTOD) fracture toughness measurement, American Society for Testing and Materials, Philadelphia, 1991. Zdravko Praunseis JET Vol. 6 (2013) Issue 1 JET Volume 6 (2013), 61 - 68 Issue 1, February 2013 http://www.fe.um.si/en/jet/e-jet.html GENERIC MODEL OF WIND TURBINE BLADES GENERIČNI MODEL LOPATIC VETRNE TURBINE Gorazd Hrenw, Andrej Predin, Ivan Žagar Keywords: CAD, wind turbine blades geometry, airfoils Abstract A wind turbine is a device that extracts energy from the wind, which has been shown to be one of the most viable sources of renewable energy. Current manufacturing technology enables the low-cost production of wind energy turbines, which is competitive with conventional sources of energy, especially in steady wind areas. The rotor blade is a key element in a wind turbine system. Most commercially available blades incorporate airfoil-shaped cross sections, which have been found to be very efficient at lower wind speeds. Computational fluid dynamics solve and analyse problems involving fluid flows. Many airfoil shapes are combined and analysed in order to achieve efficient blade design. In order to prepare various geometrical solutions for computational fluid dynamics analysis, time-consuming and error-prone work has to be done. In this paper, the automation of the process of preparing the various geometrical models for computational fluid dynamics analysis of wind turbine blade is described. Povzetek Vetrna turbina je naprava, ki pridobiva energijo iz vetra in izkazano eden izmed najbolj uspešnih obnovljivih virov energije. Aktualna tehnologija proizvodnje vetrnih turbin omogoča nizko ceno in stroške vetrne energije, ki je konkurenčna znanim konvencionalnim virom energije, posebej na območjih s stalnim vetrom. Lopatica rotorja je ključni element v sistemu vetrne turbine. m Gorazd Hren, PhD., Tel.: +386 7 220 7218, Fax: +386 7 220 7210, Hočevarjev trg 1, 8270 Krško: E-mail address: gorazd.hren@um.si Journal of Energy Technology Večina komercialnih vetrnih turbin ima lopatice v obliki profila letalskega krila, za katere je ugotovljena večja učinkovitost pri nižjih hitrostih vetra. Numerična dinamika tekočin omogoča reševanje in analizo tokov tekočin z uporabo računalniških numeričnih metod. Da bi dosegli učinkovito obliko lopatic je mogoče geometrijske oblike lopatic kombinirati in analizirati. Pripravljanje geometrije lopatic za numerične analize je zamudno in natančno delo. Članek opisuje avtomatizacijo procesa priprave različnih variant vetrnih turbinskih lopatic primernih za numerične analize dinamike tekočin. 1 INTRODUCTION One of the major challenges in this century is the efficient use of energy resources, as well as the growing production of energy from renewable sources. There are several alternative forms of energy that have been explored and developed, including geothermal, solar, wind and hydroelectric power. The affordability and performance of renewable energy technologies is the key to ensuring its availability to the market. A wind turbine uses the aerodynamic force of the lift to rotate a shaft, which in turn aids in the conversion of mechanical power to electricity by means of a generator. For large networks, modern wind turbines are connected to the grid in a wind farm in order to reduce the total electrical load. One significant downside of wind as a resource is its variable speed and fluctuating electrical output. The most common design of a wind turbine is the horizontal axis wind turbine, in which the axis of rotation is parallel to the ground, as shown in Figure 1. Figure 1: Wind farm in Lower Saxony, Germany, [1]. The rotor consists of the hub and the blades for the wind turbine. These components are often considered the most important for cost efficiency. Today, most designs have three blades, and for better performance some designs implement pitch control for the angle of rotated blade. The turbine blades are manufactured from composite materials, including fiberglass-reinforced plastic, epoxy and wood laminates. To achieve the maximum power from a given wind situation, the geometry of the blade is designed to provide the maximum power coefficient from the wind side, usually taking into consideration the average wind speed. Numerical solutions based on analytical equations are adequate tools for the design of a wind turbine blade in order to get maximum power from a given wind, [2]. Wind turbine blades' conceptual design is inherently a multi-disciplinary design process that involves numerous disciplines and forms of expertise. In this paper, how high-end CAD software can be used to automate the creation and preparation of different varieties of wind rotor blades' geometry for computational fluid dynamics analysis is investigated. The generic model that is developed in this regard is able to automate the process of the creation and modification of the wind blade geometry, based on a series of parameters that define the geometrical characteristics of wind blades using V5 CATIA software. 2 WIND TURBINE AERODYNAMICS During the operation of a horizontal axis wind turbine, the properties of the wind flowing around it are constantly changing in relation with the distance upstream or downstream of the flow field. Far upstream from the turbine, a circular boundary region starts forming, giving a cylindrical shape to the flow field. Such a boundary is marked by the rotor blades sweeping around the turbine axis, and it defines the so-called actuator disc concept, [2]. The NACA airfoil series were generated using analytical equations that describe the curvature of the geometric centre line of the airfoil section, as well as the section's thickness distribution along the length of the airfoil. Some researchers have used blended wings with two different types of airfoils in order to achieve the desired design. 3 NUMERICAL SIMULATIONS Computational Fluid Dynamics (CFD) has been developed from mathematical methods and become an essential tool for almost any problem in fluid dynamics. It is commonly recognised as a numerical solution (by computational methods) of the governing equations that describe fluid flow. As a developing science, CFD has received extensive attention throughout the international community. All CFD analysis starts by defining the geometry to be used. Numerical solutions based on analytical equations that allow the design of a wind turbine blade are available tools that can be used to extract the maximum power from a given wind situation, [2]. Usually, the geometry of the blade is obtained to provide the maximum power coefficient from a wind site, using a proper average wind speed, which is the result of adequate wind data measurements. The analysis has to consider the fact that wind power is proportional to the cube of the wind velocity. A complete design of a wind generator is a complex problem that includes the design of an optimum blade and many other turbine components, such as the tower design. While the rotational velocity of the blade can be accurately predicted with analytical models, the same rotational velocity cannot be properly predicted for wind speeds below and above the rated wind speed, [5]. The CFD method is a well-known technique applied to the solution of more complex fluid problems. A performance evaluation of a wind turbine rotor is one of these difficulties, [6]. 4 CAD AUTOMATION This article focuses on the automation of the design process of complex geometry such as wind turbine blades. The wind turbine models are creating mostly in CAD systems using Graphical User Interface to create geometry from cross-sections. The use of customised application interfaces can accelerate preparations of a model, reduce the numbers of iterations in input and make the generation of models easier to achieve. Since we consider the models to be investigated by numerical simulations, the same approximations are performed during geometry model generation to accelerate the process of mesh generation. The design of an airfoil cross-section is governed by parametric equations based on form generation standards (e.g. NACA series). The cross-sections are created by fitting spline curves through points generated from equations. The developed programs, called macros or scripts, are used to input points, create splines and finally generate loft surfaces and body. To input every single point by coordinates for every cross-section is very tedious and error-prone work. Modifications are even more difficult. Macro programming is one of several automation techniques that are usually included within CAD systems. Available techniques could be classified in three types: • Kernel APIs refers to a core functionality embedded within CAD systems. ACIS and Parasolid are two commercially available kernels that are used in CAD systems. The kernel functionality is made available through libraries of APIs, the functions of which can be called in programs written in a high-level programming language. They require knowledge of computer programming with high analytical content. • Macros are scripted instructions that are executed in interpretive way by the process running simultaneously within CAD system. Many CAD systems have supplemented their original macro programming languages with VBA; consequently, the main advantage of using VBA is utilisation integration and simplification with other applications based on the Windows framework. Macro programming provides the possibility of executing the same operations that could be performed by the GUI. Macros execute more slowly than compiled code but still faster and more error free than manual work with such a cloud of points. • Knowledge-Based automation tools have been developed to support the development of Knowledge-Based Engineering applications inside the CAD environment. The design process can be captured, leading to the reduction of redesign iterations for new designs. While macros emulate interactions with the CAD system, these tools emulate engineering processes that are related to creating new geometry or making changes to the existing one. Training and in-depth computer knowledge is required for these tools. From the three described techniques that enable automation in the design process, the use of macros is best suited for our purpose. The other two tools are more complex and require programming skills with non-procedural programming languages and deep understanding the basic theory of solid modelling. CAD automation is useful when the modelling tasks are practically impossible to do manually because the very large number of points, and great effort needed to create such geometry is time consuming, error prone and very precise and delicate work. Automation is welcome when there are more users who need to repeatedly reuse variants of the same components and when the automation process reduces the time needed for managing large numbers of models or models with such complex geometry with a large number of calculated points. To obtain the geometry of the wind turbine blade, Excel is used as the starting application. In spreadsheets, it is easy to manage and prepare all necessary data for easier use of VBA macros in a suitable manner for transferring and managing data in CAD system, [4]. Figure 2: Blade profiles from NACA 4415. In CATIA, the VBA Integrated Development Environment is used to read all data in proper order for further processing through program structures. A high level of understanding of how to access and manipulate the underlying object structure of a CATIA parametric model is needed. The base of the program is obtained using the Macro Recorder to create code snippets that can be implemented into macros to reduce the time for development. The process involves the automatic creation of a parametric model using user inputs through the VBA programming interface. The user input in CATIA is running the macro and inputs the data file. The model is built from scratch, and all the changes are performed in the Excel data file. Such automation is useful when the part's geometry is defined by physical equations. Design steps from 1 to 4 are performed in Excel, from 5 to 9 executing a macro within CATIA, and remaining steps with GUI in CATIA. 1. set up NACA profile generator; 2. generate points (X,Y) of profile on different layers (Z) from equations; 3. rearrange points to start and finished numbering with outlet edge; 4. generate names of points incrementally; 5. pass point coordinates to CATIA; 6. generate airfoil sections with spline, starting and finishing in outlet edge with no closing - the result is a sharp edge needed for numerical analyses; 7. generate guidelines, from points with the same name on different layers; 8. loft sections using guidelines; 9. close surface into solid body; 10. add the blade into the assembly with prepared shaft and perform circular pattern; 11. save result in IGES format for numerical analyses. Section geometry is based on the four-digit NACA series, the parametric equations for which are given in Figure 3. These equations are set up on an Excel sheet to create points on a section profile based on a NACA number (4415 in the example in Figure 3) and a cord length. Profile point coordinates are updated on the spreadsheet based on each cord length and passed back to CATIA to create splines on each sketch layer. Sketches are lofted (multi-section solid or surface feature in CATIA) using the leading and trailing edges as guide curves to create the final model of the airfoil. In Figure 4, versatile models that were generated using automated processes are shown. Figure 3: Parametric equations for generating NACA airfoils [3]; generated and manipulated data of wind turbine blade; wind turbine geometrical model. Figure 4: Examples of generated wind turbines. 5 CONCLUSIONS This paper presents CAD automation involving skills developed for VBA programming in Excel to encompass the use of VBA in creating and manipulating parametric 3D CAD models in CATIA. It is important to note that there is a strong correlation between programming and 3D parametric modelling. For example, variables in programming are conceptually similar to the parameters that control size and shape in a CAD model. The process of the generation of geometry of wind turbine blades is adapted for the next step in the validation process, the CFD. The experience in numerical analyses preparation, i.e. structural meshing of the domain, creating the surfaces needed to implement boundary and initial conditions, is integrated into the macro-programme. Process automation works for versatile types of airfoil designs, including wings, spoilers, propellers and wind turbines. References [1] Wikimedia Foundation, Inc., Jun 2012. [Online]. http://en.wikipedia.org/wiki/Wind farm [2] T. Burton, D. Sharpe, N. Jenkins and E. Bossanyi, Wind Energy Handbook, 1 ed., West Sussex: John Wiley & Sons, 2001. [3] I. Abbott, A. Von Doenhoff, Theory of Wing Sections. Dover, 1959. [4] J. Kelly, Excel 2003 VBA Programming. Wiley Publishing (Visual Series), 2005. [5] J.C. Menezes, M.V. Donadon, Performance prediction of wind turbine blades, 12th Pan-American Congress of Applied Mechanics , Trinidad, 2012 [6] H.V. Mahawadiwar, V.D. Dhopte, P.S. Thakare, R.D. Askhedkar, CFD Analysis Of Wind Turbine Blade, International Journal of Engineering Research and Applications (IJERA), Vol. 2, Issue 3, May-June 2012, pp. 3188-3194 Journal of H I Energy j Technology Instructions to Authors http://www.fe.um.si/en/jet.html AUTHOR INSTRUCTIONS (MAIN TITLE) SLOVENIAN TITLE Authors, Corresponding authorm Key words: (Up to 10 keywords) Abstract Abstract should be up to 500 words long, with no pictures, photos, equations, tables, only text. Povzetek (In Slovenian language) Submission of Manuscripts: All manuscripts must be submitted in English by e-mail to the editorial office at jet@uni-mb.si to ensure fast processing. Instructions for authors are also available online at http://www.fe.um.si/en/jet.html. Preparation of manuscripts: Manuscripts must be typed in English in prescribed journal form (Word editor). A Word template is available at the Journal Home page. A title page consists of the main title in the English and Slovenian languages; the author(s) name(s) as well as the address, affiliation, E-mail address, telephone and fax numbers of author(s). Corresponding author must be indicated. Main title: should be centred and written with capital letters (ARIAL bold 18 pt), in first paragraph in English language, in second paragraph in Slovenian language. Key words: A list of 3 up to 6 key words is essential for indexing purposes. (CALIBRI 10pt) Abstract: Abstract should be up to 500 words long, with no pictures, photos, equations, tables, -text only. Povzetek: - Abstract in Slovenian language. m Corresponding author and other authors: Title, Name and Surname, Tel.: +XXX x xxx xxx, Fax: +XXX x xxx xxx, Mailing address: xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx, E-mail address: email@xxx.xx Main text should be structured logically in chapters, sections and sub-sections. Type of letters is Calibri, 10pt, full justified. Units and abbreviations: Required are SI units. Abbreviations must be given in text when firstl mentioned. Proofreading: The proof will be send by e-mail to the corresponding author, who is required to make their proof corrections on a print-out of the article in pdf format. The corresponding author is responsible to introduce corrections of data in the paper. The Editors are not responsible for damage or loss of manuscripts submitted. Contributors are advised to keep copies of their manuscript, illustrations and all other materials. The statements, opinions and data contained in this publication are solely those of the individual authors and not of the publisher and the Editors. Neither the publisher nor the Editors can accept any legal responsibility for errors that could appear during the process. Copyright: Submissions of a publication article implies transfer of the copyright from the author(s) to the publisher upon acceptance of the paper. Accepted papers become the permanent property of "Journal of Energy Technology". All articles published in this journal are protected by copyright, which covers the exclusive rights to reproduce and distribute the article as well as all translation rights. No material can be published without written permission of the publisher. Chapter examples: 1 MAIN CHAPTER (Arial bold, 12pt, after paragraph 6pt space) 1.1 Section (Arial bold, 11 pt, after paragraph 6pt space) 1.1.1 Sub-section (Arial bold, 10pt, after paragraph 6pt space) Example of Equation (lined 2 cm from left margin, equation number in normal brackets (section.equation number), lined right margin, paragraph space 6pt before in after line): (1.1) Tables should have a legend that includes the title of the table at the top of the table. Each table should be cited in the text. Table legend example: Table 1: Name of the table (centred, on top of the table) Figures and images should be labelled sequentially numbered (Arabic numbers) and cited in the text - Fig. 1 or Figure 1. The legend should be below the image, picture, photo or drawing. Figure legend example: Figure 1: Name of the figure (centred, on bottom of image, photo, or drawing) References [1] Name. Surname: Title, Publisher, p.p., Year of Publication Example of reference-1 citation: In text, Predin, [1], text continue. (Reference number order!) Authors names and surname (centred, Calibri 10pt, Italic)s JET Vol. 6 (2013) Issue 1 9771855574008