熊宇軒, 葉祖洋*
(1.武漢科技大學(xué)資源與環(huán)境工程學(xué)院, 武漢 430081; 2.冶金礦產(chǎn)資源高效利用與造塊湖北省重點(diǎn)實(shí)驗(yàn)室, 武漢 430081)
隨著能源消耗問題日益嚴(yán)峻,能源結(jié)構(gòu)調(diào)整迫在眉睫,風(fēng)能作為一種清潔能源,為推動(dòng)能源可持續(xù)轉(zhuǎn)型,在能源中比重不斷增加[1]。風(fēng)力機(jī)作為一種通過能量轉(zhuǎn)換將風(fēng)能轉(zhuǎn)化為機(jī)械能的發(fā)電設(shè)備,其生產(chǎn)發(fā)電效益受風(fēng)電場(chǎng)中環(huán)境風(fēng)變化規(guī)律的影響。風(fēng)力機(jī)大部分時(shí)間在非穩(wěn)定流動(dòng)環(huán)境中運(yùn)行,受風(fēng)剪切、旋轉(zhuǎn)效應(yīng)、大氣湍流、傾斜和偏航失調(diào)等非定常因素影響的各種流動(dòng)條件影響,導(dǎo)致風(fēng)輪的不穩(wěn)定負(fù)載[2]。風(fēng)力機(jī)故障、服役壽命縮短、輸出功率降低和運(yùn)行維護(hù)增加都與此直接相關(guān)[3]。同時(shí),因?yàn)槿~片流場(chǎng)是不穩(wěn)定的三維分離流,會(huì)導(dǎo)致風(fēng)力機(jī)動(dòng)態(tài)失速,并伴隨著巨大的負(fù)載變化,致使空氣動(dòng)力學(xué)建模和非定常載荷的精確預(yù)測(cè)面臨許多挑戰(zhàn)。
風(fēng)力機(jī)葉片在復(fù)雜環(huán)境下氣動(dòng)性能分析是研究整個(gè)風(fēng)電機(jī)組發(fā)電效率的重要方法。對(duì)風(fēng)力機(jī)非定常來流作用下氣動(dòng)特性的研究具有重要意義。李歡等[4]對(duì)風(fēng)力機(jī)進(jìn)行了氣動(dòng)力和離心載荷模擬,并得出兩者對(duì)風(fēng)力機(jī)結(jié)構(gòu)特性的影響程度。張建平等[5]通過函數(shù)擬合不同槳距角下葉片位移和應(yīng)力,給出模擬計(jì)算下葉片動(dòng)力特性變化規(guī)律。代元軍等[6]分析對(duì)比了葉尖結(jié)構(gòu)對(duì)風(fēng)力機(jī)流場(chǎng)計(jì)算流體力學(xué)(computational fluid dynamics,CFD)數(shù)值模擬結(jié)果和風(fēng)洞實(shí)驗(yàn)規(guī)律,證明了兩者的一致性。李德源等[7]采用升力線理論模型和葉素動(dòng)量理論對(duì)風(fēng)力機(jī)葉片氣動(dòng)性能進(jìn)行分析對(duì)比,驗(yàn)證了計(jì)算模型的可靠性。朱呈勇等[8]建立了基于Leishman-Beddoes模型和Bak模型的動(dòng)態(tài)失速和旋轉(zhuǎn)效應(yīng)氣動(dòng)響應(yīng)預(yù)測(cè)模型,有效預(yù)測(cè)了耦合作用對(duì)非定常氣動(dòng)力的影響程度。王占洋等[9]通過考慮流固耦合作用探討風(fēng)力機(jī)應(yīng)力應(yīng)變分布特征,證明了耦合作用下風(fēng)力機(jī)輸出功率更貼近試驗(yàn)值。Cao等[10]研究了不同條件下葉素動(dòng)量理論計(jì)算作用在葉片上的氣動(dòng)荷載動(dòng)態(tài)響應(yīng)分析,驗(yàn)證了理論結(jié)果和實(shí)驗(yàn)結(jié)果的正確性。Li等[11]采用多體動(dòng)力學(xué)代碼集成CFD模擬方法,研究了在大氣非定常情況下尾流風(fēng)湍流和轉(zhuǎn)子相互作用,表明氣動(dòng)載荷和葉尖撓度在時(shí)域和頻域上具有良好的一致性。
針對(duì)風(fēng)力機(jī)非定常影響的氣動(dòng)性能研究,大多只考慮正常來流風(fēng)條件下的影響,現(xiàn)考慮均勻來流和三維旋轉(zhuǎn)效應(yīng)共同作用下非定常條件,結(jié)合ANSYS平臺(tái)流固耦合分析,通過結(jié)構(gòu)場(chǎng)接收流場(chǎng)穩(wěn)態(tài)計(jì)算結(jié)果,傳遞耦合面壓力、速度等數(shù)據(jù)研究風(fēng)輪固體結(jié)構(gòu)在流體作用下氣動(dòng)耦合特性。
流固耦合是通過數(shù)值模擬軟件求解模塊采用弱耦合分離式求解方法按照設(shè)定順序?qū)⑺⒌牧黧w域計(jì)算流體動(dòng)力學(xué)(computational fluid dyna-mics,CFD)方程和固體域結(jié)構(gòu)動(dòng)力學(xué)(computational structural dynamics,CSD)方程計(jì)算數(shù)據(jù)進(jìn)行交互傳輸處理,從而實(shí)現(xiàn)各個(gè)時(shí)間步上的流固耦合求解。流體域模擬采用ANSYS FLUENT,該模塊代碼使用有限體積法求解流體的控制方程。在本文研究中,對(duì)整個(gè)流動(dòng)區(qū)域求解不可壓縮、非定常雷諾平均Navier-Stokes(URANS)方程。選擇基于耦合壓力的求解器采用二階隱式公式,以提高精度,其中解變量采用二階迎風(fēng)離散化格式求解。
質(zhì)量守恒方程為

(1)
動(dòng)量守恒方程為

(2)
式中:t為時(shí)間;ff為體積力矢量;ρf為流體密度;v為流體速度矢量;τf為剪切力張量。
剪切力張量表達(dá)式為

(3)

固體部分系統(tǒng)守恒方程由牛頓第二定律推導(dǎo)為

(4)

計(jì)算耦合交界處存在關(guān)于應(yīng)力、位移等變量的守恒,滿足
(5)
式(5)中:τ為應(yīng)力;d為位移;n為方向余弦;下標(biāo)f表示流體;下標(biāo)s表示固體。
水平軸風(fēng)力機(jī)氣動(dòng)耦合分析采用混合k-ω/k-ε剪切應(yīng)力傳輸(shear stress transfer,SST)湍流模型,k-ω的SST模型可以考慮湍流剪切的傳輸并準(zhǔn)確預(yù)測(cè)不同壓力梯度條件下流量傳輸情況,k-ε模型則用于外部區(qū)域和自由剪切流。利用指標(biāo)殘差監(jiān)測(cè)來確定計(jì)算結(jié)果的收斂性,將流量連續(xù)性、xyz方向速度、湍流動(dòng)能和比耗散率的殘差收斂精度設(shè)置為1×10-6、1×10-3。壓力、動(dòng)量、湍流動(dòng)能、比耗散率和比能量的空間離散均近似于二階。SST模型為

(6)

(7)

混合函數(shù)F1定義為
(8)
葉片截面葉素坐標(biāo)使用Profili軟件中NACA44XX翼型坐標(biāo),通過多線段樣條曲線導(dǎo)入翼型坐標(biāo),在同一軸線上按翼型攻角縮放、扭轉(zhuǎn)葉素曲線,形成葉片基礎(chǔ)骨架,使用鋪層表面完成葉片建模。幾何結(jié)構(gòu)參考1.5 MW NREL水平軸風(fēng)力機(jī)模型,風(fēng)輪直徑為6.5 m,塔筒高度為8.75 m,塔筒頂端截面直徑為0.25 m,底部截面直徑為0.45 m,機(jī)艙通過簡(jiǎn)化設(shè)計(jì)為0.9 m×0.7 m×0.525 m立方體結(jié)構(gòu),如圖1所示。

圖1 風(fēng)力機(jī)三維幾何模型
為了模擬旋轉(zhuǎn)狀態(tài)下的風(fēng)力機(jī),將封閉整個(gè)風(fēng)力機(jī)轉(zhuǎn)子的移動(dòng)圓柱形旋轉(zhuǎn)內(nèi)域與包含機(jī)艙和塔筒的固定長(zhǎng)方體外域相結(jié)合。圓柱形內(nèi)部區(qū)域的直徑和厚度分別為6.8 m和0.6 m。葉片與內(nèi)部旋轉(zhuǎn)域通過布爾運(yùn)算完成內(nèi)流域的創(chuàng)建。模擬均勻來流風(fēng)速的非定常影響并減少邊界條件的影響,采取將外流域中塔筒和內(nèi)流域挖除的方法,整個(gè)計(jì)算區(qū)域的尺寸為40 m×30 m×21 m,如圖2所示。兩個(gè)相對(duì)運(yùn)動(dòng)的區(qū)域共享一個(gè)滑動(dòng)界面,作為數(shù)據(jù)傳輸耦合面。

圖2 風(fēng)力機(jī)計(jì)算域示意圖
2.3.1 CFD模型處理
流場(chǎng)設(shè)置由ANSYS Workbench中FLUENT模塊實(shí)現(xiàn),流場(chǎng)域計(jì)算抑制整體計(jì)算域中結(jié)構(gòu)場(chǎng)部分。葉片表面是一個(gè)無滑移剪切條件的旋轉(zhuǎn)壁面,也是兩個(gè)求解器之間傳遞數(shù)據(jù)的流體域和結(jié)構(gòu)域的邊界。用非結(jié)構(gòu)化網(wǎng)格對(duì)區(qū)域進(jìn)行離散化,并在葉片表面添加更精細(xì)的網(wǎng)格和額外的膨脹層,以提高分辨率。在Meshing中設(shè)置對(duì)稱壁面、地面單元尺寸為2.719 m,內(nèi)旋轉(zhuǎn)域設(shè)置體網(wǎng)格大小為0.1 m,采用曲率捕捉劃分方式,定義曲率法向角為18°。流場(chǎng)網(wǎng)格劃分如圖3所示。

圖3 流體域邊界層細(xì)節(jié)網(wǎng)格
上游速度入口邊界,距離整機(jī)10 m,設(shè)置為自由流風(fēng)速,湍流動(dòng)能強(qiáng)度為0.05,湍流黏度比為10。下游壓力出口邊界,距離整機(jī)30 m,壓強(qiáng)為一個(gè)標(biāo)準(zhǔn)大氣壓。外域固定對(duì)稱壁面、地面和內(nèi)部旋轉(zhuǎn)域葉片表面設(shè)置不可滑移邊界,應(yīng)用Frame motion創(chuàng)建旋轉(zhuǎn)域相對(duì)于靜止域繞z軸旋轉(zhuǎn),旋轉(zhuǎn)軸初始坐標(biāo)為(-4.956 7×10-3,2.541,0.753 5),定義轉(zhuǎn)速為200 r/min。由于模擬過程中葉片與旋轉(zhuǎn)域處于相對(duì)運(yùn)動(dòng)的狀態(tài),對(duì)葉片交界面采用同于旋轉(zhuǎn)域的壁面運(yùn)動(dòng)方式,選擇運(yùn)動(dòng)壁面,運(yùn)動(dòng)類型為旋轉(zhuǎn),相對(duì)轉(zhuǎn)速為0。
2.3.2 有限元結(jié)構(gòu)模型處理
結(jié)構(gòu)場(chǎng)的計(jì)算設(shè)置考慮受到流場(chǎng)數(shù)值計(jì)算壓力數(shù)據(jù),忽略流場(chǎng)體積部分,網(wǎng)格劃分前,為風(fēng)輪和塔筒分別賦予材料屬性,如表1所示。

表1 風(fēng)力機(jī)材料特性參數(shù)
劃分葉片網(wǎng)格時(shí),為了避免葉片前后翼緣轉(zhuǎn)角尖端在劃分網(wǎng)格時(shí)出現(xiàn)單元扭曲,嵌入表面網(wǎng)格劃分單元,分別設(shè)置葉片和塔筒單元尺寸為0.02 m和0.05 m,如圖4(a)所示。因?yàn)槿~片表面是耦合邊界,交換氣動(dòng)力和結(jié)構(gòu)位移,采用大撓度假設(shè),離心力考慮了葉片結(jié)構(gòu)的轉(zhuǎn)速,定義風(fēng)輪旋轉(zhuǎn)方式。求解設(shè)置添加風(fēng)輪旋轉(zhuǎn)軸和幅度值,選擇風(fēng)力機(jī)支架的底部,施加固定約束。在流固耦合和結(jié)構(gòu)求解時(shí)以葉片表面為耦合界面,導(dǎo)入流場(chǎng)壓力數(shù)據(jù),最大壓力值為2 448 Pa,最小壓力值為-3 559.1 Pa,如圖4(b)所示。

圖4 靜力結(jié)構(gòu)分析
為了評(píng)估CFD分析的收斂性,監(jiān)測(cè)殘差、x速度、y速度、z速度、湍流動(dòng)能k和比耗散率ω,設(shè)置迭代步數(shù)2 000次,采用標(biāo)準(zhǔn)初始化方法,并從進(jìn)口邊界計(jì)算初始值,通過判斷殘差值小于1×10-6和凈質(zhì)量百分比小于0.1%證明計(jì)算收斂[12],殘差曲線如圖5所示。

圖5 殘差
針對(duì)風(fēng)場(chǎng)處于大氣邊界層底部,獲取近地面風(fēng)能資源,受到非定常來流復(fù)雜的影響問題,采用三維數(shù)值仿真還原風(fēng)場(chǎng)模型,剖取流場(chǎng)典型截面的速度數(shù)值云圖,分析研究風(fēng)力機(jī)來流下氣動(dòng)特性。
如圖6所示,為風(fēng)力機(jī)z=1 m剖面和風(fēng)輪葉片表面速度分布云圖。從圖6分析得:風(fēng)輪剖面區(qū)域速度變化范圍為0~89.25 m/s。在旋轉(zhuǎn)域交界面輪廓處存在風(fēng)速最小的速度梯度,形成了旋轉(zhuǎn)環(huán)形流場(chǎng)分布,從葉尖擴(kuò)展至遠(yuǎn)離交界面的位置。

圖6 流場(chǎng)速度云圖
圖6(a)顯示在葉片前后緣有明顯的速度擴(kuò)張,一直延伸到葉尖,考慮到葉片結(jié)構(gòu)對(duì)氣流的阻滯作用,葉片前緣梯度范圍較小,跨越葉片表面氣流補(bǔ)充過程中速度增大,即在后緣有較好的梯度過渡。圖6(b)顯示輪轂區(qū)域處于相對(duì)靜止?fàn)顟B(tài),發(fā)展到葉根開始出現(xiàn)速度線性增長(zhǎng)趨勢(shì),并在葉片跨中處逐漸明顯發(fā)展至葉尖。
如圖7所示,為風(fēng)場(chǎng)在x=0 m和y=2.5 m剖面處風(fēng)力機(jī)速度分布云圖。從圖7(a)中可以看出旋轉(zhuǎn)域在來流風(fēng)速下產(chǎn)生阻塞,近風(fēng)輪區(qū)域速度變化比較集中,致使風(fēng)力機(jī)下游區(qū)域出現(xiàn)尾流拖滯,速度變化范圍擴(kuò)大,沿著入流方向風(fēng)力機(jī)遠(yuǎn)風(fēng)區(qū)風(fēng)速急劇降低,最終到達(dá)初始流速。

圖7 x、y剖面速度云圖
圖7(b)風(fēng)力機(jī)周圍流場(chǎng)誘導(dǎo)區(qū)風(fēng)速開始下降,氣流擾動(dòng)的影響范圍越來越大,壓力梯度隨風(fēng)速降低而增大,至近尾流區(qū)由于葉片幾何外形、塔筒和機(jī)艙的分布位置產(chǎn)生速度虧損,這是由于塔筒效應(yīng)使得其附近流場(chǎng)發(fā)生變化,引起塔架后方風(fēng)速降低。
圖8為計(jì)算域z=1 m剖面處風(fēng)力機(jī)旋轉(zhuǎn)域壓力分布云圖和葉片輪轂表面壓力云圖。從圖8(a)中可以看出風(fēng)輪旋轉(zhuǎn)剪切空氣時(shí),順著葉輪旋轉(zhuǎn)方向,葉片下表面開始產(chǎn)生負(fù)壓力,葉片跨中至葉尖壓力梯度開始增大,考慮到沿葉片相對(duì)展向旋轉(zhuǎn)速度增大,葉片上下表面壓力差導(dǎo)致氣流補(bǔ)充交換速度增大。在旋轉(zhuǎn)域內(nèi)存在明顯的壓力變化,旋轉(zhuǎn)域周圍隨著徑向距離的增大氣流擾動(dòng)作用減弱并趨于穩(wěn)定。表明正面入流風(fēng)對(duì)風(fēng)輪結(jié)構(gòu)體周圍壓力影響較小。

圖8 流場(chǎng)壓力分布云圖
圖8(b)中葉片迎風(fēng)面負(fù)壓范圍為-3 160~-217 Pa,葉片跨中部分開始有明顯的負(fù)壓區(qū)產(chǎn)生,同時(shí)在葉片表面有負(fù)壓集中而不連續(xù),因風(fēng)輪受到垂直入流,葉片前緣最先接收到風(fēng)速擾動(dòng),掠過葉片上下表面后逐漸恢復(fù)至穩(wěn)定。另外葉片受到法向合力驅(qū)動(dòng),迎風(fēng)面和背風(fēng)面均勻分布的壓力由于離心力作用不再完全附著于葉片表面,而是通過能量轉(zhuǎn)換最終有所耗散,符合能量守恒定律,也證明了風(fēng)力機(jī)氣動(dòng)特性的合理性。
圖9(a)為流固耦合作用下風(fēng)力機(jī)結(jié)構(gòu)應(yīng)力云圖。分析可得,固定約束下整體變形以懸挑梁式分布,應(yīng)力主要集中在塔筒兩側(cè),從塔筒底部至頂部應(yīng)力逐漸增大,且在與機(jī)艙連接處產(chǎn)生應(yīng)力集中,塔頂最大應(yīng)力為塔底應(yīng)力的8倍。隨著氣流的來臨葉片在定槳距角條件下,葉片前緣應(yīng)力沿翼型弦長(zhǎng)方向順梯度發(fā)展。

圖9 風(fēng)力機(jī)結(jié)構(gòu)響應(yīng)圖
風(fēng)剪切作用主要是改變?nèi)~片表面應(yīng)力梯度,順時(shí)針轉(zhuǎn)矩作用下,迎風(fēng)側(cè)表面應(yīng)力大于背風(fēng)側(cè),背風(fēng)側(cè)最小應(yīng)力為70.058 Pa,沿展向應(yīng)力先減小后增大再減小,在葉尖處突變至0。風(fēng)輪應(yīng)力集中區(qū)位于葉根同輪轂連接區(qū)域,入流風(fēng)正向剪切造成應(yīng)力主要集中區(qū)域在葉片跨中表面,且應(yīng)力峰值為5.923 9 MPa,遠(yuǎn)小于葉片材料設(shè)計(jì)強(qiáng)度。
風(fēng)力機(jī)整機(jī)應(yīng)變?cè)茍D如圖9(b)所示。分析可知,應(yīng)力應(yīng)變?cè)茍D分布情況基本一致,展現(xiàn)出應(yīng)力應(yīng)變呈現(xiàn)線性相關(guān)關(guān)系。應(yīng)變主要集中區(qū)域位于葉片迎風(fēng)側(cè),且應(yīng)變峰值為8.611 6×10-5m/m,說明葉片是風(fēng)機(jī)運(yùn)行過程中維護(hù)的重點(diǎn)。
風(fēng)力機(jī)停機(jī)工況和旋轉(zhuǎn)工況下葉片變形曲線如圖10所示。從停機(jī)位置右上方葉片開始逆時(shí)針編號(hào)葉片1、葉片2、葉片3。對(duì)比分析可得,上風(fēng)區(qū)葉片1和葉片2在不同工況下總位移曲線以正指數(shù)形式增長(zhǎng)。位移曲線起點(diǎn)表示葉片在葉根處位移且均不為零。在停機(jī)工況下,葉尖處變形達(dá)到峰值,上風(fēng)區(qū)葉片沿展向位移變化規(guī)律基本相同,在旋轉(zhuǎn)工況下也有類似的表現(xiàn)。下風(fēng)區(qū)由于重力和塔影效應(yīng)作用,葉片3整體位移沿展向先減小后增大,葉根至跨中的變形并不明顯,甚至接近于零,脫離塔影重疊后位移量突增至正常水平。

圖10 葉片變形曲線圖
旋轉(zhuǎn)工況下,葉片結(jié)構(gòu)響應(yīng)頻率變大,總體位移較停機(jī)時(shí)有較大增加,上風(fēng)區(qū)峰值增量分別達(dá)到2.143 7、0.867 4 mm。說明風(fēng)輪氣動(dòng)力和離心力作用下顯著影響運(yùn)行過程中葉片變形量。同樣,下風(fēng)區(qū)葉片3位移先減小后增大,與停機(jī)時(shí)不同的是位移突變出現(xiàn)在后半跨,而后恢復(fù)至初始值。當(dāng)葉片與塔筒位置重疊時(shí),氣流擾動(dòng)相互干涉導(dǎo)致葉片變形波動(dòng),表明風(fēng)輪方位角對(duì)結(jié)構(gòu)動(dòng)態(tài)響應(yīng)有較大影響。
為了研究均勻入流條件下不同風(fēng)速風(fēng)力機(jī)輸出功率的變化,設(shè)置入流工況為5、8、10、12、15、18、20 m/s。在FLUENT中分別模擬計(jì)算上述來流工況,繪制風(fēng)輪扭矩收斂曲線,通過數(shù)據(jù)處理獲得不同工況下風(fēng)輪扭矩值M,利用下列公式得出風(fēng)輪輸出功率P,進(jìn)而分析輸出數(shù)據(jù)。
P=Mω
(9)
式(9)中:ω為角速度。
風(fēng)輪在不同來流工況下風(fēng)力機(jī)輸出功率變化曲線如圖11所示。從圖11可以總結(jié)出,風(fēng)場(chǎng)不同均勻入流下風(fēng)輪輸出功率呈現(xiàn)近似線性變化。入流風(fēng)速在5~20 m/s增加過程中,風(fēng)輪輸出功率由68.33 kW增大到84.33 kW,在輸出功率不斷增長(zhǎng)的過程中,考慮到風(fēng)輪扭矩是隨迎風(fēng)速度的增大而增大的,風(fēng)力機(jī)葉片在較大的入流風(fēng)速下所捕獲的氣動(dòng)壓力也更多,由于葉片的結(jié)構(gòu)特征,在流場(chǎng)作用下葉片表面受載使葉片形態(tài)參數(shù)發(fā)生變化,風(fēng)輪升阻比增大,致使葉片的性能得到更好的體現(xiàn),從而輸出功率增大。

圖11 不同來流下輸出功率變化圖
利用有限元軟件ANSYS,研究了風(fēng)力機(jī)在流場(chǎng)作用下流固耦合氣動(dòng)特性,得出以下結(jié)論。
(1)針對(duì)風(fēng)力機(jī)周圍氣動(dòng)環(huán)境的復(fù)雜性,建立并驗(yàn)證了風(fēng)力機(jī)流固耦合風(fēng)場(chǎng)模型。風(fēng)力機(jī)在流場(chǎng)入流風(fēng)和旋轉(zhuǎn)作用下風(fēng)輪周圍的氣流擾動(dòng)最為顯著,沿入流風(fēng)向,在塔影效應(yīng)干擾下,近尾流區(qū)和遠(yuǎn)尾流區(qū)風(fēng)速產(chǎn)生不同程度的耗散。
(2)標(biāo)準(zhǔn)工況下葉片迎風(fēng)側(cè)壓力變化范圍較大,旋轉(zhuǎn)效應(yīng)對(duì)風(fēng)輪壓力變化影響較大,風(fēng)輪周圍壓力主要沿徑向隨氣流擾動(dòng)作用的減弱而變化。
(3)標(biāo)準(zhǔn)工況和停機(jī)工況下上風(fēng)向風(fēng)力機(jī)葉片展向位移變化規(guī)律基本一致,且呈指數(shù)形式增長(zhǎng)。下風(fēng)向由于塔筒氣動(dòng)干擾變形有所波動(dòng)。
(4)隨著均勻入流流速的增大,造成風(fēng)力機(jī)轉(zhuǎn)矩增大,輸出功率也隨之增大。