공학박사학위논문 수평축풍력발전용터빈블레이드최적설계및공력성능해석에관한연구 A Study on the Optimum Blade Design and the Aerodynamic Performance Analysis for the Horizontal Axis Wind Turbines 지도교수 : 이영호 2005 년 8 월 한국해양대학교대학원 기계공학과김범석
공학박사학위논문 수평축풍력발전용터빈블레이드최적설계및공력성능해석에관한연구 A Study on the Optimum Blade Design and the Aerodynamic Performance Analysis for the Horizontal Axis Wind Turbines 지도교수 : 이영호 2005 년 8 월 한국해양대학교대학원 기계공학과김범석
본논문을김범석의공학박사학위논문으로인준함 위원장 : 공학박사 박권하 ( 인 ) 위원 : 공학박사 정재현 ( 인 ) 위원 : 공학박사 남청도 ( 인 ) 위원 : 공학박사 김유택 ( 인 ) 위원 : 공학박사 이영호 ( 인 ) 2005 년 7 월 4 일
목 차 Abstract Nomenclature 제 1 장서론 1 1.1 풍력발전현황 1. 1.2 연구동향 4 1.3 연구목적 6 1.4 현대식풍력발전용터빈 8 1.5 현대식풍력발전용터빈의구조 10 1.5.1 로터블레이드 14 1.5.2 동력전달장치 14 1.5.3 발전기 15 1.5.4 너셀 15 1.5.5 타워 16 제 2 장수평축풍력터빈의공기역학 17 2.1 Actuator disk 이론 17 2.1.1 운동량이론 19 2.1.2 동력계수 21 2.1.3 최대동력계수 21 2.1.4 축추력계수 22 2.2 Rotor disk 이론 24 2.2.1 각운동량이론 28 2.2.2 최대동력 31
2.3 날개요소이론 33 2.4 날개요소운동량이론 37 제 3 장수평축풍력터빈블레이드최적설계 43 3.1 중ᆞ 대형 (1MW) 로터블레이드최적설계 43 3.1.1 설계풍속의결정 43 3.1.2 로터블레이드직경및정격회전수의결정 44 3.1.3 날개끝손실계수의보정 45 3.1.4 새로운흐름유도계수들의결정 48 3.1.5 무차원현의길이결정 50 3.1.6 입구유입유동각및비틀림각의결정 53 3.1.7 로터블레이드익형선정조건 55 3.1.8 2차원익형공력특성의예측 58 3.2 중형로터블레이드 (100kW) 최적설계 66 3.2.1 로터블레이드직경및정격회전수의결정 66 3.2.2 날개끝손실계수의보정 66 3.2.3 새로운흐름유도계수들의결정 68 3.2.4 무차원현의길이결정 70 3.2.5 비틀림각의결정 71 3.2.6 2차원익형공력특성 72 3.3 소형로터블레이드 (20kW) 최적설계 78 3.3.1 설계변수들의결정 78 제 4 장성능해석소프트웨어개발및평가 88 4.1 소프트웨어개발 88 4.2 실속후공력특성의예측 95
4.3 FIL-1000 성능평가 (1MW turbine) 97 4.4 FIL-100 성능평가 (100kW turbine) 103 4.5 FIL-20 성능평가 (20kW turbine) 111 제 5 장 CFD를이용한풍력터빈전산해석 121 5.1 수치해석기법 122 5.1.1지배방정식 124 5.1.2 이산화방법 125 5.1.3 난류모델링 129 5.2 계산조건 (FIL-1000) 132 5.3 계산격자및경계조건 (FIL-1000) 134 5.4 결과및고찰 (FIL-1000) 138 5.4.1 블레이드표면유선 138 5.4.2 블레이드국부단면유동특성 143 5.4.3 출력특성 151 5.5 계산조건 (FIL-20) 153 5.6 계산격자및경계조건 (FIL-20) 155 5.7 결과및고찰 159 5.7.1 블레이드표면유선 159 5.7.2 블레이드 3차원유동특성 163 5.7.3 출력특성 173 제 6 장결론 175 참고문헌 177 감사의글
A Study on the Optimum Design of the Horizontal Axis Wind Turbine and Its Aerodynamic Performance Analysis BeomSeok, KIM Department of Mechanical Engineering Graduate School of Korea Maritime University Abstract Wind energy has been emerged as one of important alternative energy from nature sources. Especially in part of Europe, for example, Denmark and Germany, wind has become a definite choice of the national power generation policy to cope with energy problem in terms of environment and energy shortage. In modern wind power system of large capacity above 1MW, horizontal axis wind turbine(hawt) is a common type. And, the optimum design of wind turbine to guarantee excellent power performance and its reliability in structure and longevity is a key technology in wind industry. In this study, mathematical expressions based upon the conventional blade element momentum theory(bemt) applying to turbine blade design was first analyzed systematically to secure accurate power prediction from the basic aerodynamic parameters such as lift and drag
coefficients, Prandtl's tip loss coefficient, tangential and axial flow induction factors and twist angle. X-FOIL open software was used to acquire lift and drag coefficients of the 2-D airfoils used in power prediction procedure. Three kinds of turbine rotor blade(named as FIL series, FIL-1000:1MW, FIL-100:100KW and FIL-20:20KW) were selected as examples of the optimum blade design and power prediction. Especially, in FIL-1000 case, three kinds of airfoil(ffa, DU and NACA series) were combined to produce maximum aerodynamic performance around blade tip and strong structure around blade root. Furthermore, optimum blade design and its power prediction software, named as "POSEIDON", was also developed from the BEMT discussed in the present study. Several input data are necessary and power curve is acquired easily to show the overall characteristics of the suggested wind turbine performance. Three kinds of blades(fil-1000, FIL-100 and FIL- 20) were again selected to apply the suggested software. In case of FIL- 1000, power coefficient was 0.47 at TSR=7. At FIL-100, power coefficients were compared between the experimental aerodynamic data and X-FOIL prediction data, resulting in fairly good agreement. In case of FIL-20, single blade shape(nrel S809) was adopted and comparison between experimental aerodynamic parameters and X-FOIL data was made with good agreement in power prediction. Finally, systematic computational fluid dynamics(cfd) analysis by commercial code CFX ver.5.7.1 was performed to study the detailed flow characteristics upon the blade surface and wake region. Turbulence model, k-ω SST was selected to guarantee stall phenomena, one of 3-D separation flow occurring in wind turbine. ICEM-CFD, reliable grid
generation commercial software was also adapted to secure good quality of grid generation necessary for the reliable CFD simulation. Three kinds of turbine were again selected to represent the complex 3-D stall flows appearing in the blade surface at various TSR condition. Centrifugal acceleration and pressure difference between hub and tip played a major role to govern the stall phenomena. Power prediction from the CFD data was also made and compared with BEMT prediction in case of FIL-20, with satisfactory agreement. In the future, aero-elastic analysis of the blade will be performed by the two-way FSI method software such as CFX-ANSYS and the more reliable blade design procedure and its performance prediction will be established.
Nomenclature A Swept area Axial flow induction factor Tangential flow induction factor C Chord length C L Lift coefficient C D Drag coefficient C L /C D Lift to drag ratio C P Power coefficient C T Thrust coefficient F Force, Blending factor f (µ) Local tip loss coefficient H Total enthalpy h Height N Blade number P Power, Pressure
Q Driven torque R Radius r Local radius S Strain rate U Wind velocity V Cut_in Cut-in wind velocity V Cut_out Cut-out wind velocity V D Design wind velocity V R Rated wind velocity W Wind velocity relative to a point on blade Greek letters α Angle of attack β Inclination of local blade chord η Efficiency of power train and generator θ Blade twist angle λ Tip Speed Ratio (TSR)
λ r Local tip speed ratio µ Non-dimensional local position, Viscosity ρ Air-density σ r Local solidity τ Stress tensor Ω Angular velocity
제 1 장 서론 1.1 풍력발전현황 1970년대에너지위기를겪은이후, 선진국을중심으로유한자원인화석에너지를대신할새로운에너지원의확보에많은관심을기울이고있으며, 정부차원의과감한지원을통해각국가별로다양한대체에너지원에대한연구개발프로젝트가진행되고있다. 최근의이라크전쟁이후, 전세계는배럴당 60 달러에육박하는고유가시대를맞이하면서자국의원활한에너지공급을위한무한에너지경쟁시대에접하게되었으며, 에너지자원의안정적인확보가국가의안보와직결되는상황에직면하게되었다. [1] 우리나라는에너지수입비중이 97.3%(439억달러, 2003년기준 ) [2] 로서세계최고수준의에너지소비국가임에도불구하고, 에너지안보및환경문제에대한인식의부재로인한산업경제활동위축및국가안보등에대한부담이가시화되고있다. 또한, 교토의정서 2차공약기간 (2013년 ~2017년 ) 중온실가스감축의무부담이가시화될전망이고, 이렇게되면국내의산업경제활동은장기적으로심각한침체기를맞이할것이라예상된다. 따라서, 종합적인국가에너지장기전략의수립을통한문제인식과적절한대응으로써, 정부주도하에과감한신ᆞ재생에너지분야에대한연구개발및보급정책의추진이필요한시점이다. 현재신ᆞ 재생에너지원으로써, 풍력, 연료전지, 태양광, 지열, 파력등다양한발전방식이제안되고있으나, 지난십수년간독일, 덴마크, 네덜란드, 영국등유럽국가들을중심으로요소기술개발및실증연구를거쳐이미보급경제성을평가받고있는풍력발전 (wind 1
power generation) 이가장현실적인대안이될수있다. 풍력발전은바람의운동에너지를이용하여블레이드 (blade) 를회전시키고이때발생하는회전력을구동장치계 (drive train) 를통해발전기로전달하여전기를생산하는방식을말한다. [3] 풍력발전은크게육상풍력발전 (on-shore) 과해상풍력발전 (offshore) 으로구분될수있으며, 설치면적의제한, 주변장애물에의한풍속감소, 소음등의문제로인하여현재점진적으로육상풍력발전단지에서해상풍력발전단지로전환되고있는추세이다. [4] 최근 5년간, 세계풍력발전시장의평균성장률은 31.7% 로써급격한보급확대가이루어지고있으며, 점차대형화, 단지화되어발전단가가 2003년기준으로 4 /kwh까지하락하여재래발전방식과비교하여충분한시장경쟁력을확보할수있게되었다. [5] Table 1.1에전세계풍력터빈설치현황및연도별설치용량변화율에관한자료를나타내었다. 2
Table 1.1 Wind turbine installation capacity growth rate in the world [6] Year Installed Increase Cumulative Increase (MW) (%) (MW) (%) 1998 2,597-10,153-1999 3,922 51 13,932 37 2000 4,495 15 18,449 32 2001 6,824 52 24,927 35 2002 7,227 6 32,037 29 2003 8,344 15 40,301 26 Averaged growth : 5 yrs 26.3 139,799 31.7 3
1.2 연구동향 최근전산해석기술의눈부신발달에힘입어선진연구기관을 중심으로 CFD 에의한로터블레이드 2 차원및 3 차원수치해석과 다양한풍력터빈로터블레이드용익형에대한수치해석적연구가 활발히진행되고있다. Wolfe [7] 는상용 Navier-Stokes solver 인 CFD-ACE 를이용해 2 차원 익형공력특성해석에관한연구를수행하였다. 그는익형의 CFD 해석시천이지점에대한정확한위치만파악하면부착류 (attached flow) 가형성되는조건에서정확한공력특성의예측이가능하다고 하였다. 또한, 일반적인상용코드의기본난류모델로적용되는 k 모델은 post-stall 영역에서의정확한공력특성예측이불가능하며, k ω 계열의난류모델이적합하다는결론을도출했다. Xu [8] 는 NREL Phase Ⅵ 로터블레이드를대상으로 in-house 코드에 의한수치해석을시도하였다. Sorensen [9] ε 등에의해 CFD 에의한 3 차원공기역학적특성파악을 위한연구가진행되었으며, solver 로써 k 적용한상용코드인 Ellipsys 3D 를이용하였다. Zuijlen [10] ω SST eddy viscosity model을 등에의해, adaptive 격자계를이용한 2 차원정상상태 조건에서의익형공력특성예측이수행되었다. 이들의연구에적용된난류모델은 1-equation 모델인 Spalart- Allmaras 모델이며, Navier-Stokes solver 로서 Hexstream 을사용하였다. 또한, 가장최근의 3 차원 CFD 해석결과로서 Sorensen [11] 등에의해 수행된결과가있으며, 상용 CFD 코드인 Edge, Ellipsys 3D 가적용 되었다. 이들의연구에서는 LES(Large Eddy Simulation) 수준의계산 정확도를확보할수있다고알려져있으며, 동일한조건에서 LES 해석시필요한계산격자수보다적은격자수가요구되는장점이 4
있는최신난류모델인 DES (Detached Eddy Simulation) 모델을 적용하였으며, post-stall 영역에서비교적정확한공력특성을예측 한다고알려져있는 k ω 계열의난류모델과비교하였다. 결과로써 로터블레이드 root 부근의극심한실속이발생하는영역에서의 공력특성예측성능에있어 DES 모델이 k 실험결과와더잘일치한다는결과를얻었다. 국내의경우, 국내최초로이영호 [12] ω 계열의난류모델보다 등에의해 Delft Technical University 의 T40/500 모델을대상으로 CFD 해석이수행되었으며, Navier-Stokes solver 로써 CFX-TASCflow 2.11.1 을이용하였다. 가장최근에발표된 CFD 연구결과에서는 BEMT(Blade Element Momentum Theory) 에의해자체적으로설계된 1MW(FIL-1000) 로터 블레이드를대상으로 3 차원유동해석및성능평가에관한연구가 수행되었고, Navier-Stokes solver 로써 CFX-5 를사용하였으며 3 차원 실속현상및후류 (wake) 의거동, 터빈출력예측성능등에관한결과 를얻었다. [13] 현재최신의고급전산해석기법인유체 ᆞ 구조연성해석 (Flow Structure Interaction, FSI) 을응용한관련연구가진행되고있으며, CFX- 5 와 ANSYS Multi-Physics 를 2-way FSI 기법으로연동하여로터 블레이드에작용하는유체력에의한진동문제및공력하중등의 해석에관한연구가진행중이다. 5
1.3 연구목적 풍력터빈구성요소중로터블레이드는바람의운동에너지를회전력으로변환하는핵심적인요소이며, 효율적인최적설계기법의확립이절실히요구되는요소기술이다. 그러나이론적최적설계에관한참고문헌의부재와선진국가에서핵심설계기술의공개를회피하는실정으로로터블레이드최적설계기법의확보가상당히어려운실정이다. 현재국내에서로터블레이드설계시 BLADED for Windows, YAWDYN 등과같은설계소프트웨어에의존하여수행하고있으나, 근본적인최적로터블레이드설계기술은확보되고있지않은실정이다. 이러한설계소프트웨어의활용은공인된소프트웨어에의한신뢰성있는설계데이터를확보할수있는장점이있지만, 설계알고리즘자체를이해하지못한상황에서더이상의발전을기대할수없으며, 향후다양한설계변수를적용한독자적인로터블레이드에대한최적설계를수행할수없다. 현재까지풍력터빈의성능평가작업은원형 (prototype) 을제작하여실운전조건하에서수행되는 field test, 대규모풍동실험을통한모형실험, BEMT(Blade Element Momentum Theory) 등에의한수치해석적방법에의존하고있다. 우선실험적기법은가장정확하고신뢰성있는데이터의확보라는큰장점이있음에도불구하고실험규모의대형화에따른제약과, 실제스케일의로터블레이드원형제작등에관한문제점으로소규모연구그룹이단독으로수행할수없는단점이있다. 해석적인성능예측방법으로써현재까지다양한성능예측모델이제안되고있으나, 일반적인설계과정에서 BEMT 를이용한성능예측이주로행해지고있으며, 풍력발전산업에서성능해석의표준 6
으로자리잡고있다. [14] 그러나, BEMT 법은정확한 2 차원익형의공력특성데이터를활용하면비교적신뢰성있는성능예측결과를도출할수있는반면, 로터블레이드주위의흐름, 블레이드표면의압력분포, 블레이드후류의유동특성등과같은어떠한유체역학적인문제도가시화할수없는단점이있다. 이를해결하기위해필드테스트를거치는방법과축소모델에의한풍동시험등과같은 2 가지방법이있으나, 필드테스트와축소모델에의한성능시험및가시화방법은상당한비용과시간을필요로하는대규모작업이므로보다효율적인방법이절실히요구된다. 2001 년이후부터 RISO, NASA, NREL 등의대규모국가연구소를중심으로 CFD 를이용한풍력발전용로터블레이드 3 차원유동해석및성능평가에관한연구가활발히수행되고있으며, 몇몇결과는상당히만족스러운결과를제시하고있다. 본연구의목적을정리하면다음과같다. 1) 날개끝손실모델을포함하는 BEMT 이론에의한최적로터 블레이드형상설계기법의확립 2) 효율적인 2 차원익형공력특성 ( 양력, 항력 ) 의예측 3) 로터블레이드설계및성능해석소프트웨어의개발 4) 소형 (20kW), 중형 (100kW), 대형 (1,000kW) 풍력터빈최적설계및 BEMT 에의한성능해석 5) CFD 를이용한로터블레이드 3 차원유동예측및성능평가 7
1.4 현대식풍력발전용터빈 풍력발전용터빈이란, 바람이가지는운동에너지를전기에너지로변환하는에너지변환장치를말한다. 이는과거바람의에너지를기계적에너지로변환하는전통적인기계장치인 windmill 의개념과는상반되는것이다. 전기에너지의생산역할을하는풍력발전기는일반송전망에연계되며, 송전시설에는배터리축전회로시스템, 거주지스케일의동력시스템, 다양한전기관련부가장치들이포함된다. 전체풍력발전용터빈의설치용량에따르면 10kW 미만의소형터빈이대다수를차지하고있으나, 전체풍력에너지생산량에따르면, 주로 500kW부터 2MW급의중ᆞ대형풍력발전용터빈들이주축을이루고있다. [15] 이러한중ᆞ대형터빈들은보통대규모송전시설들을필요로하게되며, 대다수의풍력터빈이미국과유럽등지에설치되어있다. 풍력발전용터빈의용도를이해하기위해터빈의작동개념에대한몇몇기본적인원리를이해하는것이중요하다. 현대식풍력발전용터빈에서실제바람에너지의변환과정은회전하는축 (shaft) 에작용하는전체토오크 (torque) 의발생원인이되는공기역학적힘 (aerodynamic force), 양력 (lift force) 의발생으로부터시작된다. 양력에의해발생하는토오크는기계적동력을발생시키고, 발생되는동력은발전기를통해전기에너지로재생산된다. 바람에너지는저장할수없으며, 차후재사용할수도없기때문에풍력터빈은바람에너지밀도변화에신속히대응할필요가있으며, 이러한점이대다수기타발전방식과의가장큰차이점에해당한다. 역사적으로밀과같은곡물을분쇄하기위한목적으로 windmill과같은장치가등장하였으며, 이러한장치들은현대에이르러그사용 8
목적에대한관점이변화되었다. 오늘날의광범위한송전망을통한전기에너지수송능력은풍력 발전의발전가능성을확대해주는중요한역할을하고있다. 9
1.5 현대식풍력발전용터빈의구조 오늘날대다수의상용풍력발전용터빈의설계형식은로터의회전축이지면에대해평행한구조를갖는수평축로터시스템 (horizontal rotor system) 을따르고있다. [16] 보통 HAWT(Horizontal Axis Wind Turbine) 로터는로터의방향 (upwind or downwind), 허브디자인 (rigid or teetering), 로터제어시스템 (pitch or stall), 블레이드개수 ( 보통 2~3개 ), 그리고바람의방향성에대한로터블레이드의추종방식 (free yaw or active yaw) 에따라구분된다. Fig. 1.1에전형적인 upwind 시스템과 downwind 시스템의개략도를나타내었다. 10
Wind direction Fig. 1.1 HAWT rotor configurations (upwind - left, downwind - right) 11
일반적인수평축풍력발전용터빈을구성하는하위시스템장치 구성도를 Fig. 1.2 에나타내었다. l 로터 : 블레이드와지지허브로구성되어있음 l 동력전달장치 : 풍력발전용터빈의회전부를포함하고있음보통축, 기어박스, 커플링, 기계적브레이크장치, 발전기등으로구성됨 l 너셀 (nacelle) : 메인프레임하우징, 베드플레이트, 요우장치등을포함 l 타워와기초부 l 제어장치 풍력발전용터빈블레이드설계및제작시다음과같은선택사항 을고려한다. l l l l l l l l l 블레이드개수 : 보통 2~3 개로터설치방향 : downwind or upwind 블레이드의재질, 제작공법, 형상허브설계 : rigid, teetering, hinged 출력제어 : 실속제어형, 가변피치제어형회전속도제어 : 고정식, 가변식풍향변화에대한추종성 : free yaw, active yaw 발전기 : 동기형, 유도형회전력전달 : 기어박스, 직접연결식 12
Fig. 1.2 Configuration of wind turbine components 13
1.5.1 로터블레이드로터는풍력발전용터빈에서허브 (hub) 와블레이드로구성되는요소이며, 성능에직접적인영향을미치는가장핵심적인요소이다. 현재몇몇 2 블레이드형태의 downwind 형식풍력터빈이운전되고있기도하지만, 대다수의현대식풍력터빈은 3 블레이드 upwind 형식으로설계및제작된다. 과거단일블레이드를장착한풍력터빈을제작하기도했으나밸런싱 (balancing), 소음 (noise), 진동 (vibration), 출력저하등의문제로인해더이상생산되고있지않다. [17] 대다수의중형풍력터빈은블레이드피치각을고정하고출력제어방법으로서실속제어방식 (stall control) 을사용한다. 이는특히덴마크에서주로사용하는방법이며, 몇몇미국제조업체에서는출력제어방식으로피치제어방식 (pitch control) 을적용하고있으며, 현재중형급이상의터빈에대해서는구조적으로실속제어방식에비해다소복잡하지만정확한출력제어가가능한피치제어식출력제어방법을적용하는추세이다. [18] 로터블레이드는간혹 wood/epoxy 적층판으로제작되기도하나, 대부분이주로 GRP(fiberglass reinforced plastics) 등의복합재료로제작된다. 1.5.2 동력전달장치동력전달장치는풍력터빈의회전부를포함하는로터에연결되는저속축 (low-speed shaft), 저속의회전속도를발전기의정격회전속도로증속시켜주는역할을담당하는기어박스그리고발전기에직접연결되는고속축 (high-speed shaft) 으로구성된다. 기타기계요소들로써지지베어링과, 다수의커플링, 브레이크시스템, 그리고발전기의회전부등이있다. 관련기계요소중기어박스는저속으로회전하는로터의회전속도 14
( 수십 rpm) 를일반발전기의구동에적합한회전속도 ( 수백혹은수천 rpm) 까지증속시키는역할을한다. 기어박스의종류로서병렬축 (parallel shaft) 형과유성기어형 (planetary) 의두가지형태가있다. 대형풍력터빈 (500kW급이상 ) 에서는가볍고소형인유성기어형기어박스를주로사용하며, 몇몇풍력터빈은저속발전기를사용함으로서기어박스를사용하지않고로터블레이드와발전기가직접연결되어구동하기도한다. 이러한풍력터빈용동력전달장치의설계는기계공학에서의전형적인동력전달장치설계법을따르고, 가변적인바람에의한하중의변화와직경이큰로터가회전할때발생하는동역학적문제등은동력전달장치구성요소에급격한하중의변화를초래하므로설계과정에서충분한검토가필요하다. 1.5.3 발전기거의모든풍력터빈들은유도형 (induction) 발전기를사용하거나, 동기식 (synchronous) 발전기를사용한다. 이러한설계방식은발전기들이전력계통에직접연계되어있을때발전기가항상일정한회전속도를유지해야하는문제를수반한다. 현재대다수의풍력터빈발전기로서유도형전동기가사용되며, 동기식발전기에비해내구성이좋고, 가격이낮으며전력계통에연결하기쉬운장점을가진다. [19] 1.5.4 너셀너셀 (nacelle) 은 풍력터빈의 하우징 (housing) 과 베드플레이트 (bed plate), 메인프레임 (main frame), 방향제어시스템 (yaw control system) 등을포함한다. 메인프레임은동력전달장치의장착및정확한고정 을위한장치이며, 너셀커버는외부환경으로부터주요기계요소들을 보호하는역할을한다. 15
방향제어시스템은로터축을항상바람이불어오는방향에일치하도록제어하는장치이다. 방향제어시스템을구성하는주된장치는메인프레임과타워를연결하는대형베어링이다. Active yaw 구동시스템은일반적으로단일혹은이중 yaw 모터를장착한 upwind 형수평축풍력발전기에적용되며, yaw 베어링에장착된불기어 (bull gear) 에대해피니언기어 (pinion gear) 가맞물려구동한다. [20] 이러한방향제어시스템은풍력터빈의너셀후방에장착된풍향센서신호에의해자동으로제어되며, 너셀의방향을제자리에고정하기위한목적으로사용되기도한다. Free yaw 시스템의경우통상 downwind 시스템에적용된다. 1.5.5 타워일반적으로타워는강관 (steel tube), 격자구조물 (lattice structure), 콘크리트등을사용하여 free standing 형식으로설계된다. 소형터빈의경우가이드역할을하는타워들이종종설치되기도한다. 타워의높이는보통로터직경의 1~1.5배의크기를가지며, 어떠한경우에도최소한 20m 높이를확보하도록설계된다. 타워의선정은설치지형의특성에따라상당한영향을받는다. 회전하는로터와타워는상호작용에의한공진을발생시킬수있는위험요소가있으므로, 타워의강성은풍력발전시스템에서상당히중요한고려대상이다. 또한, downwind 형터빈에서는타워에의해발생하는후류의영향인타워그림자효과 (tower shadow effect) 가발생하며, 이는불균일한출력의발생및소음, 진동등의문제를수반하므로신중히고려해야한다. [21] 16
제 2 장수평축풍력터빈의공기역학 2.1 Actuator disk 이론 Fig. 2.1에서디스크상류의유관 (stream tube) 은디스크단면적보다작고, 디스크하류는디스크보다큰단면적을갖는다. 모든지점에서질량유량 (mass flow rate) 이같기때문에하류로갈수록유관의단면적변화가증가한다. 주어진유관내의임의의 위치에서단면적을통과하는질량유량은 ρ AU 로표현되고, ρ 는 밀도 (density), A 는단면적, U 는유속을나타낸다. 질량보존의법칙에따라질량유량은유관내의어느단면에서측정되더라도항상일정한값을유지하므로식 (2.1) 과같은관계가성립한다. ρ A U ρ A U ρ A U = = disk disk wake wake (2.1) 하첨자 는로터디스크상류로부터무한한거리에위치함을의미하고 disk 는디스크단면에서의조건을의미하며 wake 는디스크단면후류로멀리떨어진곳에서의조건을나타낸다. 일반적으로 actuator disk에의해속도의변화가발생된다고생각되고, 디스크에의해유도되는흐름의유선방향속도성분을 au 라한다. 여기서 a 를축흐름유도계수 (axial flow induction factor) 또는유입변수라한다. 디스크단면을중심으로총유선방향속도는식 (2.2) 와같이정리된다. U = U au (2.2) disk 17
Fig. 2.1 Idealized flow through an actuator disc 18
2.1.1 운동량이론 디스크단면을통과하는기류의속도는가변적이며, 유속의 변화량은 U U wake 와같이표현된다. 운동량변화율은전체속도 변화량에질량유량을곱한형태이므로, 이들의관계는식 (2.3) 과 같이표현된다. 운동량변화율 = ( U U ) ρ A U (2.3) wake disk disk 운동량변화에의해발생되는힘은디스크단면을가로지르는 압력의변화로부터발생한다. 따라서, 디스크에작용하는힘은식 (2.4) 와같이표현될수있다. + ( P P ) A = ( U U ) ρ A U (1 a) (2.4) disk disk disk wake disk 디스크단면에서발생하는압력차를알기위해유관의상류와하류에각각 Bernoulli 방정식을적용한다. 흐름방향으로의전체에너지는정상상태조건하에서운동에너지, 정압에너지, 위치에너지로표현된다. 먼저디스크상류에대해 Bernoulli 방정식을적용하면, 식 (2.5) 와같다. 1 2 1 2 + ρ U + P + ρ gh = ρdisku disk + Pdisk + ρdisk ghdisk (2.5) 2 2 흐름이비압축성이고, 수평높이가일정하다고가정하면식 (2.5) 는 식 (2.6) 과같이표현된다. 19
1 2 1 2 + ρu + P = ρu disk + Pdisk (2.6) 2 2 후류에대해서도이와동일한방법을적용하면식 (2.7) 과같이 표현됨을알수있다. 1 2 1 2 ρu wake + P = ρu disk + Pdisk (2.7) 2 2 식 (2.6) 에서식 (2.7) 을소거하면식 (2.8) 과같이정리된다. 1 = (2.8) 2 + 2 2 ( Pdisk Pdisk ) ρ( U U wake) 식 (2.8) 을식 (2.4) 에대입하면식 (2.9) 와같이정리된다. 1 ( 2 2 ρ U U wake ) Ad ( U U wake ) ρ AdiskU (1 a ) 2 = (2.9) 그리고, U wake = (1 2 a) U (2.10) 즉, 축방향손실중절반은디스크상류에서발생하고절반은 디스크하류에서발생한다. 20
2.1.2 동력계수 흐름에의해디스크에발생하는힘의변화는식 (2.4) 로부터구해질 수있으며, 식 (2.11) 과같다. F = P P A = A U a a (2.11) + 2 ( disk disk ) disk 2 ρ disk (1 ) 이힘은로터디스크면에집중되어발생하는힘이므로힘에의해 행해진일률은 FU disk 동력은식 (2.12) 와같이생각될수있다. 로표현된다. 따라서, 바람으로부터추출되는 Power = FU = 2 ρ A U a(1 a) (2.12) disk disk 3 2 따라서, 동력계수는식 (2.13) 과같이정의될수있다. C P Power = (2.13) 1 3 ρu Adisk 2 2.1.3 최대동력계수 바람으로부터이론적으로추출할수있는최대동력계수 C P,max 를구하기위해서는디스크단면을가로지르는바람의축방향손실 이없다는조건을만족해야하므로, 식 (2.14) 와같이표현될수있다. dc P 4(1 a)(1 3 a) 0 da = = (2.14) 21
식 (2.14) 의해를구하면축흐름유도계수는다음과같다. a = 1, 0.333 여기서, 축흐름유도계수값이 1 인것은물리적의미가없으므로 0.333 만해로취할수있다. 이를식 (2.14) 에대입하고정리하면식 (2.15) 를얻는다. 16 C P,max = = 0.593 (2.15) 27 따라서, 이론적으로추출할수있는풍력터빈의최대동력계수는 0.593을넘을수없으며, 이는독일의항공역학자인 Albert Betz [22] 에의해유도되었으며, 실제로 Betz의한계값을넘는풍력터빈을설계할수없다. 이를 Betz limit라한다. 2.1.4 축추력계수 식 (2.11) 은축추력계수 차원화될수있다. CT 를계산하기위해식 (2.16) 과같이무 C T Power = (2.16) 1 2 ρu Adisk 2 C = 4 a(1 a) (2.17) T 여기서, a 0.5 인경우후류속도 (1 2 a) U wake 가 0이되거나음의값을가질수도있기때문에문제가발생한다. 이러한경우더 22
이상운동량이론은적용될수없으며다음장에서논의될날개 요소 - 운동량이론에서와같은실험적보정이반드시필요하게된다. 23
2.2 Rotor disk 이론 바람으로부터유용한에너지를얻기위해서는보다효율적인특정터빈의설계가중요하다. 대다수의풍력발전용터빈들은각속도 Ω 로회전하는다수의블레이드로구성된로터를장착하며, 바람이불어오는방향에대해평행하게위치한다. 블레이드는바람의운동에너지에의해발생되는공기역학적힘에의해회전력을얻게되고후류의운동량손실에의한영향으로로터전 후단면에압력차가발생하게된다. 축방향운동량손실은바람의운동에너지손실을의미하며, 손실에너지는로터에의해흡수되어결과적으로동력전달장치저속축의회전력으로변환되게된다. 이때발생하는축회전력은증속기를통해고속축에부착되어있는발전기를구동시키게되며결국, 로터에서발생하는공기역학적힘, 토오크가전기적에너지로변환된다. 보다효율적인회전력의발생을위한로터의공기역학적최적설계부분은제 3 장에서자세히논하기로한다. 바람은로터디스크면을통과하면서로터의회전력으로변환되고, 동시에반대방향의회전력이로터디스크에작용하게된다. 이와같은반작용토오크에의한영향으로바람이로터하류로빠져나갈때반대방향의회전성분을가지게된다. 이때바람은각운동량을갖게되고, 로터디스크후류에서바람의단일미소입자에축방향속도성분과함께로터의회전면에대한접선방향으로의속도성분이동시에존재한다. 바람이로터블레이드를빠져나가면서접선방향속도성분을갖는다는것은로터후류에서압력감소에의한손실상쇄효과로써운동에너지의증가를의미한다. Fig. 2.2에로터디스크를통과하는입자의궤적에대해간략히 24
나타내었다. 로터디스크로유입되는유동은회전운동을하지않지만, 로터 디스크로부터빠져나오는흐름은회전성분을갖게되고이회전 성분은로터후류방향으로진행하는흐름속에서일정한값으로 존재하게된다. 접선방향속도의변화는접선방향흐름유도계수인 a 의항으로 표현된다. 일반적으로로터디스크상류에서의접선방향속도성분은없다. 로터디스크하류에서의접선방향속도성분은 2Ωra 이며, 디스크 중앙에서의접선속도변화는허브로부터임의의반경까지의거리를 r 이라할때 Ωra 으로표현된다. 앞서언급한바와같이접선방향 속도성분은토오크의반작용에의해발생하므로접선방향속도성분 의방향성은로터의회전방향과반대이다. Fig. 2.3 은블레이드사이에서서서히가속되는접선방향흐름의 가속성분을나타내고있다. 25
r Wind particle trajectory Ω Rotor disk plane Fig. 2.2 Concept of a wind particle trajectory through the rotor disk plane 26
1 (2 ) 2 Pdisk ρ a Ωr 2 2a Ωr P + disk U (1 a) U (1 a) a Ωr U (1 a) Fig. 2.3 Tangential velocity accelerations through the rotor disk 27
2.2.1 각운동량이론 접선방향속도와축방향속도는블레이드반경방향위치마다 모두다르다. 이두가지유도속도성분에대한검토를위해로터 디스크가반경 r 을가지고반경방향미소길이 δ r 을가지는환형 고리라고생각한다. 환형고리에작용하는로터토오크의증가에의해 접선방향속도성분은가속되는반면, 환형고리에작용하는축방향 힘은축방향속도를감소시키는결과를초래하게된다. 전체로터디스크는무수히많은 δ r 의적분형태로표현되며, 각각의미소환형링은상호간의운동량전달이전혀없고, 단지링을통과하는흐름에의해서만운동량의전달이이루어지는독립적인프로세스를갖는다고가정한다. [23] 미소환형링에작용하는토오크는환형링을통과하는흐름의각운동량변화율과동일하며식 (2.18) 과같이표현된다. δq ρδ A U a a r 2 = disk (1 )2Ω (2.18) 여기서, δ Adisk 는환형링의면적을의미한다. 로터축에대한구동토크또한 증분은식 (2.19) 와같이표현된다. δq 로표현되며로터축동력의 δ P = δqω (2.19) 따라서, 바람으로부터추출할수있는전체동력에너지는식 (2.12) 에의해주어진축운동량변화율에의해결정될수있다. δ P ρδ A U a a 3 2 = 2 disk (1 ) (2.20) 28
축추력에의해유도된동력과토오크에의해유도된동력을서로 같다고두면식 (2..21) 과같다. 2 ρδ A U a(1 a) = A U (1 a)2ω a r (2.21) disk 3 2 2 2 ρδ disk 정리하면, 식 (2.22) 와같다. U a a = Ω a r (2.22) 2 2 2 (1 ) 여기서, Ωr 은회전하는환형링의접선속도이므로국부속도비를 식 (2.23) 과같이정의한다. λ = Ω r / U (2.23) r 디스크끝단에서 r = R 이므로, 날개끝속도비 (tip speed ratio, TSR) 는식 (2.24) 와같이정의된다. λ = Ω R / U (2.24) 따라서, 식 (2.22) 는식 (2.25) 와같이표현된다. (1 ) = λ r (2.25) 2 a a a 환형링의면적은 2π rδ r 이므로축동력의증분은식 (2.18) 을 적용하면식 (2.26) 과같이표현된다. 29
1 δ P = dqω = ρu 2π rdr 4 a (1 a) λ 2 3 2 r (2.26) 식 (2.26) 에서괄호내의항은환형고리를통과하는동력의변화를의미하고, 괄호밖의항들은바람으로부터동력을추출하는블레이드요소의효율을의미한다. 따라서, 로터블레이드의효율은식 (2.27) 과같이정리된다. η = 4 a (1 a) λ (2.27) r 2 r 이를동력계수의형태로표현하면식 (2.28) 과같이정리된다. d C dr P 3 2 2 4 πρu (1 a) a λr r 8(1 a) a λr r = = 2 1 3 2 ρu π R R 2 d C P dµ 8(1 a ) a 2 3 = λ µ (2.28) 여기서, µ = r / R 이다. 식 (2.28) 을적분하면날개끝속도비의변화 에따른 a, a 의반경방향변화율과로터디스크전체동력계수를 계산할수있다. 30
2.2.2 최대동력 최대출력을발생시킬수있는 a, a 은식 (2.27) 을각각의변수에 대해양변을미분형태로취하면얻을수있으며그결과값을 0 이라 둔다. 식 (2.27) 을 a 에대하여미분하면, a λ r = 이되고 a 에대해 2 4 0 미분하면 2 4(1 ) 0 a λ r = 이된다. 두식을같다고놓고정리하면식 (2.29) 와같다. 1 a = 0 a (2.29) 식 (2.25) 도위와같은방법으로표현하면식 (2.30) 과같다. 2 d λr = da 1 2a (2.30) 식 (2.29) 와식 (2.30) 을같다고놓고정리하면식 (2.31) 과같다. λ = (1 )(1 2 ) (2.31) 2 a a a r 식 (2.31) 을 a 에대해정리하면식 (2.32) 와같다. a(1 a) a = (2.32) 2 2 λ µ 식 (2.32) 를식 (2.31) 에대입하여정리하면 a 를구할수있으며, 31
최종적인 a, a 은식 (2.33) 과같이표현된다. 1 a(1 a) a =, a = (2.33) 2 2 3 λ µ 식 (2.33) 에서축흐름유도계수 a 는운동량이론에서와같은 1 3 로표현되고, 이는전체디스크에대해항상균일한값이다. 그러나, 회전흐름유도계수 a 은반경방향에따라값이변한다. 식 (2.28) 로부터유도되는최대동력계수는식 (2.34) 와같다. 1 2 3 CP = a a d 8(1 ) λ µ µ (2.34) 0 식 (2.34) 에식 (2.33) 을대입해적분하면식 (2.35) 와같다. 2 16 CP = 4 a(1 a) = (2.35) 27 이값은운동량이론으로부터유도된최대출력계수와정확히일치 한다. 32
2.3 날개요소이론 풍력터빈로터블레이드에서스팬 (span) 방향반경 r 과길이 δ r 에 대한공기역학적양 항력은환형면적을통과하는모든바람의축 운동량, 각운동량의변화를초래하고, 후류내회전속도변화와 관련된압력강하의원인이되는힘을발생시킨다. 로터로접근하는바람의회전성분은없으므로, 후류의회전성분에 의해로터후방에서계단형압력강하가나타나게된다. 로터후방 으로멀리떨어진곳에서도후류는회전성분을가지고있으므로, 압력강하는여전히존재하나, 축운동량변화의원인이되지는않는 다. [24] 블레이드요소에작용하는힘은요소의단면으로유입되는축 방향과회전방향속도성분에대한합성속도가블레이드전연과 이루는각도인받음각을이용한 2 차원익형의공력특성에의해계산 될수있다는가정을한다. 이때스팬방향으로진행하는속도성분과 3 차원속도성분은무시한다. 각각의블레이드반경방향위치에서속도성분은풍속, 흐름유도 계수, 로터회전속도의항으로설명되며받음각을결정할수있다. 익형의공력특성으로대표되는양력계수와항력계수 ( C, C ) 의 받음각변화에따른변화량을알면주어진 a, a 값에서로터에 작용하는힘을계산할수있다. 블레이드현의길이가 C, 팁반경이 R, 블레이드수가 N 인로터 블레이드를생각해보면, 높은효율을갖는블레이드를설계하기 위해서스팬방향으로현의길이, 피치각등이적절하게변화하게될 것이다. 이블레이드가각속도 Ω 로회전하고있으며, 유입풍속을 U 라한다. L D Fig. 2.4(a) 에서와같이블레이드요소의접선속도는 Ωr 이며, 33
후류의접선속도는 Ωra 속도의합이 (1 + a ) Ωr 과같다는의미이다. 이다. 이는, 블레이드요소에서전체유동 Fig. 2.4(b) 는반경 r 에서블레이드현의길이에관계되는모든 속도성분과힘의관계를나타내고있다. Fig. 2.4(b) 로부터블레이드에작용하는합속도는식 (2.36) 과같이 표현된다. W = U a + Ω r + a (2.36) 2 2 2 2 2 (1 ) (1 ) 합속도는로터의회전면에대해각도 φ 를갖는다. U (1 a) sinφ =, W Ω r(1 + a ) cosφ = (2.37) W 받음각 α 는식 (2.38) 과같이정의된다.. α = φ β (2.38) 블레이드요소길이 δ r 에작용하는양력은식 (2.39) 와같다. 1 2 δ L = ρw CCLδ r (2.39) 2 항력은식 (2.40) 과같다. 1 2 δ D = ρw CCDδ r (2.40) 2 34
(a) Blade element sweeps out an annular ring Fig. 2.4 Relationship between blade element forces and velocities (continued) 35
(b) Blade element velocities and forces Fig. 2.4 Relationship between blade element forces and velocities 36
2.4 날개요소운동량이론 날개요소이론의기본적인가정은블레이드요소에서발생하는힘이, 요소의면적을통과하는바람의운동량변화를초래하는유일한성분이라는것이다. 따라서, 반경방향요소들사이에서의상호작용은없다고가정할수있으며, 축방향흐름계수의반경방향변화는없다고생각할수있다. 실제로는, 축흐름유도계수가반경방향으로항상일정한값을유지하지는않지만 1934년 Lock [25] 의프로펠러실험에의하면, 이와같은반경방향에대한독립성은타당한가정이라고생각된다. N 개의블레이드에작용하는축방향공기역학적힘은식 (2.41) 과같다. 1 2 δl cosφ + δdsin φ = ρw NC( CL cosφ + CD sin φ) δr (2.41) 2 요소면적을통과하는기류의축방향운동량변화율은 mu & 이 므로블레이드를통과하는바람의크기는 ( U U wake ) 이고, U wake = (1 2 a) U 이므로환형국부회전면적을통과하는바람의축 방향운동량변화율은식 (2.42) 와같다. ρ π δ = πρ δ (2.42) 2 U (1 a)2 r r 2aU 4 u a(1 a) r r P a r 2 drop, wake = 1/ 2(2 Ω ) (2.43) 후류의회전성분에의해발생되는압력강하를식 (2.43) 과같다고 37
하면, 축방향으로작용하는힘은식 (2.44) 와같다. = Ω (2.44) 2 F 1/ 2(2 a r) 2π rδ r 그러므로, 공기역학적인축방향힘을나타내는식 (2.41) 은바람이통과할때축방향운동량의변화와로터블레이드후류압력강하에의해발생되는힘의합으로표현될수있으며, 이를수식으로표현하면식 (2.45) 와같다. 2 2 2 1/ 2 ( L cos + D sin ) = 4 [ (1 ) + ( Ω ) ] ρw Nc C φ C φ δr πρ U a a a r rδr (2.45) 식 (2.45) 를정리하면식 (2.46) 과같이표현할수있다. W U C 2 2 N ( C cos sin ) 8 [ (1 ) ( ) ] 2 L CD a a a R φ + φ = π + λµ µ (2.46) 또한, 블레이드는축방향힘외에도접선방향회전력즉, 토오크 에의한힘을발생하며식 (2.47) 과같이정리된다. 2 ( δ sinφ δ cos φ) 1/ 2 ρ ( L sinφ D cos φ) L D r = W NC C C rdr (2.47) 각운동량변화율은 m& U 이므로, 식 (2.48) 과같다. tangential 38
ρ Ω π δ = πρ Ω δ (2.48) 2 U (1 a) r2a r2 r r 4 U ra (1 a) r r 정리하면, 식 (2.49) 와같이표현된다. W u c 2 2 N ( C sin cos ) 8 (1 ) 2 L φ CD φ πλµ a a R = (2.49) 축방향운동량변화에의해유도된식 (2.46) 과접선방향회전력으로부터유도된식 (2.49) 의양력계수와항력계수에관계되는항을식 (2.50) 과같이치환하면, 식 (2.51), 식 (2.52) 와같은단순화된식으로표현할수있다. 최종적으로식 (2.51), 식 (2.52) 를통해반복계산법과블레이드각각의국부위치에해당하는 2차원익형의공력특성데이터를이용하여축방향흐름유도계수 a 와회전방향흐름유도계수 a 을각각계산할수있다. C cosφ + C sinφ = C, C sinφ C cosφ = C (2.50) L D x L D y a σ r σ = 1 a 4sinφ 4sin φ r 2 [ C ] 2 x C 2 y (2.51) a σ rcy = 1+ a 4sinφ cosφ (2.52) 여기서, 블레이드솔리디티 (solidity) σ 는로터블레이드의전체 39
회전면적 (swept area) 중로터블레이드면적이차지하는비를말하며, 로터블레이드성능및제작단가에영향을미치는변수이다. σr NC N C = = (2.53) 2π r 2πµ R 주어진식 (2.51) 과 (2.52) 에대한해를구하면적절한블레이드의기하학적형상을설계할수있으며, 동력, 토오크등의관계값들을계산할수있다. 1974년 Wilson과 Lissaman [26] 은익형의후연 (trailing edge) 으로부터협소후류영역에대해서만항력에의한속도감소영향이존재하므로식 (2.51) 과식 (2.52) 에서항력항은무시해도된다는주장을하였다. 더욱이, 항력에의한후류속도감소는단지후류발생특징에제한될뿐로터블레이드상류영역의속도감소에는어떠한영향도미치지못한다. 또한, 표면부착류 (attached flow) 에있어항력은표면마찰에만의존하는요소이므로로터블레이드단면에서발생하는압력강하에대한영향을미치질못한다. 하지만로터블레이드에서표면부착류가형성되지못하고실속상태에놓여있을때항력계수는로터블레이드압력강하에영향을미치는인자로써작용하게된다. 이러한항력계수항의무시에관한논의는로터블레이드초기설계단계에서적용되는수식의단순화및계산과정의단순화를위해많은도움이되는부분이며, 또한최적설계조건에서의로터블레이드표면흐름은부착류가형성되는특징을가지므로초기형상설계시항력항을무시하고설계를수행하는방법이효과적이다. 기본설계를바탕으로한로터블레이드형상은차후항력항이포함된수식을통해최적화되며그성능평가가진행된다. 40
로터의토오크와동력을계산하기위해서는식 (2.51) 과식 (2.52) 의해를구함으로써얻어지는흐름유도계수값이필요하다. 흐름유도계수값의계산은익형의공력특성데이터가받음각의변화에대해비선형함수로표현되므로반복계산법에의해수행될수밖에없다. 식 (2.49) 로부터반경방향길이 δ r 을가지는블레이드요소에서 작용하는토오크를계산할수있다. δ = πρ Ω δ (2.54) 2 Q 4 U ra (1 a) r r 만약, 흐름유도계수들을구할때항력에의한영향을고려하지 않았다면식 (2.54) 는항력에대한관계식이포함된식 (2.55) 와같이 수정되어야한다. 2 1 2 δq = 4 πρu Ωra (1 a) r δr ρw NCCD cosφδr (2.55) 2 따라서, 전체토오크는식 (2.56) 에의해결정된다. C 1 R N 2 3 2 W Q U R 8 a (1 a) R = ρ π λ µ µ C (1 ) 2 0 D + a d µ (2.56) U π 로터의동력은 P 된다. = QΩ 로표현되므로동력계수는식 (2.57) 로계산 41
C P P = (2.57) 1 3 2 ρu π R 2 BEMT에의한설계및성능평가는초기에 로터블레이드반경방향으로진행하는흐름은없다 라는균일순환조건의가정하에진행된다. 즉, 로터블레이드상류에서유입되는흐름의축방향성분은축방향으로항상일정하게유지된다는것이다. 그러나, 불균일순환조건에서는로터블레이드반경방향의상호작용과블레이드요소를통과하는흐름사이에서운동량교환이발생하므로, 로터블레이드단면압력강하의발생원인이되는블레이드요소를통과하는흐름에대해반경방향흐름에대한영향이없고축방향흐름만이작용한다는초기가정은문제가있다. 그러나, 실제로날개끝속도비가 3 이상의범위를가지는경우에서는실제실험데이터와 BEMT에의해계산된성능특성의오차범위가상당히작기때문에초기가정의적용성에대한문제가없다고알려져있으며, 현대식풍력발전용터빈설계및성능평가의표준으로자리매김하고있다. [27] 42
제 3 장수평축풍력터빈블레이드최적설계 3.1 중 대형 (1MW) 로터블레이드최적설계 3.1.1 설계풍속의결정풍력발전용로터블레이드설계는풍력터빈이설치될장소에서다년간의풍향 풍속데이터측정을통한신뢰성있는정격풍속 (rated wind speed, V ) 및설계풍속 (design wind speed, V ) 의결정으로부터 R 시작된다. 정확한설계풍속의결정과신뢰성있는형상설계를위해서는설치입지에서의다년간기상데이터를바탕으로한풍황자원조사가필수적으로선행되어야만하나, 국내실정상원하는입지에서다년간의풍향 풍속데이터를쉽게확보할수없는실정이다. 특히, 해상풍력발전단지의경우로터블레이드의표면도색및코팅작업공정의결정을위해설치위치에서공기의염분함유량에관한데이터까지필요한경우가있으나, 현재까지의국내기상관련측정데이터로부터이러한관련정보를얻기란사실상불가능한일이다. 따라서, 본연구에서는설계자료확보를위한풍력터빈의입지를한국해양대학교내방파제근방으로가정하고정격풍속을결정하기위해설동일 [28] 등에의해 12개월동안측정된한국해양대학교내풍황자원데이터를이용하여풍력터빈로터허브설치높이에서의풍속데이터보정을행하였다. 측정된풍속데이터로부터 Weibull 분포를적용하여얻은지상으로부터 45m 높이에서의정격풍속은 10.13m/s이며, 설계된로터블레이드는로터허브의높이가지상으로부터 60m 정도의위치에설치된다고가정하였다. 이때, 지면으로부터수직방향으로의속도구배 D 43
영향을고려하여허브위치에서의보정풍속을결정하면약 12m/s 이다. 따라서, 형상설계에적용될정격풍속은 12m/s로결정하였다. 설계풍속은일반적으로정격풍속보다 1.3~1.5배낮은범위에서결정되며이경우 10m/s 이다. Table 3.1에설계풍속계산결과를나타내었다. Table 3.1 Decision of the wind velocity 정격풍속 V R = 12m/s 설계풍속 시동풍속 V D = 10m/s V = 3m/s Cut _ in 정지풍속 V = 25m/s Cut _ out 3.1.2 로터블레이드직경및정격회전수의결정로터블레이드직경의결정을위하여식 (3.1) 을이용하였으며, 추정동력계수는 0.45, 동력전달계통장치 (power train) 와발전기 (genera -tor) 의효율은 0.9로가정하였다. D = 8Pr 54.5m ηc ρπv = (3.1) p 3 D 여기서, η ( 동력전달장치및발전기의효율 ) 는 0.9 로가정하였고, C P ( 출력계수 ) 는 0.45, ρ ( 공기밀도 ) 는 1.225kg/m 3 의조건을갖는다. 44
로터블레이드의정격회전수를결정하기위하여식 (3.2) 를적용 하였으며, 설계 TSR 은 7.5 로하였다. N rpm VD = 60 26rpm π D λ = (3.2) 3.1.3 날개끝손실계수의보정회전하는로터블레이드는날개끝와류 (tip vortex) 의발생으로인해순환 (circulation) 분포가감소하게되고, 이로인한날개끝손실 (tip loss) 이발생한다. 이를예측하기위한해석적인방법으로 1919년에 Ludwig Prandtl [29] 이제안한식 (3.3) 과같은날개끝손실예측모델을적용하였다. f ( µ ) 2 2 2 cos 1 [ [(( N / 2)(1 µ ) / µ ) 1 + ( λµ ) /(1 a) e ] = (3.3) π 여기서, µ 는허브로부터팁까지블레이드국부위치를무차원화한변수를나타내며, 0.05(hub) 부터 1(tip) 까지등간격으로 0.05씩증가하도록설정되었다. 따라서, 로터블레이드는허브부터팁까지모두 20 요소로구획된다. 식 (3.3) 에서사용되는축흐름유도계수 a 는초기설계단계부터 BEMT에의한반복계산을통해결정하기가불가능하므로공기역학적항력에의한손실이없고블레이드날개수가무한하다는가정으로부터유도되는운동량이론에서풍력터빈이이론적최대효율을얻을수있는경우의축흐름유도계수 a 를설계에반영한다. 이값은 Betz의이론으로부터유도되는이상적인값과동일하다. 45
블레이드각국부위치에서계산된손실계수값을 Fig. 3.1 에나타 내었다. 46
1.2 Prandtl's Tip Loss Coefficient 1.0 0.8 0.6 0.4 0.2 0.0 0.0 0.2 0.4 0.6 0.8 1.0 r/r Fig. 3.1 Prandtl s tip loss coefficient 47
3.1.4 새로운흐름유도계수들의결정 3.1.3 절에서블레이드국부위치변화에따른날개끝손실을예측하기위해사용된식 (3.3) 의변수 a 는운동량이론으로부터구할수있는이론적최대값이다. 날개끝손실의발생은블레이드국부위치에서의축흐름유도계수변화를초래하므로, 이를반영한축흐름유도계수와회전흐름유도계수를계산하기위하여 BEMT로부터유도된식 (2.51), 식 (2.52) 에서항력항을무시하고식 (3.3) 의날개끝손실계수를적용하면식 (3.4) 와같이표현된다. a a(1 ) 1 1 1 2 f = + 1 +, = 2 2 (3.4) a f f f a 3 3 3 λ µ 식 (3.4) 를이용하여블레이드국부위치변화에따른축흐름유도계수및회전흐름유도계수를다시계산하면, Fig. 3.2(a), Fig. 3.2 (b) 와같으며, 설계 TSR이 3이상인경우식 (3.4) 를이용하여비교적정확한결과를예측할수있다고알려져있다. 48
1.8 Tangential Flow Induction Factor 1.6 1.4 1.2 1.0 0.8 0.6 0.4 0.2 0.0 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 r/r (a) Tangential flow induction factor 0.3336 0.3334 Axial Flow Induction Factor 0.3332 0.3330 0.3328 0.3326 0.3324 0.3322 0.3320 0.3318 0.3316 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 r/r (b) Axial flow induction factor Fig. 3.2 Estimation of flow induction factors as a function of r/r 49
3.1.5 무차원현의길이결정반경방향국부위치에서무차원현의길이를각각결정하기위해서는우선블레이드날개끝부분에사용될익형의선정이필요하다. 그후, 선정된날개끝익형의공력특성데이터를확보하고익형의최대양항비를나타내는받음각에서의양력계수, 축흐름유도계수, 회전흐름유도계수, 블레이드개수, 설계 TSR 등을이용하여식 (3.5) 를통해국부위치에서의무차원현의길이를각각계산할수있다. 2 C 2π 4λµ a = R NλCl 1 a + [ λµ (1 + a )] ( ) 2 2 (3.5) 본연구에서사용된날개끝익형은 NACA 63(2)-415이며최대양항비를나타내는지점에서의양력계수와항력계수값들을 Table 3.2 에나타내었다. 식 (3.5) 에의해계산된현의길이분포는블레이드허브부분으로갈수록급격히증가하는특징을나타내고있으나, 로터블레이드성능을크게좌우하는영역은허브로부터 70% 위치에서블레이드팁까지, 약 30% 정도의영역에지나지않으므로, 이영역에대한정확한형상설계데이터의반영만이루어지면된다. 허브로부터 30% 까지의영역은성능을고려한설계법보다는구조적으로안정된블레이드지지를위한설계방법을고려하여야하는부분이다. 따라서우수한성능을확보하기위한설계적관점을고려하면, 실제로는식 (3.5) 로부터얻은현의길이분포를정확하게반영한블레이드를제작하는것이이상적이다. 50
블레이드제작용이성및제작비용절감을위해성능에큰영향을미치는영역인허브로부터 70% ~ 90% 사이의영역을기준으로 1 차방정식에의한선형근사화를통해전체블레이드국부위치에서의현의길이분포를재구성한다. 선형근사화를통해계산된최종로터블레이드현의길이분포를 Fig. 3.3에나타내었다. Table 3.2 Aerodynamic characteristics - NACA 63(2)-415, X-FOIL RE α C L C D C L/ D 3,000,000 3.5 0.7704 0.00568 135.6338 51
6000 5000 Optimum Chord Length Original Chord Length Chord Length [mm] 4000 3000 2000 1000 0 0.2 0.4 0.6 0.8 1.0 r/r Fig. 3.3 Distributions of chord length at each section 52
3.1.6 입구유입유동각및비틀림각의설정 입구유입유동각 (inflow angle) φ 는식 (3.6) 을통해결정될수 있으며비틀림각 (twist angle) 의계산에사용된다. 1 a tanφ = λµ (1 + a ) (3.6) 비틀림각은앞서계산된 φ 와 α 의관계를통해식 (3.7) 에의해 간단히계산될수있다. θ = φ α (3.7) 계산된비틀림각분포를 Fig. 3.4에나타내었다. 최종블레이드형상은식 (3.8), 식 (3.9) 에의해결정되며, 이수식들은항력에대한영향과날개끝손실의영향을고려한 BEMT 모델로부터유도된수식이다. af σ r σ r 2 1 a = ( C ) 2 x C 2 y 1 a 4sin φ 4sin φ 1 af (3.8) a f σ rcy 1 a = 1+ a 4sinφ cosφ 1 af (3.9) 53
20 18 Twist angles in degrees 16 14 12 10 8 6 4 2 0 0.2 0.4 0.6 0.8 1.0 Local position Fig. 3.4 Twist angle distribution 54
3.1.7 로터블레이드익형선정조건풍력터빈로터블레이드용익형의설계는우수한성능을가지는풍력터빈의연속적인개발을위해필수적인분야이다. Fuglsang과 Madsen [30] 은익형의최적설계관련연구에서적절한공력성능을가진익형은에너지생산단가저감에중요한역할을한다고발표하였다. 오늘날풍력터빈블레이드설계에사용되는익형의종류는, 고전적인비행기날개익형으로사용되는 NACA 시리즈와효율적인에너지추출을위해풍력터빈용으로설계된다수의전문익형시리즈로구분된다. 풍력터빈익형은설계점과탈설계점에서의공력성능, 구조적특징등의부분에서항공기용익형과구분된다. 풍력터빈용익형의개발에대한연구는 1980년대중반부터시작되어현재까지활발히진행되고있으며, Tangler, Somers, Bjork, Timmer, van Rooy, Hill, Garrad, Chaviatopiulos 등과같은연구자들에의해많은연구결과들이발표되었다. 현재까지대다수의풍력터빈용익형은역설계기법을이용하여개발이진행되어왔으며, 익형에관한수많은설계법은 Henne과 Dulikravich [31] 의연구결과에자세히언급되어있다. 원리적으로가장이상적인풍력발전용익형의공력성능은설계되는로터블레이드의운전조건에따라가변적이다. 그러나, 대다수의풍력발전용터빈설계시에일반적으로다음과같은특성이요구된다. 최대출력을발생하기위한조건으로블레이드외부 (outer part) 에적용되는익형의양항비가우수해야한다. 출력제어방식으로피치제어방식이나, 능동실속제어를사용하는경우, 블레이드외부영역에적용되는익형의양항비는설계운전조건에서가장높은값을갖도록설계되어야만한다. 55
실속제어형출력제어방법을적용하는경우양항비는전체터빈의작동범위에걸쳐높은값을유지하도록설계되어야만한다. 이는즉, 받음각의변화가최대양력계수를발생케하는범위이하로결정되어야한다는의미이다. 로터블레이드의내부영역 (in-board part) 에서는양항비가외부영역에비해크게중요하지는않으나블레이드가차지하는면적을줄이기위한목적으로최대양력계수가높은익형의선정이필요하다. 또한, 풍력터빈설계시블레이드를구성하는익형의양력계수가가장높은조건에서운전되도록최적운전설계점의신중한결정이필요하다. 보통광범위한운전영역을가지는풍력터빈로터블레이드는받음각의변화폭이크기때문에, 탈설계점에서우수한공력특성을확보하는것또한중요하다. 만약, 익형의전연으로부터박리가진행되면급격한실속현상이초래되며, 실속유도진동 (stall induced vibration) 등의주된원인이된다. 이에따른위험을최소화하기위한목적으로, 실속제어형풍력터빈의경우, 최대양력이발생하는시점으로부터발생하는흐름의박리현상이후연으로부터서서히진행되도록설계되어야하고, 깊은실속영역에서양력곡선의변화가급격히발생하지않는익형의사용이터빈안정성확보에도움이된다. 자연적인운전조건에서벌레와오염물질들이로터블레이드전연에부착되면, 익형전연부의미소한형상및표면거칠기변화가발생하고, 공력특성감소에의한운전효율이저감되기도한다. 블레이드전연의표면거칠기는층류에서난류로확장되는천이영역의발생시점을앞당기게되고경계층두께의급격한확장을초래한다. 최대양력의감소와익형표면마찰계수의증가는, 결론적으로에너지생산량의감소로이어지기때문에, 특히실속형출력제어방식을채용한터빈의경우표면마찰계수의변화에따른최대양력계수의변화폭이크지않은익형을주의깊게선정할필요가있다. 56
풍력터빈블레이드의허브근방영역에배치되는익형은제한된블레이드재질의무게와팁에서발생하는처짐및비틀림변위를충분히견딜수있을정도의강성을확보할수있어야한다. 이러한강성을확보하기위해보통양항비와같은공력특성을희생하고익형의최대두께를증가시킬수있는두꺼운익형시리즈를채택한다. 일반적으로풍력발전용터빈블레이드를구성하는익형은높은양항비를가지는두꺼운형상, 전연부의표면거칠기에상관없이높은최대양력계수를유지할수있을것, 모든운전영역에걸쳐높은양항비를유지할수있을것, 실속후급격한공력특성의저하현상을보이지않을것등의요구조건을만족해야한다. 하지만, 이러한모든조건을충족시키는익형을설계하기란여간까다로운일이아닐수없다. 57
3.1.8 2차원익형공력특성의예측정확한로터블레이드형상설계및성능예측을위해서는로터블레이드를구성하는익형의신뢰성있는공력특성데이터의확보가필수적이나, 수많은종류의익형에대한신뢰성있는실험공력특성데이터의확보는현실적으로거의불가능한실정이다. 신속한피드백 (feed back) 에의한효율적인설계의진행을위해설계자는다양한익형의공력특성데이터를쉽게확보할수있어야한다. 풍력선진국가에서는일반적으로 X-Foil(MIT Aero&Astro) [32] 등과같은수치해석소프트웨어를이용하여 2차원익형공력특성데이터를확보하고있으며, 최근들어전산해석기법의발달에힘입어 CFD 기법을활용한사례들도활발히보고되고있다. X-Foil은 2차원 panel method를기본으로하고, CFD에비해상당히신속한해석이가능하며, 실속영역이전의받음각변화에대해서비교적신뢰성있는공력특성데이터의확보에상당한도움이된다. 본연구에사용된익형의종류는날개끝단으로부터 NACA 63(2)- 415, NACA 63(2)-418, DU 93W2-210, DU 91W2-250, FFAW-301 이다. 보통풍력발전용익형은날개끝부분에서 30% 위치까지 NACA 5, 6자리계열익형을사용하고, 그이하의영역에서는구조적안정성, 높은받음각, 익형표면거칠기등의변화를적절히고려하여공력특성의변화가급격히발생하지않는두꺼운익형시리즈를주로사용한다. [33] 일반적으로팁에서허브부분으로갈수록익형두께가두꺼워지는형태의프로파일을공력특성데이터의분석을통해적절히선정하는것이중요하다. 본설계에사용된익형의받음각변화에따른양력계수그래프를 Fig. 3.5에나타내었다. 58
Lift Coefficient 1.8 1.7 1.6 1.5 1.4 1.3 1.2 1.1 1.0 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 NACA 63(2)-415 (X-foil, Re=3,000,000) NACA 63(3)-418 (X-foil, Re=3,000,000) DU-93-W2-210 (X-foil, Re=3,000,000) DU-91-W2-250 (X-foil, Re=3,000,000) FFA-W-301 (X-foil, Re=3,000,000) 0 1 2 3 4 5 6 7 8 9 10111213141516171819202122232425 AOA Fig. 3.5 Estimation of the lift coefficients for several airfoils by X-Foil 59
Fig. 3.5에나타낸양력계수분포는모두 X-Foil에의해계산된데이터이며, 실속이후특정받음각범위에서실험데이터와비교할때다소오차가발생한다. 일반적으로실험데이터와수치해석데이터의오차범위는익형의표면경계층박리지점및천이영역의부정확한예측에의한수치해석기법상의특징에의해발생되며, 이로인해수치해석데이터는실속현상이지연되는특징을나타냄에따라, 실속영역이후의공력특성데이터들은실험데이터와비교했을때비교적정확하지못한값을예측하는것이일반적이다. Fig. 3.6에수치해석결과와실험결과를비교해나타내었다. Fig. 3.6에서두결과가실속전해당받음각영역내에서상당히잘일치하고있음을알수있으며, 로터블레이드설계및성능해석과정에서 X-Foil에의한수치해석데이터를사용하여도충분한신뢰성을확보할수있을것이라사료되며, 이를통해신속하게다양한익형의공력특성데이터를확보할수있으고설계자는상당히효율적으로설계를진행할수있다. 보통익형의양력계수값이상승한다고하더라도양력계수의증가분보다항력계수의증가분이더크다면양력과항력의비로나타내어지는익형의양항비는양력계수가상승함에도불구하고오히려감소하는경우가발생한다. 따라서정확한익형의공력성능을예측하기위해서는블레이드국부위치에서의익형공력특성을양항비로나타낼필요성이있으며, 이때블레이드성능에큰영향을미치는허브로부터 70% ~ 100% 까지범위내의양항비는반드시 100 이상에서결정될수있도록하는점이우수한출력성능확보를위해중요하다. 60
Lift Coefficient 1.8 1.7 1.6 1.5 1.4 1.3 1.2 1.1 1.0 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 Good agreement with Exp. NACA 63(3)-418(X-foil, Re=3,000,000) NACA 63(3)-418 (Experiment, Re=3,000,000) NACA 63(2)-415(Experiment, Re=3,000,000) NACA 63(3)-415(X-foil, Re=3,000,000) 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 AOA Fig. 3.6 Comparison of lift coefficients calculated by X-FOIL with experiments 61
Fig. 3.7에익형의양항비를그래프로나타내었다. 각각의블레이드국부위치에서의양항비는 FFA 익형을제외하고설정받음각의범위내에서모두 100을넘는값을가진다. FFA 익형의양항비가 100을넘지못하는이유는블레이드의공력성능에기여하는영역에적용될익형이아니라, 블레이드와허브의구조적인연결부분의지지익형으로적용되기때문에양항비가다른익형들에비해다소낮은경향을나타낸다. 하지만, 지지익형인 FFA의양항비가거의 100에근접하며, 날개끝 70% ~ 100% 범위에서익형의양항비는거의 140에육박하는수치를나타내고있다. 따라서, 이러한조합의익형을사용하면전체적으로상당히우수한성능을가지는이상적인블레이드설계가가능할것이라사료된다. 상기기술한 BEMT 모델의적용에의해최종적으로설계된 3차원로터블레이드형상을바탕으로가시적이해를돕기위해너셀, 타워, 스피너 (spinner), 허브등의구조적요소를포함하는 1MW 풍력터빈전체시스템 3차원모델링을수행하였으며, 본풍력터빈시스템을 FIL-1000로명명하였다. Table 3.3에 FIL-1000의자세한설계제원을요약정리하였고, Fig. 3.8에 FIL-1000의 3차원시스템설계결과를나타내었다. 62
Lift to Drag Ratio 160 150 140 130 120 110 100 90 80 70 60 50 40 30 20 10 0-1 0 1 2 3 4 5 6 7 8 9 101112131415161718192021222324 AOA NACA 63415 NACA 63418 DU93W2-210 DU91W2-250 FFA W-301 Fig. 3.7 Comparison of lift to drag ratios, calculated by X-Foil 63
Table 3.3 Summarization of 1MW HAWT(FIL-1000) design parameters 64
(a) Isometric view of 3D rotor model (b) Full 3D FIL-1000 wind turbine system modeling Fig. 3.8 3D configurations of FIL-1000 HAWT 65
3.2 중형로터블레이드 (100kW) 최적설계 3.2.1 로터블레이드직경및정격회전수의결정본연구에서는 3.1절에언급한 MW 급대형터빈설계법을근거로 100kW 풍력터빈로터블레이드설계를수행하였다. 로터블레이드직경은식 (3.1) 에의해결정되었고, 추정동력계수 0.45, 동력전달계통장치와발전기효율을 0.9로가정하였다. 공기밀 도는 1.225 / 3 kg m 이다. 설계풍속은 8m/s이며, 이때설계 TSR은 6으로 하였다. 로터블레이드정격회전속도는식 (3.2) 에의해결정되었다. 자세한설계초기조건을 Table 3.4 에요약정리하였다. Table 3.4 General design parameters (100kW) Rated power(kw) 100 Rated wind velocity(m/s) 10 Design TSR 6 Rotor blade diameter(m) 22.65 Rated rotational speed(rpm) 40.46 3 Air density( kg / m ) 1.225 3.2.2 날개끝손실계수의보정날개끝순환분포의감소에의해발생하는날개끝손실의예측을위해식 (3.3) 을적용하였다. 블레이드국부위치변화에따른손실계수변화값을 Fig. 3.9에나타내었다. 66
1.2 Prandtl's Tip Loss Coefficient 1.0 0.8 0.6 0.4 0.2 0.0 0.0 0.2 0.4 0.6 0.8 1.0 r/r Fig. 3.9 Estimation of Prandtl s tip loss coefficient as a function of r/r 67
3.2.3 새로운흐름유도계수들의결정초기설계단계에적용될흐름유도계수들의결정을위해 BEMT로부터유도된식의항력항을무시하고정리한식 (3.4) 를사용하였다. Fig. 3.10에블레이드국부위치변화에따른흐름유도계수변화를그래프로나타내었다. 0.34 0.33 Axial Flow Induction Factor 0.32 0.31 0.30 0.29 0.28 0.27 0.26 0.25 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 r/r (a) Axial flow induction factor Fig. 3.10 Estimation of flow induction factors as a function of r/r (continued) 68
3.0 Tangential Flow Induction Factor 2.5 2.0 1.5 1.0 0.5 0.0 0.0 0.2 0.4 0.6 0.8 1.0 r/r (b) Tangential flow induction factor Fig. 3.10 Estimation of flow induction factors as a function of r/r 69
3.2.4 무차원현의길이결정무차원현의길이를결정하기위해식 (3.5) 를적용하였다. Table 3.5에블레이드설계시적용될익형의공력특성설계기준값을나타내었으며, Fig. 3.11에결정된블레이드현의길이분포를나타내었다. 설계에적용된익형은 NACA63(2)-415V이다. Table 3.5 NACA 63(2)-415V, Wind tunnel test RE α C L C D C L/ D 1,600,000 4.3949 0.7988 0.011 72.6182 3500 3000 Original Chord Length Optimum Chord Length Chord Length [mm] 2500 2000 1500 1000 500 0 0.2 0.4 0.6 0.8 1.0 r/r Fig. 3.11 Estimation of chord length 70
3.2.5 비틀림각의결정 비틀림각의계산을위해식 (3.6), 식 (3.7) 을이용하였으며, Fig. 3.12 에비틀림각분포를나타내었다. 24 22 20 Twist Angles in Degrees 18 16 14 12 10 8 6 4 2 0 0.2 0.4 0.6 0.8 1.0 r/r Fig. 3.12 Estimation of twist angle variation 71
3.2.6 2차원익형공력특성본설계에서반영한익형은 NACA 63 (2) 415V 이며, Fig. 3.13, Fig. 3.14, Fig. 3.15에양력계수, 항력계수, 양항비에관한실험데이터를나타내었다. 최종적인설계변수들을 Table 3.6에나타내었으며, 블레이드 3차원형상모델링결과를 Fig. 3.16, Fig. 3.17에나타내었다. Lift Coefficient 1.5 1.4 1.3 1.2 1.1 1.0 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0.0-0.1-0.2-0.3-0.4-0.5-0.6-0.7 NACA 632415V, Velux wind tunnel test RE=1,600,000-10 -8-6 -4-2 0 2 4 6 8 10 12 14 16 18 20 22 24 26 28 30 32 34 Angle of Attack Fig. 3.13 Lift coefficient of NACA 63(2)-415V 72
0.45 0.40 NACA 632415V, Velux wind tunnel test RE=1,600,000 0.35 Drag Coefficient 0.30 0.25 0.20 0.15 0.10 0.05 0.00-12-10-8 -6-4 -2 0 2 4 6 8 10 12 14 16 18 20 22 24 26 28 30 32 34 AOA Fig. 3.14 Drag coefficient of NACA 63(2)-415V 73
80 70 60 50 Lift to Drag Ratio 40 30 20 10 0-10 -20-30 -40 NACA 632415V, Velux wind tunnel test RE=1,600,000-50 -12-10 -8-6 -4-2 0 2 4 6 8 10 12 14 16 18 20 22 24 26 28 30 32 34 AOA Fig. 3.15 Lift to drag ratio of NACA 63(2)-415V 74
Table 3.6 Summarization of 100kW HAWT(FIL-100) design parameters 75
Fig. 3.16 Isometric view of 3D rotor model (FIL-100) 76
Fig. 3.17 Isometric view of 3D rotor model (FIL-100) 77
3.3 소형로터블레이드 (20kW) 최적설계 3.3.1 설계변수들의결정본연구에서는 3.1절과 3.2절에언급한 BEMT에의한풍력터빈설계법을근거로 20kW 풍력터빈로터블레이드설계를수행하였다. 관련설계제반수식에의해결정된로터블레이드직경은 10.13m 이며, 공기밀도는 1.225 / 3 kg m 이다. 정격풍속은 10m/s로결정하였으며 설계풍속은 8m/s이다. 설계 TSR은 6이며, 로터블레이드정격회전속도는 90.49rpm으로결정하였다. 이미, 3.1절, 3.2절에서대형및중형터빈의설계방법에대해상세히논하였으므로소형터빈의자세한설계과정은생략한다. 이상의초기설계사항을요약하면 Table 3.7과같다. 날개끝손실계수예측결과를 Fig. 3.18에나타내었으며, 관련흐름유도계수들의예측결과를 Fig. 3.19에나타내었다. Table 3.7 General design parameters (20kW) Rated power(kw) 20 Rated wind velocity(m/s) 10 Design TSR 6 Rotor blade diameter(m) 10.13 Rated rotational speed(rpm) 90.49 3 Air density( kg / m ) 1.225 78
1.2 Prandtl's Tip Loss Coefficient 1.0 0.8 0.6 0.4 0.2 0.0 0.0 0.2 0.4 0.6 0.8 1.0 r/r Fig. 3.18 Prandtl s tip loss coefficient 79
0.34 0.33 Axial Flow Induction Factor 0.32 0.31 0.30 0.29 0.28 0.27 0.26 0.25 0.0 0.2 0.4 0.6 0.8 1.0 r/r 3.0 (a) Axial flow induction factor Tangential Flow Induction Factor 2.5 2.0 1.5 1.0 0.5 0.0 0.0 0.2 0.4 0.6 0.8 1.0 r/r (b) Tangential flow induction factor Fig. 3.19 Flow induction factors 80
무차원현의길이를결정하였으며, 20kW 로터블레이드형상화를위해 NREL의 S809 시리즈를사용하였다. Table 7에익형의공력특성설계기준값을나타내었으며, Fig. 3.20에블레이드스팬방향변화에따른현의길이변화를나타내었다. Table 3.8 Aerodynamic characteristics S-809, Wind tunnel test (DUT) RE α C L C D C L/ D 1,000,000 6.16 0.851 0.0095 87.58 1800 1600 1400 Original Chord Length Optimum Chord Length Chord Length [mm] 1200 1000 800 600 400 200 0 0.0 0.2 0.4 0.6 0.8 1.0 r/r Fig. 3.20 Comparison of chord length distributions 81
블레이드비틀림각분포계산결과를 Fig. 3.21에나타내었으며, 본설계에사용된익형의공력특성변화그래프를 Fig. 3.22, Fig. 3.23. Fig. 3.24에각각나타내었다. 익형의공력특성실험데이터는 1997년 Sommers에의해 Delft University of Technology의풍동에서레이놀즈수 10 6 의조건에서수행된결과이다. Table 3.9에 20kW 로터블레이드설계결과를요약하여나타내었으며, Fig. 3.25에로터블레이드 3차원형상설계결과를나타내었다. 22 20 18 Twist Angles in Degrees 16 14 12 10 8 6 4 2 0-2 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 r/r Fig. 3.21 Twist angle distributions 82
1.2 1.1 1.0 0.9 Lift Coefficient 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0.0-0.1-3 -2-1 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 AOA S809, DUT Wind tunnel test RE=1,000,000 Fig. 3.22 Lift coefficient of S809 83
Drag Coefficient 0.20 0.19 0.18 0.17 0.16 0.15 0.14 0.13 0.12 0.11 0.10 0.09 0.08 0.07 0.06 0.05 0.04 0.03 0.02 0.01 0.00 S809, DUT Wind tunnel test RE=1,000,000-3 -2-1 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 AOA Fig. 3.23 Drag coefficient of S809 84
100 90 80 S809, DUT Wind tunnel test RE=1,000,000 Lift to Drag Ratio 70 60 50 40 30 20 10 0-3 -2-1 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 AOA Fig. 3.24 Lift to drag ratio of S809 85
Table 3.9 Summarization of 100kW HAWT(FIL-100) design parameters ROOT TIP Local position Twist angle(deg) Position(mm) Chord length(mm) Airfoil series 0.20 19.54 1014.00 689.23 S809 0.25 15.86 1267.50 662.54 S809 0.30 12.95 1521.00 635.86 S809 0.35 10.65 1774.50 609.17 S809 0.40 8.81 2028.00 582.48 S809 0.45 7.31 2281.50 555.80 S809 0.50 6.07 2535.00 529.11 S809 0.55 5.03 2788.50 502.43 S809 0.60 4.15 3042.00 475.74 S809 0.65 3.39 3295.50 449.05 S809 0.70 2.72 3549.00 422.37 S809 0.75 2.13 3802.50 395.68 S809 0.80 1.59 4056.00 369.00 S809 0.85 1.06 4309.50 342.31 S809 0.90 0.52 4563.00 315.62 S809 0.95-0.12 4816.50 288.94 S809 1.00-0.12 5070.00 262.25 S809 86
Fig. 3.25 Isometric view of 3D rotor model (FIL-20) 87
제 4 장성능해석소프트웨어개발및평가 4.1 소프트웨어개발 국내의로터블레이드설계관련기술에대한명확한기준이없으며, 현재국내에서는로터블레이드를설계하는데있어영국의 Garrad Hassen 사에서개발한 BLADED FOR WINDOWS [34] 등과같은국외상용소프트웨어에의존한설계및성능평가가이루어지고있는실정이다. 또한, 현재까지풍력발전시스템요소기술에대한명확한설계기준및성능평가기준이마련되어있지않은국내의현실을감안할때, 국내에서도블레이드설계및성능평가기술에대한명확한기준등이제시될필요가있다. 따라서, 본연구에서는다양한손실모델을적용한 BEMT 이론을활용하여로터블레이드설계를수행하고, 성능평가를위해국산성능해석소프트웨어의개발을진행하였다. 개발된성능평가소프트웨어는국내해상풍력발전활성화를의미하는뜻에서 POSEIDON 으로명명하였다. POSEIDON 의성능평가알고리즘은 Fig. 4.1 과같은절차로구성된다. 제 2 장에서설명한 BEMT 모델의성능해석과정을반영하여개발된소프트웨어의화면구성을 Fig. 4.2, Fig. 4.3, Fig. 4.4, Fig. 4.5 에나타내었다. Fig. 4.2는성능해석을위한기본설계정보를입력하는화면이며, 사용자는로터블레이드초기형상설계정보를확보해야한다. 사용자로부터초기입력변수로써블레이드형상정보를입력받은후날개끝속도비, 블레이드반경, 블레이드개수, 밀도, 유입속도등의제반정보를입력받는다. 88
기본적인설계변수의입력이끝나면, 블레이드각각의단면을정의하는익형의공력특성데이터를입력하며, Fig. 4.3에익형공력특성입력화면을나타내었다. 익형공력특성데이터입력화면은모두 15 단면에서각각의공력특성데이터를입력받을수있도록구성하였으며, 이는블레이드전체공력특성을예측함에있어적분오차를줄이고자하는목적이다. 익형공력특성데이터의입력후메인페이지에서실행버튼을클릭하면 Fig. 4.1과같은순서도를통해관련된모든정보에관한계산을수행하며, 계산된정보는 Fig. 4.4의화면과같이나타난다. Fig. 4.5는입력받은공력특성데이터를그래프로나타낸화면이며, 실제성능계산시참조되어사용된공력특성범위를표시하여보여준다. 이는특정 TSR 영역에대한출력성능을계산한후블레이드실속발생여부등의판단기준으로활용될수있다. 89
Fig. 4.1 Flow chart of POSEIDON 90
Fig. 4.2 Main input screen design 91
Fig. 4.3 Aerodynamic characteristics loading 92
Fig. 4.4 Result output 93
Fig. 4.5 Aerodynamic characteristics plotting screen 94
4.2 실속후공력특성의예측 풍력발전용로터블레이드는낮은범위의 TSR에서운전되는경우에유입풍속의급격한증가에따라블레이드특정영역, 혹은전체영역에대한실속이발생된다. 항공기의경우깊은실속이발생한상태에서는운항이불가능하지만, 풍력터빈의경우출력제어등의목적으로실속상태에서운전되기도하므로, BEMT에의한정확한성능예측을위해익형의실속후공력특성데이터의확보는필수적인요소이다. 사실상다양한익형의실속후공력특성을포함하는실험데이터의확보는설계자들에게가장해결하기힘든부분이다. 대부분의경우공력특성데이터의확보를위해 X-FOIL등의익형공력특성예측소프트웨어를활용하고있으나, 이는실속후공력특성의보정을위해 Viterna-Corrigan [35] 등이제안한모델에의한수치해석데이터의재가공이필요하다. BEMT에대한적용을위해 Viterna는실속후공력특성예측모델을적용하였으며, MOD-0, MOD-1 터빈의성능해석을위해풍동실험으로부터획득한실험데이터와실속후모델에의해예측된후보정데이터를적용한연구를수행하였다. 95
Cd, Max = 1.11+ 0.018AR ( α = 90 deg) (4.1) C = B sin α + B cosα ( α = 15deg to 90deg) (4.2) d 2 1 2 여기서, B1 = C d, Max (4.3) 2 Cd, stall Cd, Max sin αstall B2 = (4.4) cosα stall 2 cos α C1 = A1 sin 2α + A2 ( α = 15deg to 90deg) (4.5) sinα 여기서, B 2 1 A 1 = (4.6) sinαstall A2 = ( Cl, stall Cd, Max sinαstall cos αstall ) (4.7) cos 2 α stall 96
4.3 FIL-1000 성능평가 (1MW turbine) 본연구를통해설계된 1MW 풍력발전용로터블레이드 (FIL- 1000) 의성능평가를위해 4.1절에언급한 POSEIDON을사용하였다. FIL-1000 로터블레이드를구성하는익형은허브로부터팁부분까지 FFA-W-301, DU-91W2-250, DU-93W2-210, NACA63(2)-418, NACA63(2)-415 이다. 익형공력특성은 X-FOIL로부터획득하였으며, 실속후영역에대한공력특성보정은 Viterna-Corrigan의실속후예측모델을사용하여받음각범위 60 까지양력계수, 항력계수에대한데이터를 0.5 간격으로획득하여초기입력데이터로사용하였다. Fig. 4.6(a) 에블레이드국부위치에서의받음각변화를 TSR 변화에따라나타내었다. FIL-1000에사용된익형들의실속각은약 13 부근에서형성되고있으며, TSR=5, 6, 7, 8, 9의경우실속각범위내에서거의대부분의받음각이결정되고있으므로, 로터블레이드출력은설계시예측된적정성능을유지할것이라고예상할수있으나, TSR=2, 3, 4의경우받음각이실속각을상당히벗어난범위에서결정되고있으므로, 이영역에대한시스템출력특성은다소낮게나타날것임을예측할수있다. Fig. 4.6(b) 에는블레이드국부위치에서의 Prandtl의날개끝손실변화를 TSR에따라비교한그래프를나타내었다. 대부분의 TSR 영역에대해날개끝손실은허브로부터날개팁까지진행될수록점점크게나타나고있음을알수있으며, TSR=9인경우에날개끝손실이가장낮게나타나고있음을알수있다. 날개끝손실은시스템의전체출력에영향을미치는변수이며, Fig. 4.6(a) 의결과와마찬가지로 TSR=5 이하인영역에서날개 97
끝손실이블레이드전영역에걸쳐발생하고있음을알수있다. FIL-1000은설계 TSR=7.5의경우에서최적의공력특성을확보할수있는받음각의영역이결정되고있으며, 날개끝손실도비교적낮은범위에서형성되고있다. Fig. 4.7(a), Fig. 4.7(b) 에축흐름유도계수와회전흐름유도계수변화를 TSR에따라각각나타내었다. 항력의효과를무시한운동량방정식으로부터유도되는가장이상적인조건에서의축흐름유도계수의값은 0.333 이다. 따라서, 최종적인성능해석을통한축흐름유도계수결과값이이상적인조건범위에근접할수록설계된로터블레이드의출력특성은최대출력치에근접하게된다. Fig. 4.7(a) 의결과에서, TSR=2, 3, 4, 5의경우이상적인축흐름유도계수값에대해상당히낮은분포의값을가지는특징을나타내고있으며, 이범위에서는높은출력특성을기대하기어렵다고판단되며, TSR=6, 7, 8의경우대체적으로결과값의범위가이상적인조건부근에형성되고있음을알수있으므로, 이범위의 TSR에서우수한출력특성을나타낼것임을알수있다. TSR=9의경우다소높은범위의결과분포를나타내고있으므로, TSR=6, 7, 8의결과보다다소낮은출력특성을나타낼것이라예상되나, TSR=5 이하의범위보다는상대적으로우수한출력특성을나타낼것임을알수있다. Fig. 4.8에 FIL-1000의성능곡선을 TSR의변화에따라그래프로나타내었다. FIL-1000은 TSR=6, 7에서가장좋은성능특성을나타내고있음을알수있으며, TSR=6 이전영역은 TSR의정의에따라블레이드상류로부터접근하는유속의범위가정격출력을발생하는풍속범위보다높은값을가지므로, 블레이드전연으로유입되는상대유입 98
각도가증가하고, 동시에실속각이상으로받음각이증가하면서급격한실속에의한공력특성저하가초래된다. 따라서, 이범위의 TSR에서상당히낮은출력특성을나타내고있음을알수있다. 또한, TSR=7 이후는입구유속이정격출력을내는범위보다감소함에따라상대적으로블레이드전체영역에걸쳐낮은범위의받음각이형성되고받음각의범위가낮아짐에따라적정받음각의분포에비해익형공력특성이저하되므로출력특성은점점감소하는경향을보인다. FIL-1000 로터블레이드는 TSR=7에서 0.47 정도의출력계수를나타내고있음을알수있으며, 입구풍속이 12m/s이고, 표준대기압공기밀도조건에서약 1,133.8kW의정격출력을발생시킬수있다. 99
(a) Comparison of Angle of Attack according to TSR (b) Comparison of Prandtl s tip loss according to TSR Fig. 4.6 AOA and tip loss distributions for a range of TSR 100
(a) Comparison of axial flow induction factors (b) Comparison of tangential flow induction factors Fig. 4.7 Distribution of the flow induction factors for a range of TSR 101
Fig. 4.8 Power coefficient TSR performance curve 102
4.4 FIL-100 성능평가 (100kW turbine) FIL-100은 100kW 수평축풍력터빈이며, 블레이드를구성하는익형으로써 NACA 63 (2) -415V를사용하였다. 본연구에서는실험에의해측정된공력특성값과수치해석에의해예측된공력특성값의적용에대한성능계수값의변화를분석하고자 Velux wind tunnel test [36] 를통해얻은실험공력특성과 X- FOIL에의해예측된공력특성두가지를비교대상으로하였다. NACA 63 (2) -415V의양력계수변화에대한비교결과를 Fig. 4.9에나타내었다. 제공받은풍동실험결과의경우, 실속영역이전의경우에해당하는공력특성데이터만을포함하고있어, 실속후공력특성보정식을적용하였다. 수치해석데이터의경우, 실속영역이전공력특성데이터는 X-FOIL에의해예측되었으며, 실속후공력특성데이터는실험데이터와마찬가지로 Viterna-Corrigan의보정식을사용하였다. 실속이전영역에대한공력특성데이터의경우 X-FOIL에의해예측된데이터와실험에의해예측된데이터의결과가정량적으로매우잘일치하고있음을알수있으며, 실속이발생한후부터 X- FOIL에의해예측된데이터가실험데이터보다다소높은양력계수값을나타내고있음을알수있다. Fig. 4.10에항력계수변화에대한비교결과를나타내었다. 항력계수도마찬가지로, 실속이전의영역에대해실험결과와수치해석결과가정량적으로잘일치하고있음을알수있으며, 실속이발생한후부터미소한값의차이가나타나고있음을알수있다. 103
Lift Coefficient 1.80 1.65 1.50 1.35 1.20 1.05 0.90 0.75 0.60 0.45 0.30 0.15 0.00-0.15-0.30-0.45-0.60-0.75-0.90 Measured, Velux wind tunnel Predicted, X-FOIL RE=1,600,000-16-12-8 -4 0 4 8 12 16 20 24 28 32 36 40 44 48 52 56 60 64 68 72 76 80 84 88 92 96 Angle of Attack Fig. 4.9 Comparison of lift coefficient (NACA63 (2) -415) 104
1.2 1.1 1.0 0.9 Drag Coefficient 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0.0 Measured, Velux wind tunnel Predicted, X-FOIL RE=1,600,000-16-12-8 -4 0 4 8 12 16 20 24 28 32 36 40 44 48 52 56 60 64 68 72 76 80 84 88 92 96 Angle of Attack Fig. 4.10 Comparison of drag coefficient (NACA63 (2) -415) 105
Fig. 4.11(a), Fig. 4.11(b) 에실험데이터를사용한경우의블레이드반경방향위치변화에따른받음각의변화를나타내었다. 실험공력특성의경우실속이발생하는받음각은 15.5455 이며, TSR 4.5부터 TSR 10까지의경우실속이전영역에서받음각이결정되고있음을알수있다. Fig. 4.11(b) 는수치해석에의해예측된공력특성데이터를사용한경우를나타내며, 받음각의형성범위가 Fig. 4.11(a) 의경우와유사한특징을나타내고있으나, TSR 4의경우까지비실속영역에포함되고있음을알수있다. 대체로수치해석에의해예측된경우의받음각이미소하게나마최적공력성능을나타내는범위에근접하고있다. 이러한차이는최종적으로예측된터빈의성능특성에영향을미칠것이라사료되며, 그범위는추후성능비교그래프를통해논의하고자한다. Fig. 4.12(a), Fig. 4.12(b) 에실험과수치해석에의해예측된블레이드반경방향위치변화에따른축흐름유도계수의변화를나타내고있다. 실험결과에의해예측된축흐름유도계수의경우최적성능을나타내는설계 TSR 영역부근의 TSR=5, 5.5, 6에서이상적인조건에비해다소낮은영역에걸쳐형성되고있으며, 수치해석에의해예측된공력특성을사용한경우, 이상적인조건범위에더욱근접하여그래프가형성되고있음을알수있다. 이러한차이는익형공력특성데이터의차이에의해발생하는것이라사료되며, 결과적으로풍력터빈의성능변화에영향을미칠것이라예상된다. Fig. 4.13(a), Fig. 4.13(b) 에실험과수치해석에의해예측된블레이드반경방향위치변화에따른회전흐름유도계수의변화를나타내고있다. 실험과수치해석데이터를사용한결과가대체로잘일치하고있으나, 공력특성데이터의차이로인해발생하는미소한값의차이가나타나고있다. 106
최종적으로 FIL-100의출력성능특성그래프를 Fig. 4.11에나타내었다. 실선은실험공력특성데이터를사용해예측된결과를나타내고, 파선은수치해석공력특성데이터를사용해예측된결과를나타낸다. 전체적으로수치해석에의해예측된풍력터빈의성능이실험공력특성을사용한경우보다다소높게나타나고있다. 이는수치해석에의한공력특성데이터값이실험값에비해다소높게예측되었기때문에발생하는문제이다. 따라서, BEMT 해석에있어정확한익형공력특성데이터의사용은시스템의출력특성을예측함에있어서중요한변수가됨을알수있다. 107
Axial Flow Induction Factor 0.70 0.65 0.60 0.55 0.50 0.45 0.40 0.35 0.30 0.25 0.20 0.15 0.10 0.05 0.00 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 1.1 1.2 1.3 r/r TSR=2.5 TSR=3.0 TSR=3.5 TSR=4.0 TSR=4.5 TSR=5.0 TSR=5.5 TSR=6.0 TSR=7.0 TSR=8.0 TSR=9.0 TSR=10.0 (a) Comparison of axial flow induction factors (by measured lift&drag) Axial Flow Induction Factor 0.70 0.65 0.60 0.55 0.50 0.45 0.40 0.35 0.30 0.25 0.20 0.15 0.10 0.05 0.00 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 1.1 1.2 1.3 r/r TSR=2.5 TSR=3.0 TSR=3.5 TSR=4.0 TSR=4.5 TSR=5.0 TSR=5.5 TSR=6.0 TSR=7.0 TSR=8.0 TSR=9.0 TSR=10.0 (b) Comparison of axial flow induction factors (by predicted lift&drag) Fig. 4.11 Distributions of axial flow induction factors 108
0.16 Tangential Flow Induction Factor 0.14 0.12 0.10 0.08 0.06 0.04 0.02 0.00 TSR=2.5 TSR=3.0 TSR=3.5 TSR=4.0 TSR=4.5 TSR=5.0 TSR=5.5 TSR=6.0 TSR=7.0 TSR=8.0 TSR=9.0 TSR=10.0 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 1.1 1.2 1.3 r/r 0.14 (a) Comparison of tangential flow induction factors (by measured lift&drag) Tangential Flow Induction Factor 0.12 0.10 0.08 0.06 0.04 0.02 0.00 TSR=2.5 TSR=3.0 TSR=3.5 TSR=4.0 TSR=4.5 TSR=5.0 TSR=5.5 TSR=6.0 TSR=7.0 TSR=8.0 TSR=9.0 TSR=10.0 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 1.1 1.2 1.3 r/r (b) Comparison of tangential flow induction factors (by predicted lift&drag) Fig. 4.12 Distributions of tangential flow induction factors 109
0.50 0.45 0.40 Power Coefficient 0.35 0.30 0.25 0.20 0.15 0.10 0.05 CP(BEMT), Measured airfoil data CP(BEMT), Predicted airfoil data by X-FOIL 1 2 3 4 5 6 7 8 9 10 11 TSR Fig. 4.13 Comparison of the power coefficient on the FIL-100 110
4.5 FIL-20 성능평가 (20kW turbine) FIL-20은 20kW 수평축풍력터빈이며, 블레이드를구성하는익형으로써 NREL의 S-809를사용하였다. 로터블레이드성능평가를위해 Delft university of technology의풍동에서 Sommers [37] 에의해수행된공력특성결과를사용하였으며, X-FOIL에의해예측된수치해석데이터를비교목적으로사용하였다. 실험및수치해석의레이놀즈수는모든경우에서 10 6 이다. 수치해석공력특성의경우실속후공력특성예측을위해 Viterna- Corrigan 모델을사용하였으나, 실험데이터에대한보정은행하지않았다. S809 익형은 21% 의두께비를가지며, 높은최대양력계수를가지도록설계되었으며, 공력특성변화가표면거칠기변화에대해상대적으로덜민감한형태의수평축풍력발전용익형이다. Fig. 4.14, Fig. 4.15에 X-FOIL에의해예측된익형의공력특성변화를나타내었다. Fig. 16과 Fig. 4.17에실험공력특성과수치해석공력특성의비교결과를나타내었으며, 실속이전영역에대한비교를수행하였다. Fig. 4.16에나타낸양력계수의경우실험과수치해석데이터모두정량적으로잘일치하고있음을알수있으나, 받음각의범위 7 부터 9.5 까지의결과는서로미소한차이를나타내고있다. Fig. 4.17의항력계수비교데이터의경우전체적으로서로잘일치하고있으나, 마찬가지로받음각 7 이상의범위에서부터미소한차이를나타내고있다. 상대적으로항력계수의오차범위가그래프상에서크게나타나는것은항력계수의값이 10-3 승오더에서결정되기때문이다. Fig. 4.18(a), Fig. 4.18(b) 에실험공력특성및수치해석데이터를사 111
용해예측된 FIL-20의축흐름유도계수변화를나타내었다. Fig. 4.19(a), Fig. 4.19(b) 에회전흐름유도계수의변화를나타내었다. 두경우모두실험결과와수치해석결과가비슷한분포를나타내고있으나국부적으로미소한차이를나타내고있다. 이는실험과수치해석사이에발생하는공력특성데이터의오차에기인한다. Fig. 4.20에 FIL-20의성능계수비교결과를나타내었다. 실험데이터의경우실속후보정을행하지않았으므로, 수치해석데이터와 TSR 6, 7, 8, 9의경우에한해비교되었다. 비교결과 TSR 6, 9의경우실험결과와수치해석결과가잘일치하고있음을알수있었으나, TSR=8, 9의경우미소한차이가나타나고있다. 이는 Fig. 4.16의양력계수비교결과에서받음각의범위 7 부터 10 까지의영역에서공력특성데이터가서로미소한차이를나타내는것과일치한다. 실제로, TSR=8, 9의경우 BEMT에의한해석과정에서참조값으로사용되는공력특성데이터의범위가약 6.5 부터 10 사이에서결정되고있기때문에이범위의공력특성데이터오차의영향이최종성능예측에반영되었다고할수있다. TSR=6, 10의경우실험결과와수치해석결과가서로잘일치하는특징을나타내고있으며, 이범위에서성능예측에참조되는익형공력특성예측정도가높기때문이라사료된다. 112
Lift Coefficient 1.35 1.20 1.05 0.90 0.75 0.60 0.45 0.30 0.15 0.00-0.15-0.30-0.45-0.60-8 -4 0 4 8 12 16 20 24 28 32 36 40 44 48 52 56 60 64 68 72 76 80 84 88 92 Angle of Attack Predicted, X-FOIL Post stall correction RE=1,000,000 Fig. 4.14 Predicted lift coefficient (S809) 113
1.2 1.1 1.0 0.9 Drag Coefficient 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0.0 Predicted, X-FOIL Post stall correction RE=1,000,000-8 -4 0 4 8 12 16 20 24 28 32 36 40 44 48 52 56 60 64 68 72 76 80 84 88 92 96 Angle of Attack Fig. 4.15 Predicted drag coefficient (S809) 114
1.1 1.0 0.9 0.8 Lift Coefficient 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0.0 Predicted, X-FOIL Measured, Delft Univ. of Tech., Sommers RE=1,000,000-3 -2-1 0 1 2 3 4 5 6 7 8 9 10 11 12 Angle of Attack Fig. 4.16 Comparison of lift coefficient (S809) 115
0.030 0.028 0.026 0.024 Measured, Delft Univ. of Tech., Sommers Predicted, X-FOIL RE=1,000,000 Drag Coefficient 0.022 0.020 0.018 0.016 0.014 0.012 0.010 0.008 0.006-2 -1 0 1 2 3 4 5 6 7 8 9 10 11 12 Angle of Attack Fig. 4.17 Comparison of drag coefficient (S809) 116
0.65 Axial Flow Induction Factor 0.60 0.55 0.50 0.45 0.40 0.35 0.30 0.25 TSR=6.0 TSR=7.0 TSR=8.0 TSR=9.0 0.20 0.15 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 r/r (a) Comparison of axial flow induction factors (by measured lift&drag) 0.60 0.55 Axial Flow Induction Factor 0.50 0.45 0.40 0.35 0.30 0.25 0.20 0.15 0.10 TSR=2.5 TSR=3.0 TSR=3.5 TSR=4.0 TSR=4.5 TSR=5.0 TSR=5.5 TSR=6.0 TSR=7.0 TSR=8.0 TSR=9.0 0.05 0.00 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 1.1 1.2 1.3 r/r (b) Comparison of axial flow induction factors (by predicted lift&drag) Fig. 4.18 Distributions of axial flow induction factors 117
0.08 Tangential Flow Induction Factor 0.06 0.04 0.02 TSR=6 TSR=7 TSR=8 TSR=9 0.00 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 r/r 0.10 (a) Comparison of tangential flow induction factors (by measured lift&drag) Tangential Flow Induction Factor 0.08 0.06 0.04 0.02 TSR=2.5 TSR=3.0 TSR=3.5 TSR=4.0 TSR=4.5 TSR=5 TSR=5.5 TSR=6 TSR=6.5 mu vs 7.0 TSR=7.5 TSR=9 0.00 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 r/r (c) Comparison of tangential flow induction factors (by predicted lift&drag) Fig. 4.19 Distributions of tangential flow induction factors 118
0.50 0.45 0.40 Power Coefficient 0.35 0.30 0.25 0.20 0.15 0.10 0.05 0.00 1 2 3 4 5 6 7 8 9 10 TSR CP(BEMT), Predicted airfoil data CP(BEMT), Measured airfoil data TSR 6,7,8,9 Fig. 4.20 Comparison of the power coefficient on the FIL-20 119
제 5 장 CFD 를이용한풍력터빈전산해석 시뮬레이션을통해복잡하고다양한물리현상을정확하게구현할수있는전산해석환경이보편화되면서산업계전반에걸쳐전산유체역학 (CFD, Computational Fluid Dynamics) 의적용이활발히진행되고있으며, CAE 산업의한분야로서중요한역활을담당하고있다. 본격적인풍력발전산업에대한전산유체역학적용은 2001년부터이며, 주로 NREL, RISO, NASA 등의유수한연구기관을통해연구가진행되고있다. 현재풍력발전산업에서설계및성능해석의표준으로서 BEMT를이용하는것이보편화되어있으나, 블레이드요소간의독립성을가정하는등의불완전한가정들에의한문제점을안고있다. 또한, 단순한성능해석결과값에대한정보만을제공할뿐이므로, 성능향상및감소에영향을미치는 3차원유동장의정확한모사가불가능하고, 정확한성능예측을위해서는정확한 2차원익형공력특성의확보가필수적이라는것이다. 이에비해, 상대적으로복잡한 Reynolds Averaged Navier-Stokes 방정식에기반을둔 CFD의경우이론적으로명확한해석이가능하고, 다양한물리적변수를표현할수있다는장점이있기때문에미소한형상의변화및유동장의변화에따른성능변화를예측할수있다. 또한, 최근적용가능성을보이고있는 FSI(Fluid Structure Interaction) 기법의도입을통해로터블레이드의비정상하중특성및소음예측까지가능하게되었다. CFD는실제복잡한해석과정과비정상외부유동장에대한방대한계산량등의문제로적용범위가제한되어있지만, 향후, BEMT를대신하여풍력발전용터빈설계및성능평가의표준으로자리매김할것으로예상된다. 120
5.1 수치해석기법 본연구에서는유체기계해석에탁월한성능을보이는 CFX-TASC flow 의솔버기술과다상유동및연소, 화학반응등의해석에탁월 한 CFX-4 의다양한물리모델을결합하여새롭게개발된 fully implicit pressure based AMG coupled solver 인 CFX-5 를사용하였다. 일반적으로범용의상용코드들은 SIMPLE 또는 SIMPLEC, Rhie and Chow 방법과같은압력에기초한방정식으로부터만들어져있다. 이 들압력에기초한코드들은일반적으로다양한물리적모델들과경 계조건을제공하고다른 CAE 도구들과연동을포함하는복잡한 multi-physics 문제등에적용될수있다. 압축성유동의해석에많 이사용되고있는코드들은일반적으로오일러방정식의해법을발 전시킨방법으로부터개발되어져온밀도에기초한공식으로부터만 들어지며충격파해석등에최적화되어있다. 유체기계의정확한해석을위해서는점성저층 (viscous sublayer) 영역 을안정적으로처리할수있는최적화된난류모델들이필요하다. 현 재상용코드에서오랫동안사용되어온 k ε 모델과벽함수의조합 은유체기계에서요구되는높은해의정밀도를만족시킬수없다. 보 다정확한해를구하기위해점성저층에대한해석이요구되고, 이를 만족시키기위해서는벽면근처에종횡비 (aspect-ratio) 가매우큰높은 격자밀집도를확보해야만한다. 이러한요구조건을만족하는 CFX- 5 의수치기법의핵심은질량과운동량방정식의연동화 (coupled formation) 이다. CFX-5 는압력기반유한체적법 (finite volume method) 을 fully implicit 이산화하여얻어지는방정식을 algebric multigrid coupled solver 를이용하여해석한다. SIMPLE 등고전적인 segregated 접근방법 에비해 implicit coupling 방법은수렴을가속시키고, 압축성유동에 있어서수렴성의난점을피할수있고, 높은종횡비의격자를다룰 121
수있기때문에필수적인기법이다. 유체기계에서난류모델의적용에있어아주간단한난류모델도충 분한정확도를제공한다는입장과해석정확도의확보를위해서는 가장발달된난류모델과천이모델을사용해야한다는견해가있다. 이러한견해의차이는엔지니어가해석하는유체기계가서로다르고 요구하는정밀도에대한만족범위가다르기때문이다. 많은경우, 특히유체기계의설계점영역에서는점성과난류효과는 단지전체손실에추가적인기여를하며, 이러한유동장에서는난류 모델이나천이모델의정밀도가상대적으로중요하지않다. 그러나, 1- 방정식모델이나, 2- 방정식모델은계산시간을크게증가시키지않으 므로, 해석의일관성을위해적절히사용하는것이좋다. 반면, 풍력터빈과같이탈설계점영역에서운전되는경우의해석 을수행하는경우 tip vortex, wake, blade stall 등과같은 2 차유동에대 한전체시스템의출력이큰영향을받게되고경계층천이에의한 블레이드공력특성변화가발생하므로이들유동장에대한정확한 예측을위해완전한물리적모델과고난도의수치해석기법이요구 된다. 이러한 2 차유동의예측을위해 2 방정식모델중 k ω 계열 의모델이유동박리현상의지연된예측을피할수있고상대적으로 정확한예측이가능하다고알려져있다. 그러나, k ω 계열의난류 모델도천이영역에대한해석이불가능하기때문에블레이드의전연 (leading edge) 부분의천이가문제가되는해석에는 RSM 모델과 SST 모델에기초하여새롭게발전된 DES(Detached Eddy Simulation) 모델의 적용이필수적이나, 비정상해석을통한적용이가능한단점이있다. 따라서, 본연구에서는다양한난류모델중블레이드주위로발달하 는 2 차유동에대한비교적정확한예측이가능하다고알려져있는 k ω SST 모델을적용하였다. 122
5-1-1 지배방정식 일반적인유동에서다루는운동방정식은질량, 운동량, 에너지방 정식이다. ρ + ( ρu j ) = 0 t x j (5.1) P τ ij ( ρu ) + ( ρu ρu ) = + + S t x x x j j i ui j i i (5.2) T ( ρh P) + ( ρu H ) = ( k ) ( u τ ) + S t x x x x j j ij T j j j j (5.3) 여기서, ρ = 밀도, u i = 속도, P = 압력, µ = 점성계수, H = 전엔탈피, h = 정적엔탈피, T = 온도, τ ij = 응력텐서 상태방정식 ρ = ρ( P, T ), 점성계수와변형률의함수로서응력텐서 의응력-변형률관계, 그리고, h = h( P, T ) 가이들방정식에추가된 다. 난류유동은와점성 (eddy viscosity) 이추가되고방정식은레이놀즈 평균물리량에대해푼다. 방정식의일반적인형태는동일하다. 123
5-1-2 이산화방법 CFX-5 는 implicit pressure based 방법을사용하며, 여기서사용되는 주요독립변수들은 ( P, ui, H ) 이다. 일반적으로, 범용상용코드는위의 변수에대해해석을수행하며, 이는현장에서일반적으로부딪히는 비압축성유동의해석이용이하기때문이다. Fig. 5.1 에서실선으로나타난것은일반적으로알려져있는격자 즉, cell 이다. 실선으로나타낸 cell, 즉 element 를다시나누어점선으 로표현된 sub-element 를구성하며제어체적은음영으로나타낸부분 과같이 node 를둘러싼 sub-element 들로구성되며 hex, tet, wedge, pyramid 등모든 element 형태에대해동일하게적용된다. 모든변수 값과유체의물성치는이노드에저장된다. [38] 수치해석의정확도는노드값으로표현되는적분점 (IP) 들에서의표 면적분 (fluxes) 값들의정확도에의해결정된다. 계산을통해얻어진 해는격자노드에저장되지만, 방정식의대류항, 확산항, 압력구배항 등의다양한항들은적분점에서의해나해의구배값을필요로하며 따라서, element 내부에서의해의변화를계산하기위해 finite element shape function 이사용된다. 이러한방식을 FEM based FVM 혹은 element based FVM 이라한다. Fig. 1 과같이제어체적면에서의적분점 의개수가 2 차원인경우일반적인 FVM 의 4 개에비해 8 개로 2 배가 많은것을알수있다. 3D 육면체격자의경우 6 개에서 24 개로, 사면 체의경우 4 개에서평균 60 개로적분점이많아지므로비교적성긴 격자에대해서도해의정확도가뛰어난장점이있다. 식 (5.1), (5.2), (5.3) 의방정식들을제어체적에걸쳐적분함으로써 질량, 운동량, 그리고에너지방정식에대한이산화적분식은다음과 같다. 0 ρ ρ ρv ( ) + ( ρu j n j ) ip = 0 (5.4) t ip 124
U U u u ρ V ( ) m ( u ) ( P n ) ( µ ( ) n ) S V 0 i i i j + ip i ip = i ip + eff + j ip + ui t ip ip ip x j xi (5.5) 0 0 ( H P / ρ) ( H P / ρ) T ρv ( ) + mip Hip = ( keff n j ) ip + SφV t x ip ip j (5.6) 여기서, ( ) n j ip 는적분점위치에서국부표면벡터이다. 그리고, 제 한체적의적분점표면을통과하는 m ip 는질량유동이다. 모든방정식들 은시간간격의제한을피하기위하여 implicit하게다루어지며, 비정상항에는 1차와 2차 backward Euler 방정식이사용된다. 확산항은 element shape function의미분형태로각적분점의위치에서구배계수를계산함으로써결정된다. 대류항은 Upwind, Quick 등몇몇기법에의해평가될수있으나, 기본설정된기법인 high-resolution 기법을사용한다. High-resolution 기법은대류항에대한 2차정확도의 upwind biased approach에기초한기법이며 Barth와 Jesperson에의해기술된방법과유사하다. φ = φ + β ( φ) r uuv (5.7) ip P ip ip 식 (5.8) 과같이 divergence 형태에서모든항들에대해질량 divergence 항은표면적분의형태로변환된다. m = ρ u n (5.8) ip ip j, ip j, ip 125
밀도는다른대류항처럼표준 high resolution 스킴을적용하여계 산된다. ρ = ρ + β ( ρ) r uuv (5.9) ip P ip ip 이 upwind biased 평가는운동량과에너지방정식의다른대류량과 마찬가지로유동이상당히압축성이어도안정적이며, 2 차의정확도를 가진다. Implicit 방법에서중요한것은 ρu 의선형화이다. 먼저 ρu 는 Newton-Raphson 선형화에의해확정된다. ( ρu) ρ u ρ u ρ u n 0 0 n 0 0 + (5.10) 여기서위첨자 n 은새로운값 (implicit) 을의미하고 0 는예전 ( 지연 된값 ) 시간레벨이다. 이러한선형화는전영역에걸친마하수의신 뢰성있는수렴을보장한다. 마지막으로, 식 (5.11) 과같이밀도에대한상태방정식은압력의항 으로구성된 된상태방정식에서미분항 n ρ 의 implicit 표현을얻기위하여차분되며앞서제공 n 0 ρ n 0 ρ = ρ + ( P P ) P ρ P 를계산한다. (5.11) 126
Fig. 5.1 Mesh arrangement and terminology for dual mesh 127
5-1-3 난류모델링 유동해석을수행할때가장큰에러의원인중의하나는난류모델 의부적절한사용이라할수있으며, 특히, 벽면근처의격자생성에 있어모든영역에 y + 를일정한수준으로유지한다는것은 3 차원유 동장의경우상당히어려운작업이다. Wilcox model 의벽근처방정식 에는부가적인 viscous sub-layer damping 함수가필요치않다. 일반적으로 Wilcox model 의단점으로자유유선에민감한결과를보 이는것을들수있는데 CFX-5 에서는이러한단점을보완하여벽면 근처에서는 k ω 모델을사용하고바깥쪽은 k ε 모델을사용하는 BSL(Baseline Model) 과 SST(Shear Stress Transport) 모델을지원한다. ω 모델의또다른장점은쉽게자동처리벽처리법 (automatic wall treatment) 로확장이가능하다는것이다. 이는가능한격자의 y + 에무 관하게해의정확성을확보하기위한것이다. 표준 viscous sub-layer model 들이벽면전단응력을정확히해석하기위해 y + 1 의수준을 요구하는반면자동벽면처리기법은성긴벽면격자를처리할수 있는장점이있다. 유체기계유동장은상당히복잡한형태이므로이 러한자동벽면처리조건은상당히유용한기능이다. k ω SST 모델은난류전단응력의수송 (transport) 을계산하기때문 에역압력구배에의해발생하는유동박리크기와발생시점을정확 히예측할수있다. Wilcox 모델과 k ε 모델의장점만을취해 BSL 모델이개발되었으나, 매끄러운표면에발생하는유동박리시점및 크기에대한정확한예측에실패하였다. 이러한원인에대한상세한 내용은 Menter 의연구결과에상세히기술되어있다. 가장주된원인 으로서이전의난류모델들은모두난류전단응력의수송에대한고려 를하지않았기때문이며, 그결과 eddy-viscosity 에대한과다예측을 하였다. 수송항은식 5.12 와같이 eddy-viscosity 형태의방정식에대한 제한으로얻어질수있다. 128
ν a k 1 t = (5.12) max( a1ω, SF2 ) 여기서, ν = µ / ρ, F 2 =blending function, S=strain rate t t Blending function은난류모델의성공을위해매우중요한요소이다. 이방정식의형태는표면과의가장가까운거리와유동변수를기반으로한다. F = tanh(arg ) (5.13) 4 1 1 k 500ν 4ρk arg1 = min max,, ' 2 2 β ω y y ω CDkwσ ω 2 y (5.14) 여기서, y 는벽면으로부터가장가까운곳까지의거리를의미한다. ν 는동점성계수이다. CD kω 1 = max 2 ρ k ω,1.0 10 σω 2ω 10 (5.15) F = tanh(arg ) (5.16) 2 2 2 2 k 500ν arg 2 = max, ' 2 β ω y y ω (5.17) SST 모델이나 BSL 모델은 k ε 과 k ω 사이의 blending을위해 129
벽면과가장가까운거리에위치한노드의거리정보를필요로한다. Wall scale 방정식은식 (5.18) 과같은단순한형태의방정식으로부터 구할수있다. 2 φ = 1 (5.18) 여기서, φ 는 wall scale 값을의미한다. 벽면거리는식 (5.19) 에의해 wall scale 로부터계산되어진다. Wall Distance = 2 ( φ + 2 φ) φ (5.19) 130
5.2 계산조건 (FIL-1000) 본연구에서는컴퓨터하드웨어의제한적인문제로인하여 3 차원 정상상태유동장에대한해석을수행하였으며, Navier-Stokes solver 로서 CFX 5.7 을이용하였다. 로터블레이드는 120 간격의주기조건을만족하기때문에하나의블레이드만을해석영역으로설정하였다. 너셀과타워에의한블레이드의공력특성변화는무시하였다. 입구조건으로서, 균일속도유입조건을적용하였으며, 입구속도의변화에따른계산을수행하였다. 회전속도는정격회전속도를유지하고입구유입속도를 5.76 m/s 에서 36.48 m/s 까지변화시키며계산을수행하였다. Table 5.1 에자세한계산조건을요약정리하였다. 131
Table 5.1 Simulation conditions TSR Wind speed (m/s) Rotating speed (rpm) 3 24.98 26.28 4 18.74 26.28 5 14.99 26.28 6 12.49 26.28 7.5 9.99 26.28 9 8.33 26.28 10 7.49 26.28 11 6.81 26.28 12 6.24 26.28 13 5.76 26.28 132
5.3 계산격자및경계조건 (FIL-1000) 계산격자의생성은적용난류모델의특성에따라원활한수렴및신뢰성있는결과를확보하기위해 y+, 경계층격자밀집도, 격자형태, aspect ratio 등을신중히고려해야만한다. 따라서, 우수한품질의계산격자의확보가 CFD 에서첫번째필수적인요소라할수있다. 그러나, 단일 CPU 에의한계산환경에서는하드웨어성능의제한에의한충분한격자공간해상도를확보하기가상당히까다롭다. 본연구에서는효율적인계산격자의생성을위하여블레이드를포함하는내부영역과블레이드를포함하지않는외부영역으로격자계를이중분할하였다. 블레이드에서발생하는토오크의정확한예측을위해블레이드주변영역을 multi-block unstructured hexa 격자계로구성하였고, 나머지영역은 tetra-prism 격자계로구상하였으며, 두영역의경계면은 pyramid 격자로처리하였다. 전체계산격자는 hexa-tetra-prism-pyramid 형태의다양한격자를사용해이루어진 hybrid type 격자계로구성되었다. 일반적으로 k-ε 모델은유동박리현상이지배적인유동장의예측에있어해석결과의정도가낮다고평가되고있으므로, 로터블레이드 표면으로부터발생되는실속현상등을포함하는복잡한 3 차원 유동현상을파악하기위한적용난류모델로서적합하지않다. 따라서, 본연구에서는 k-ω SST 난류모델 (Menter, 1993) 을적용하였다. [39] 블레이드벽면경계로부터첫번째지점격자까지거리의척도인 y + 는 20 이하로제한하였으며, 익형현의길이방향을따라모두 180 노드의격자를분포시키고로터블레이드반경방향으로 42 노드를분포시켰다. 133
로터블레이드를포함하는도메인의전체격자수는약 480,000 노드이며, 로터블레이드를포함하지않는외부도메인의격자수는약 250,000 노드이다. 전체계산격자수는약 730,000 노드를사용하였으며, ICEM-CFD 5.0 을이용하여생성하였다. 경계조건으로외부도메인출구영역에 averaged static pressure 조건을부여하였고, 단일블레이드계산을위하여양쪽주기경계면에 periodic 조건을부여하였다. 로터블레이드로부터입구까지의거리는로터직경의 3배, 윗면으로 5배, 후방으로 7배를확보하였으며, 로터블레이드벽면은 no-slip 조건으로처리되었고, 상대회전조건을부여하였다. 입구유입조건으로서균일한흐름이유입된다고가정하였다. 계산은단일 CPU (Pentium4 3.05GHz) 로수행되었으며, 2GB RAM 을사용하였다. 모든계산은 100 회반복계산이전에수렴이되었으며, 수렴판단조건은 10-5 이다. 단일케이스에대한계산수행시간은최적화된상태로부터대략 5~7 시간정도소요되었다. Fig. 5.2, Fig. 5.3 에계산격자의다양한형태를각각나타내었다. 134
Fig. 5.2 O type grid distributions for a change of local blade radius 135
(a) Hexa mesh(near by the rotor) (b) Tetra-prism mesh(outer domain) Fig. 5.3 Computational grids 136
5.4 결과및고찰 (FIL-1000) 5.4.1 블레이드표면유선로터블레이드주위로발생하는복잡한 3 차원유동현상등의정확한이해를위해 TSR 의변화에따라표면압력분포, 표면유선, 다양한중첩된가시화결과를나타내었다. Fig. 5.3, Fig. 5.4, Fig. 5.5 에 TSR 3, 7.5, 10 에서로터블레이드흡입면 (suction side) 의표면유선 (streamline) 을나타내었다. 항공기날개에서는받음각의변화에따라양력, 항력의변화가발생하게되고, 특정받음각이상으로유지되는경우, 날개의전연부터박리가진행되며깊은실속에빠지게된다. 대부분의실속을발생시키는주된원인이받음각의변화이다. 그러나풍력발전용로터블레이드는항공기날개와는달리빠른속도로회전하는상태이므로, 원심가속력 (centrifugal acceleration) 과로터허브와팁사이의압력차가발생하게된다. 로터블레이드허브근방영역은비교적높은받음각에서운전되기때문에박리에의해흐름이표면으로부터이탈하기쉽다. 일단표면으로부터이탈된흐름은블레이드후방으로빠져나가지못하고, 로터의회전으로부터발생하는원심가속력과압력차에의해블레이드표면을따라팁방향으로진행하는특징을보인다. 이를반경류 (radial flow) 라하며받음각의변화와함께풍력터빈용블레이드에서실속을초래하는중요한요소가된다. [40] Fig. 5.4 는입구속도가 24.98 m/s 인경우에서블레이드표면의흐름을나타내고있다. 로터블레이드의회전방향속도성분과유입방향속도성분의합벡터는바람이블레이드전연과이루는각도인상대유입각도라생각할수있다. 137
계산조건에서로터회전속도는고정되어있으므로유입풍속이증가함에따라상대유입각도또한증가하게된다. 따라서, 비교적입구유입풍속이높은조건에서받음각의증가에의한실속현상이발생되고있음을알수있다. Fig. 5.5 는 Fig. 5.4 에비해비교적안정된형태의유동장을나타내고있으며, 상대적으로높은받음각이형성되는허브근방을제외한영역에서는실속현상을관찰할수없다. 이때의입구풍속은 9.99 m/s 이며설계 TSR 조건이다. Fig. 5.6의경우입구풍속이 7.49 m/s 일때의결과를나타내고있으며, 반경류에의한영향이지속적으로감소하고있음을제외하고는 Fig. 5.5의결과와거의유사한형태의유동장을관찰할수있다. 138
Fig. 5.3 Surface streamline V in =24.98 m/s 139
Fig. 5.4 Surface streamline V in =9.99 m/s 140
Fig. 5.5 Surface streamline V in =7.49 m/s 141
5.4.2 블레이드국부단면흐름특성 Fig. 5.6, Fig. 5.7, Fig. 5.8 에로터블레이드국부단면에서의흐름변화상태를입구풍속의변화에따라나타내었다. 절단단면은허브로부터각각 5m, 10m, 15m 이다. Fig. 5.6은허브로부터 5m 위치단면의결과를나타낸다. 입구풍속이 24.98 m/s인경우, Fig. 5.3의표면유선결과로부터이미블레이드전체에걸쳐실속이진행되었음을알수있다. 입구풍속이 9.99 m/s인경우전연박리현상은없으며, 박리점이후연쪽으로상당부분후퇴되었음을알수있다. 입구풍속이 7.49 m/s 인경우거의후연끝단에이르러서야미소한박리가발생하고있다. 익형의전연박리는깊은실속을뜻하며, 양력의증가보다항력의증가폭이급격히상승하기때문에공력특성은현저히저하되고, 이는결국전체시스템의출력저하로확장됨을예상할수있다. 이러한현상은로터블레이드가정격회전속도이상으로회전할때적절한출력을유지하기위한출력제어도구로사용될수있으며, 이를실속제어방식이라한다. Fig. 5.7은허브로부터 10m 위치단면의결과를나타낸다. 입구풍속이 24.98 m/s인경우, 5m 단면에서와마찬가지로블레이드는깊은실속에빠져있다. 9.99 m/s의경우에서는후연끝단에미소한박리가관찰되고있으나, 7.49 m/s의경우블레이드주위의흐름이완전한부착류를형성하고있음을알수있다. Fig. 5.8은허브로부터 15m 위치에서의결과를나타낸다. 입구풍속이 24.98 m/s인경우, 5m, 10m 단면에서와마찬가지로블레이드는깊은실속에빠져있으나, 상대적으로흡입면에서관찰되는 실속의크기는앞의두경우보다작다. 9.99 m/s 의경우에서는후연 끝단에미소한박리가관찰되고있으나, 7.49 m/s 의경우블레이드 주위의흐름이완전한부착류 (attached flow) 를형성하고있음을알 142
수있다. 일반적으로받음각의범위가낮아질수록익형주위흐름은부착류를형성하는특징이나타내기때문에완전한부착류의형성이최고의공력특성을나타내는지점이라생각할수없다. 실속은실속유도진동등의불필요한위험요소를발생케하는문제가있으므로, 가급적억제하는것이바람직한방법이나, 출력제어용으로실속현상을이용하기도한다. Fig. 5.9 에입구풍속이 24.98 m/s 일때의로터블레이드표면에서속도벡터와유선등의결과를동시에나타내었다. 블레이드표면을따라반경방향으로진행하는반경류의형성을표면유선으로부터확인할수있으며, 이반경류는블레이드표면을파고들면서표면으로부터흐름을이탈시키는역할을하고있음을알수있다. 로터 블레이드 2 차원단면의속도벡터를나타낸부분을살펴보면 블레이드회전방향속도성분이지배적임을알수있으나, 블레이드표면근처에서는반경방향속도성분또한크게나타나고있음을알수있다. Fig. 5.10은입구풍속 24.98 m/s의날개끝와류발생을나타내고있다. 날개끝와류의발생은순환분포의감소에의한손실의형태로표현되므로이러한현상을줄이기위하여다양한형태의날개끝형상의변형이나, 보조장치의장착등을고려할필요가있다. Fig. 5.11에입구풍속 9.99 m/s인경우의후류분포를나타내었다. 143
(a) V in =24.98 m/s (b) V in =9.99 m/s (c) V in =7.49 m/s Fig. 5.6 Sectional streamline 5m from root 144
(a) V in =24.98 m/s (b) V in =9.99 m/s (c) V in =7.49 m/s Fig. 5.7 Sectional streamline 10m from root 145
(a) V in =24.98 m/s (b) V in =9.99 m/s (c) V in =7.49 m/s Fig. 5.8 Sectional streamline 15m from root 146
Fig. 5.9 Overlapped visualization results on the rotor surface 147
Fig. 5.10 Tip vortices generation - V in =24.98 m/s 148
Fig. 5.11 Visualization of 3D wake geometry - V in =9.99 m/s 149
5.4.3 출력특성 Fig. 5.12에입구풍속의변화에따른출력특성변화를나타내었다. 출력성능의비교를위해 BEMT 해석결과와비교하였으며, 두결과는비교적잘일치하고있음을알수있다. BEMT 해석의경우 POSEIDON 에서 TSR의변화에따라구한출력계수와 CFD의계산조건에해당하는입구풍속을식 2.13에대입하여출력값을구하였다. CFD 해석의경우입구풍속의변화에따라계산이수행되기때문에각각의계산조건에서블레이드에작용하는토오크를구한후식 2.19의적분형태로출력을계산하였다. 150
Fig. 5.12 Comparison of the power between CFD and BEMT 151
5.5 계산조건 (FIL-20) CFD 를이용한소형풍력발전용터빈의 3 차원유동해석및성능평가에관한정상상태해석을수행하였다. Navier-Stokes solver 로써 CFX 5.7 을이용하였으며, FIL-1000 의해석과마찬가지로효율적인계산격자계의구성을위해 120 간격으로주기조건을적용하고단일블레이드에대한해석영역을설정하였다. 너셀과타워에의한블레이드의공력특성변화는무시하였다. 입구조건으로서, 균일속도유입조건을적용하였으며, 입구속도를 6 m/s 에서 14 m/s 까지변화시켰으며, 회전속도는설계시설정된정격회전속도인 90.49 rpm 로고정되었다. Table 5.2 에자세한계산조건을요약정리하였다. 152
Table 5.2 Simulation conditions FIL-20 TSR Wind speed (m/s) Rotating speed (rpm) 8 6 90.49 6 8 90.49 4.8 10 90.49 4 12 90.49 3.42 14 90.49 153
5.6 계산격자및경계조건 (FIL-20) 계산격자의생성을위해 ICEM-CFD 5.1을이용하였으며, 블레이드를포함하는내부영역과외부영역으로 2분할하여격자생성의효율성을도모하였다. 계산격자는내부영역에약 505,000 노드를분포시켰으며 hexa 타입의격자이다. 외부영역은 tetra-prism 타입의격자계로구성되었으며약 200,000 노드의계산격자로구성된다. 내부영역과외부영역의경계면에대한원활한격자매칭을위해 pyramid 타입격자를적절히분포시켰다. 전체계산격자수는약 705,000 노드이며 hybrid 타입이다. Fig. 5.13에블레이드를포함하는내부영역격자구성도를나타내었으며, Fig. 5.14에외부영역에대한격자구성도를나타내었다. 블레이드벽면경계로부터첫번째지점격자까지거리의척도인 y + 는 20 이하로제한하였으며, 익형현의길이방향을따라모두 180 노드의격자를분포시키고로터블레이드반경방향으로 42 노드를분포시켰다. 로터블레이드를포함하는도메인의전체격자수는약 480,000 노드이며, 로터블레이드를포함하지않는외부도메인의격자수는약 250,000 노드이다. 전체계산격자수는약 730,000 노드를사용하였으며, ICEM-CFD 5.0 을이용하여생성하였다. 경계조건으로서외부도메인출구영역에 averaged static pressure 조건을부여하였고, 단일블레이드계산을위해주기경계면에대해 periodic 조건을부여하였다. 입구경계면으로부터로터블레이드압력면까지의거리는로터직경의 3배, 블레이드팁부분부터윗면까지 5배, 블레이드흡입면으로부터출구경계면까지 7배의거리를확보하였으며, 벽면은 no-slip 조건으로처리되었고, 상대회전조건을부여하였다. 입구유입조건으로서균일한흐름이유입된다고가정하였다. 154
계산은단일 CPU (Pentium4 3.05GHz) 로수행되었으며, 2GB RAM 을사용하였다. 모든계산은 70 회반복계산이전에수렴이되었으며, 수렴판단조건은 10-5 로설정하였다. 단일케이스에대한계산수행시간은최적화된상태로부터대략 5~7 시간정도소요되었다. 155
Fig. 5.13 Inner computational domain of the FIL-20 156
Fig. 5.14 Outer computational domain of the FIL-20 157
5.7 결과및고찰 5.7.1 블레이드표면유선 로터블레이드주위로발생하는복잡한 3 차원유동현상등의 정확한이해를위해유입풍속의변화에따라표면압력분포, 표면유선, 다양한중첩된가시화결과등을나타내었다. Fig. 5.15, Fig. 5.16, Fig. 5.17 에입구풍속 6m, 8m, 10m 인경우로터블레이드흡입면의표면유선변화를나타내었다. Fig. 5.15 는입구속도가 6 m/s 인경우에서블레이드표면의흐름을나타내고있다. 전체적으로낮은입구풍속에대한영향으로받음각의변화폭이적기때문에허브근방영역을제외하고는반경류의형성에의한실속현상이나타나고있지않다. Fig. 5.16 의경우입구풍속이 8 m/s 로증가하였으며, 6 m/s 의경우에비해비교적넓은영역에걸쳐반경류가형성되고있음을뚜렷이확인할수있다. 그러나, 실제로터블레이드의출력대부분을발생시키는허브로부터 70% ~ 100% 의영역에대해서는큰영향을미치질못하고있다. Fig. 5.17 의경우입구풍속이 10 m/s 로증가한경우를나타내며, 거의전영역에걸쳐반경류가형성되고있음을알수있다. 이러한반경류에의해블레이드전영역에걸쳐후연박리가진행되고있음을알수있다. 158
Fig. 5.15 Surface streamline V in = 6 m/s 159
Fig. 5.16 Surface streamline V in = 8 m/s 160
Fig. 5.17 Surface streamline V in = 10 m/s 161
5.7.2 블레이드 3차원유동특성 Fig. 5.18, Fig. 5.19, Fig. 5.20 에입구풍속의변화에따라로터블레이드국부단면에서의흐름변화상태를나타내었다. 절단단면은허브로부터각각 1m, 2m, 5m 의위치이다. Fig. 5.18은허브로부터 1m 위치단면의결과를나타낸다. 1 입구풍속이 6 m/s인경우, Fig. 5.15의표면유선결과로부터이미블레이드허브근방특정영역을제외하고반경류의형성이확인되지않았으므로전체적인흐름은부착류의형태를보이고있다. 입구풍속이 8 m/s인경우전연박리현상은없으며, 후연부근에박리가시작되고있음을알수있다. 이는 Fig. 5.16의결과의다소넓은영역에걸친반경류의형성과관계가있다. 입구풍속이 10 m/s 인경우, 이전의결과와동일하게후연박리가진행되고있으며, 박리점의위치가 6 m/s, 8 m/s의결과보다더전연쪽에근접하고있다. Fig. 5.19는허브로부터 2m 위치에서의단면결과를나타낸다. 입구풍속이 6 m/s인경우, 1m 단면에서와마찬가지로블레이드는거의완전한부착류를형성하고있으며. 반경류의형성이거의없다. 8 m/s의경우에서는후연끝단에반경류와받음각의변화에따른미소한박리가관찰되고있으나, 10 m/s의경우여전히후연박리가강하게진행되고있음을알수있다. Fig. 5.20은허브로부터 5m 위치에서의결과를나타낸다. 입구풍속이 6 m/s인경우, 1m, 2m 단면에서와마찬가지로블레이드주위의흐름은완전한부착류를형성하고있으며, 8 m/s, 10 m/s의경우도마찬가지로반경류에의한영향을거의받고있지않음을알수있다. Fig. 5.21에입구풍속이 8 m/s일때로터블레이드표면단면 1m, 2m, 3m, 4m에서의단면유선분포를나타내었다. 162
Fig. 5.22에는반경류의형성및 2차원단면에의영향을살펴보고자허브로부터 1.5m위치에서 2차원단면의속도벡터와블레이드표면유선을동시에나타내었다. 속도벡터는 2차원단면에대해 3차원성분을갖도록표현하였다. 익형의후연근방에서시작되는와류영역을살펴보면블레이드팁방향으로의 3차원속도성분이존재하고있음을확인할수있다. 이 3차원속도성분은블레이드표면에의유선의방향과크기가정확히일치하며, 결국표면을따라흐르는반경류의영향이받음각의변화와함께블레이드실속의주된원인이됨을잘알수있다. Fig. 5.23은입구풍속이 10 m/s인경우에서의 3차원유선을나타내고있다. 블레이드허브끝단으로부터두갈래의강한소용돌이와류가형성되고있으며, 이흐름이블레이드표면을따라진행하게된다. 이두갈래의와류는블레이드지지형상이원형에가까운형상이므로원주후류에발생하는대칭형와류의 3차원성분이다. Fig. 5.24에블레이드흡입면에서의압력분포를나타내었다. 회전유동장이기때문에블레이드팁부분에서의흡입면압력이상대적으로낮게형성되고있다. Fig. 5.25에입구풍속 8 m/s인경우의 3차원후류진행형태를나타내었다. 163
(a) V in = 6 m/s (b) V in = 8 m/s (c) V in = 10 m/s Fig. 5.18 Sectional surface streamline 1m from root 164
(a) V in = 6 m/s (b) V in = 8 m/s (c) V in = 10 m/s Fig. 5.19 Sectional surface streamline 2m from root 165
(a) V in = 6 m/s (b) V in = 8 m/s (c) V in = 10 m/s Fig. 5.20 Sectional surface streamline 5m from root 166
Fig. 5.21 Sample of section view (V in = 8 m/s, 1m, 2m, 3m, 4m sections) 167
Fig. 5.22 Sample of section view (V in = 10 m/s, 1.5m section from root) 168
Fig. 5.23 3D stall on the rotor surface (V in = 10 m/s) 169
Fig. 5.24 Suction side surface pressure distribution 170
Fig. 5.25 Visualization of 3D wake geometry - V in =8 m/s 171