李松林,朱衛(wèi)軍,孫振業(yè),陶秋晗,曾明伍,3
(1.東方電氣風(fēng)電股份有限公司,四川德陽 618000;2.揚(yáng)州大學(xué) 電氣與能源動(dòng)力工程學(xué)院,江蘇揚(yáng)州 225127;3.西安交通大學(xué) 航天航空學(xué)院,西安 710049)
新一代風(fēng)力機(jī)正趨于大型化發(fā)展,無論是內(nèi)陸低風(fēng)速風(fēng)力機(jī)的推進(jìn)和海上風(fēng)電的技術(shù)創(chuàng)新,風(fēng)力機(jī)葉片尺寸的不斷增長給氣動(dòng)載荷設(shè)計(jì)提出了很高的要求。風(fēng)力機(jī)葉片受到湍流脈動(dòng)載荷的影響、風(fēng)切變、塔架陰影和偏航誤差。這些影響因素組成復(fù)雜非定常氣動(dòng)載荷使葉片產(chǎn)生變形和振動(dòng)。高疲勞載荷和極端負(fù)荷條件的下運(yùn)行的風(fēng)力機(jī)設(shè)計(jì)成本更高,并且維護(hù)成本和系統(tǒng)可靠性更差。為了調(diào)節(jié)功率輸出和降低系統(tǒng)載荷,現(xiàn)代風(fēng)力機(jī)普遍采用了自動(dòng)變槳和自動(dòng)偏航的功能。然而,傳統(tǒng)的變槳和偏航控制系統(tǒng)需要很大的驅(qū)動(dòng)扭矩,面對(duì)大型風(fēng)力機(jī),系統(tǒng)調(diào)節(jié)的靈敏度受到葉片的大質(zhì)量和慣性的限制,系統(tǒng)在實(shí)際上無法快速補(bǔ)償時(shí)間尺度較小的湍流運(yùn)動(dòng)帶來的湍流脈動(dòng)載荷,因此實(shí)際上疲勞載荷很難有效地消除。
對(duì)于未來大型風(fēng)力機(jī)的設(shè)計(jì),智能葉片被認(rèn)為是大幅度降低風(fēng)能成本和提高風(fēng)能利用效率的最優(yōu)前途的方法[1]。尾緣襟翼廣泛應(yīng)用與飛機(jī),在起飛時(shí)增大整體升力,降落時(shí)擺動(dòng)襟翼角度增加飛行阻力。圖1 中顯示采用動(dòng)尾翼(在風(fēng)能領(lǐng)域經(jīng)常被稱為動(dòng)尾翼)的風(fēng)力機(jī)葉輪示意圖。通過改變?nèi)~片后緣的主動(dòng)控制機(jī)構(gòu),使得葉片在主動(dòng)變槳的同時(shí)可以針對(duì)較小的湍流變化相應(yīng)的控制尾翼的角度,從而減少葉片所承受的空氣動(dòng)力載荷??梢灶A(yù)見,未來大型葉片將出現(xiàn)同時(shí)進(jìn)行變槳和變尾緣襟翼的運(yùn)行形式。余畏等[2]在FAST/Aerodyn 代碼中加入柔性尾緣襟翼控制器,成功實(shí)現(xiàn)了葉片在湍流中的降載。變尾緣襟翼包含了動(dòng)尾翼、動(dòng)尾翼與變槳同時(shí)運(yùn)動(dòng)的復(fù)雜模式。采用尾緣襟翼實(shí)現(xiàn)升阻力的實(shí)時(shí)調(diào)節(jié)功能,在失速之前,尾翼向下運(yùn)動(dòng)時(shí)翼型整體升力提高,對(duì)應(yīng)的零升力功角變小,這種效應(yīng)近似于增大了翼型的曲面弧度。與之相反,當(dāng)襟翼向上運(yùn)動(dòng)時(shí),翼型段的氣動(dòng)性能下降。通過實(shí)驗(yàn)觀測(cè),尾緣襟翼在一些工況下還能夠減少翼型的阻力[3]。
圖1 動(dòng)尾翼風(fēng)力機(jī)葉片
可應(yīng)用計(jì)算流體方法(Computational fluid dynamics,CFD)和風(fēng)洞實(shí)驗(yàn)研究俯仰運(yùn)動(dòng)對(duì)翼型的影響。關(guān)于俯翼型仰運(yùn)動(dòng)CFD 仿真方面,國內(nèi)學(xué)著進(jìn)行了充分而精確的研究。謝凱等[4]對(duì)直升機(jī)二維翼型SC1095的俯仰、揮舞、擺振耦合運(yùn)動(dòng)進(jìn)行仿真,采用ICEM軟件生成嵌套網(wǎng)格并求解雷諾平均的NS 方程(RANS)。嵌套網(wǎng)格可以保持翼型邊界層的貼體網(wǎng)格,實(shí)現(xiàn)翼型的俯仰旋轉(zhuǎn)運(yùn)動(dòng),避免了壁面邊界層網(wǎng)格重構(gòu)。徐佩等[5]采用Star-CCM+和重疊網(wǎng)格技術(shù),對(duì)NACA0020翼型的俯仰和升沉耦合運(yùn)動(dòng)進(jìn)行RANS 仿真。張一楠等[6]對(duì)有仿鯨魚鰭前緣的翼型段的俯仰運(yùn)動(dòng)等進(jìn)行了風(fēng)洞試驗(yàn),結(jié)果表明仿生設(shè)計(jì)改善了翼型動(dòng)態(tài)失速特性,降低葉片極限載荷的產(chǎn)生。
關(guān)于動(dòng)尾翼對(duì)風(fēng)力機(jī)翼型的影響,國內(nèi)學(xué)著也進(jìn)行了一些研究。郝文星等[7]將瞬態(tài)風(fēng)速導(dǎo)致的升力與給定升力的差作為控制信號(hào),控制襟翼在Fluent軟件中擺動(dòng),通過Fluent 軟件的網(wǎng)格彈性變形實(shí)現(xiàn)襟翼偏轉(zhuǎn),研究了襟翼最大偏轉(zhuǎn)速度、襟翼動(dòng)作延遲時(shí)間對(duì)控制效果的影響,發(fā)現(xiàn)延遲時(shí)間的影響較大。李傳峰等[8]采用CFD 研究了動(dòng)尾緣襟翼的偏轉(zhuǎn)角頻率對(duì)氣動(dòng)性能的影響,發(fā)現(xiàn)尾緣襟翼降載的能力隨著角頻率的增大而減??;折合頻率的大小可作為非定常氣動(dòng)特性出現(xiàn)的判斷準(zhǔn)則。未來智能風(fēng)力機(jī)葉片會(huì)實(shí)現(xiàn)邊槳距和動(dòng)尾翼的聯(lián)合控制,然而,同時(shí)考慮翼型的俯仰運(yùn)動(dòng)和動(dòng)尾翼偏轉(zhuǎn)的耦合數(shù)值模擬較少。這是由于傳統(tǒng)的貼體網(wǎng)格CFD 方法,在應(yīng)對(duì)翼型整體俯仰外加尾緣局部變形時(shí)存在一些挑戰(zhàn)。為了精確地研究后緣襟翼翼型的物理效應(yīng),本文將采用傳統(tǒng)貼體網(wǎng)格結(jié)合浸入邊界(Immersed boundary,IB)的計(jì)算方法[9-14]。該方法采用標(biāo)準(zhǔn)有限體積法求解翼型主固定部分周圍的繞流,而采用邊界浸入法在曲線網(wǎng)格上模擬可運(yùn)動(dòng)的后緣襟翼。其中,翼型的俯仰運(yùn)動(dòng)通過旋轉(zhuǎn)網(wǎng)格坐標(biāo)系來實(shí)現(xiàn),尾緣襟翼同時(shí)參與俯仰運(yùn)動(dòng)和相對(duì)的柔性擺動(dòng)。
IB 方法專門用于處理任意運(yùn)動(dòng)中復(fù)雜物體的流動(dòng)問題。這個(gè)想法最早是由Peskin[9]應(yīng)用于研究心血管系統(tǒng)中的低雷諾數(shù)流動(dòng)彈性邊界。在高雷諾數(shù)計(jì)算方面,文獻(xiàn)[10]采用S-A 單方程湍流模型和k-ω Wilcox 兩方程模型進(jìn)行了比較計(jì)算。本文介紹雷諾平均法(Reynolds averaged navier stokes,RANS)結(jié)合浸入邊界的數(shù)值表達(dá)方法。NS 方程的動(dòng)量方程和連續(xù)方程為
方程區(qū)別于傳統(tǒng)NS 方程的表達(dá)形式,其中動(dòng)量方程中添加了一個(gè)源項(xiàng)fi。該源項(xiàng)表示了非貼體網(wǎng)格中,浸入邊界部分表達(dá)的虛擬外力作用。如圖2和圖3 所示,全貼體網(wǎng)格計(jì)算中,物體邊界的受力通過求解速度與壓力耦合的方程得到,主要體現(xiàn)為壓強(qiáng)力和黏性力。在非貼體網(wǎng)格中,由于缺失物體壁面黏性邊界條件,如方程不做額外變化,該計(jì)算區(qū)域?qū)⒌韧谄渌耐饬鲌?chǎng)區(qū)域。
圖2 全貼體網(wǎng)格
圖3 貼體網(wǎng)格與浸入邊界網(wǎng)格的混合
式(1)中,將速度的時(shí)間導(dǎo)數(shù)離散化,將其余項(xiàng)移到方程右側(cè)合并,并把右側(cè)合并項(xiàng)表示為RHS,可得
假設(shè)迭代次數(shù)n=1 時(shí)速度場(chǎng)為已知,即初始條件給定,并且浸入邊界的各個(gè)幾何點(diǎn)上的速度uib已知,則臨近浸入邊界網(wǎng)格節(jié)點(diǎn)上的的速度由流場(chǎng)中的速度和邊界運(yùn)動(dòng)速度插值求得,即
因此,式(2)與式(3)建立了速度和外力的迭代關(guān)系。
流體計(jì)算中常常有物體自身做復(fù)雜運(yùn)動(dòng)的情況,理論上可以采用上述浸入邊界的方法求解整個(gè)物體周邊的流動(dòng)場(chǎng),比如當(dāng)前情況下翼型的復(fù)雜運(yùn)動(dòng)。但是全部采用浸入邊界方法將大大增加計(jì)算量。采用旋轉(zhuǎn)坐標(biāo)系結(jié)合尾緣襟翼位置浸入邊界方法將大大提高計(jì)算效率和精度。在旋轉(zhuǎn)坐標(biāo)系下,式(1)成立必須滿足前提條件,即流體質(zhì)點(diǎn)的加速度必須是相對(duì)于慣性坐標(biāo)系下的加速度。在旋轉(zhuǎn)坐標(biāo)系下,式(1)的速度導(dǎo)數(shù)項(xiàng)(即質(zhì)點(diǎn)加速度)必須包含額外的加速度項(xiàng),即
式中:|?|為旋轉(zhuǎn)速度;r為質(zhì)點(diǎn)與旋轉(zhuǎn)中心的距離矢量;V為質(zhì)點(diǎn)運(yùn)動(dòng)的速度矢量。對(duì)于當(dāng)前的二維翼型計(jì)算,展開各個(gè)向量積可得:
式(5)各項(xiàng)相加作為動(dòng)量式(1)的添加項(xiàng)放入右側(cè),此方程將結(jié)合IB 方法求解翼型俯仰運(yùn)動(dòng)與動(dòng)尾緣襟翼的疊加效應(yīng)。
計(jì)算程序采用速度壓力耦合方程,應(yīng)用 SIMPLEC/PISO 方法和多重網(wǎng)格策略。在動(dòng)量方程中首先用預(yù)猜測(cè)的壓力作為預(yù)測(cè)因子求解方程,然后,用連續(xù)方程為約束,得到壓力修正方程。時(shí)間和空間差分均采用2 階精度,其中空間差分采用有限體積法,而湍流擴(kuò)散項(xiàng)采用中心差分格式求解。在校正步驟,采用Rhie Chow 插值技術(shù)是用來抑制速度壓力耦合求解帶來的數(shù)值振蕩問題。在結(jié)構(gòu)化網(wǎng)格中結(jié)合IB 方法,方程添加外力項(xiàng)如式(1)所示,在每個(gè)迭代步長中更新速度場(chǎng),從而更新外力項(xiàng)。其具體的計(jì)算流程如圖4 所示,首先初始化流場(chǎng),應(yīng)用經(jīng)典流體力學(xué)中的邊界條件和初始值,然后僅在第一個(gè)計(jì)算步長讀入IB 外形幾何數(shù)據(jù),該組數(shù)據(jù)通常需要根據(jù)點(diǎn)的密集程度重新插值組合。在每一時(shí)刻,判斷浸入邊界的內(nèi)外計(jì)算節(jié)點(diǎn),如浸入邊界相對(duì)于翼型本地坐標(biāo)的運(yùn)動(dòng)速度為零,則該過程僅需執(zhí)行一次。由于RANS 方法求解k-ω方程時(shí)采用的壁面邊界條件與垂直距離有關(guān),因此在虛擬的IB 邊界上同樣需要求解壁面垂直距離。根據(jù)給定的IB 運(yùn)動(dòng)方程(采用主動(dòng)控制方法的除外),求得該時(shí)刻IB 上每個(gè)點(diǎn)在水平方向、垂直方向的位移和速度,如尾緣襟翼為靜止?fàn)顟B(tài),則位移和速度均為零。判斷下一步尾緣襟翼是否運(yùn)動(dòng),如運(yùn)動(dòng)則返回,重新區(qū)分并標(biāo)識(shí)IB 分割線的內(nèi)外計(jì)算網(wǎng)格節(jié)點(diǎn)。在IB 附近的網(wǎng)格點(diǎn)上進(jìn)行速度的插值運(yùn)算,采用IB 上的運(yùn)動(dòng)速度與流場(chǎng)內(nèi)部相鄰點(diǎn)的速度進(jìn)行插值。根據(jù)該速度值更新外力項(xiàng),將此外力合并至動(dòng)量方程的源項(xiàng)中,進(jìn)一步更新流場(chǎng)。下一個(gè)步長,如尾緣襟翼持續(xù)改變位置,則重復(fù)進(jìn)行迭代計(jì)算,直至設(shè)定的總迭代步數(shù),最后積分求得升阻力。
圖4 結(jié)合浸入邊界方法的NS 方程求解流程圖
當(dāng)前的模擬中,尾緣襟翼部分由函數(shù)表示,該函數(shù)能夠較好地描述翼型尾緣襟翼的柔性變化。動(dòng)尾翼部分各個(gè)坐標(biāo)點(diǎn)沿著x軸和y軸的位移變化式為:
如圖5 所示,a、b兩點(diǎn)始終固定,其余各點(diǎn)按照式(6)和式(7)的定義變化位置。
圖5 尾緣襟翼彈性運(yùn)動(dòng)示意圖
當(dāng)?shù)介L為n+1 時(shí),新位置坐標(biāo)等于前一步的位置坐標(biāo)加上單位時(shí)間 ?t內(nèi)產(chǎn)生的新位移,即和中幾何邊界的運(yùn)動(dòng)速度為時(shí)間的函數(shù),即考慮存在任意運(yùn)動(dòng)形式的情況。方程中的冪指數(shù)p描述了尾緣襟翼的柔度,當(dāng)p=1 時(shí),尾緣襟翼做剛性運(yùn)動(dòng),當(dāng)采用冪指數(shù)為p=2 時(shí)的運(yùn)動(dòng)形態(tài)如圖5 所示。
實(shí)踐中,由于貼體網(wǎng)格的計(jì)算效率更高,靜態(tài)翼型無需采用IB 方法進(jìn)行模擬。然而作為IB 方法的驗(yàn)證,貼體網(wǎng)格算法可以作為一種參考方法進(jìn)行對(duì)比。圖6 顯示了NACA0012 翼型的表面壓力系數(shù)分布,雷諾數(shù)為5×105,來流攻角 α=4?。由圖6 可知,在靜止?fàn)顟B(tài)下,采用貼體網(wǎng)格和浸入邊界方法計(jì)算所得到的壓力分布結(jié)果幾乎完全對(duì)應(yīng)。采用IB 方法在尾緣部分需要的網(wǎng)格要多于貼體網(wǎng)格,若采用IB 方法模擬整體翼型,其計(jì)算效率將低于傳統(tǒng)貼體網(wǎng)格算法。
圖6 尾翼靜止時(shí)貼體網(wǎng)格和混合網(wǎng)格對(duì)比
對(duì)于流動(dòng)中包含復(fù)雜幾何外形或物體做復(fù)雜運(yùn)動(dòng)的情況,IB 方法展現(xiàn)了其強(qiáng)大的靈活性。對(duì)于運(yùn)動(dòng)的尾緣襟翼流場(chǎng)計(jì)算,傳統(tǒng)的貼體網(wǎng)格算法必須在每個(gè)時(shí)間步長上根據(jù)翼型新的形體重新生成新的網(wǎng)格,如此大大增加了計(jì)算的時(shí)間。采用IB 方法將保持固定的網(wǎng)格,尾緣襟翼部分則在網(wǎng)格中做相對(duì)運(yùn)動(dòng),其運(yùn)動(dòng)的形式可以由任意方程提前給定。比如,當(dāng)前的尾緣襟翼的運(yùn)動(dòng)采用正弦函數(shù)描述,即
式中:δo為初始角度;A為最大位移;k1為折算頻率(Reduced frequency);φ 為相位角。結(jié)合式(6)和式(7)關(guān)于位移的函數(shù)描述可以得到尾緣襟翼每個(gè)點(diǎn)的運(yùn)動(dòng)速度:
式中:r(t)為第i個(gè)邊界點(diǎn)距離轉(zhuǎn)動(dòng)中心的距離,在每個(gè)時(shí)刻,尾緣襟翼產(chǎn)生的線速度分解為水平方向和垂直方向的速度分量
圖7 中顯示了某一時(shí)刻N(yùn)ACA0012 翼型尾緣襟翼運(yùn)動(dòng)時(shí)的水平速度分量云圖和速度流線圖,圖8中可見升力的變化同樣近似于正弦,這也與主觀預(yù)期一致。當(dāng)前計(jì)算采用入流功角為4°,雷諾數(shù)為1.63×106,翼型尾緣襟翼部分大小為整體弦長的22.6%。由圖8 的升力變化結(jié)果可以看到,翼型升力的變化主要取決于幾個(gè)條件:尾緣襟翼相對(duì)于整體翼型的大小、擺動(dòng)的幅度和擺動(dòng)的頻率。翼型襟翼的運(yùn)動(dòng)方式可以提前給定也可以根據(jù)升力的大小自動(dòng)調(diào)節(jié),后者涉及到升力的主動(dòng)控制,即給定控制算法在每個(gè)時(shí)間步長中判斷尾緣襟翼下一步的位移量,從而使升力不斷接近目標(biāo)量。當(dāng)前關(guān)于翼型與運(yùn)行條件的選取是為了與風(fēng)洞實(shí)驗(yàn)進(jìn)行真實(shí)對(duì)比。下文將模擬翼型的復(fù)雜運(yùn)動(dòng)并與風(fēng)洞實(shí)驗(yàn)進(jìn)行比較。
圖7 尾翼運(yùn)動(dòng)時(shí)水平速度分量云圖
圖8 翼型升力曲線
為了驗(yàn)證動(dòng)尾翼計(jì)算的精度,本文將對(duì)比Krzysiak 等[15]的風(fēng)洞實(shí)驗(yàn)結(jié)果。在文獻(xiàn)[15]中,作者對(duì)NACA0012 翼型進(jìn)行了風(fēng)洞氣動(dòng)實(shí)驗(yàn),并與半經(jīng)驗(yàn)的理論模型進(jìn)行了對(duì)比,由于模型建立在無粘流動(dòng)基礎(chǔ)之上,在許多工況的對(duì)比中發(fā)現(xiàn)理論模型與實(shí)際相差較大,由此本文展開數(shù)值計(jì)算方法模擬翼型復(fù)雜的耦合運(yùn)動(dòng)。如圖9 所示,該翼型的俯仰旋轉(zhuǎn)中心位于距離翼型前緣位置0.35C處,尾緣襟翼的長度為0.226C,其旋轉(zhuǎn)中心位于0.2C處。俯仰與尾緣襟翼運(yùn)動(dòng)函數(shù)組合如表1 所示。
表1 俯仰與尾緣襟翼運(yùn)動(dòng)函數(shù)組合
圖9 風(fēng)洞實(shí)驗(yàn)中NACA 0012 翼型和可動(dòng)襟翼[15]
入流攻角采用正弦周期性函數(shù)表示為
組合1~組合3 中,入流攻角的初始位置為 4?,俯仰的最大振幅為 6?;尾緣襟翼初始位置為 0?,擺動(dòng)最大幅度為5.4?。實(shí)驗(yàn)中分別采用了3 組相位差:148?,298?,357?。圖10a)、圖11a)和圖12a)中分別顯示了計(jì)算模擬采用的動(dòng)態(tài)函數(shù),其中橫坐標(biāo)為攻角,縱坐標(biāo)為尾緣襟翼角。由于相位不同,3 組情況下運(yùn)動(dòng)形態(tài)具有顯著差異。圖10b)、圖11b)和圖12b)中分別比較了單獨(dú)俯仰運(yùn)動(dòng)和組合運(yùn)動(dòng)的動(dòng)態(tài)升力。單獨(dú)俯仰運(yùn)動(dòng)時(shí),翼型的升力呈經(jīng)典的“O”型閉合循環(huán),數(shù)值模擬與實(shí)驗(yàn)整體擬合較好。當(dāng)尾緣襟翼一起運(yùn)動(dòng)時(shí),升力曲線在初始攻角附近形成交叉,呈“∞”型閉合循環(huán),結(jié)果整體吻合。需要指出的是,由于實(shí)驗(yàn)機(jī)構(gòu)的振動(dòng)和伺服機(jī)構(gòu)的缺陷,實(shí)際的俯仰和尾緣運(yùn)動(dòng)趨近于正弦函數(shù),但是仍有明顯的差異[15](圖13 和圖14),這些也是導(dǎo)致計(jì)算與實(shí)驗(yàn)存在一定誤差的原因之一。其次,根據(jù)計(jì)算流體的模擬經(jīng)驗(yàn),對(duì)于雷諾平均方法而言,模擬翼型的大攻角和失速區(qū)間結(jié)果精確度通常較低,不同湍流模型之間往往也差別較大。適當(dāng)?shù)慕档透┭龊臀簿壗笠磉\(yùn)動(dòng)的幅值可以改善計(jì)算模擬的準(zhǔn)確性。組合4和組合5 中變化了運(yùn)動(dòng)函數(shù)組合,其中兩種運(yùn)動(dòng)的幅值均有所降低。圖13 中,動(dòng)態(tài)升力的比較與實(shí)驗(yàn)基本吻合。由于伺服機(jī)構(gòu)的誤差,組合4 中實(shí)驗(yàn)所測(cè)的攻角范圍超出了最小極限值?0.5°較多。圖14中,動(dòng)態(tài)升力的變化無論從趨勢(shì)上還是幅值上都實(shí)現(xiàn)與實(shí)驗(yàn)較好的吻合。
圖10 組合1 計(jì)算結(jié)果
圖11 組合2 計(jì)算結(jié)果
圖12 組合3 計(jì)算結(jié)果
圖13 組合4 計(jì)算結(jié)果
圖14 組合5 計(jì)算結(jié)果
本文針對(duì)風(fēng)力機(jī)葉片發(fā)展大型化和智能化的趨勢(shì),從翼型變槳控制結(jié)合動(dòng)尾翼控制的策略出發(fā),探索了一種基于浸入邊界方法的實(shí)踐應(yīng)用。為了與風(fēng)洞試驗(yàn)進(jìn)行比較,本文的研究中,變槳運(yùn)動(dòng)和尾緣襟翼運(yùn)動(dòng)均采用正弦函數(shù)表達(dá)方式。葉片的實(shí)際的運(yùn)動(dòng)形式應(yīng)完全由給定的升力、阻力等目標(biāo)量進(jìn)行反饋和輸出控制。本文采用了5 種不同組合的槳距角變化和尾緣的變化形式,結(jié)果顯示出兩者的組合形式帶來的翼型升力變化是顯著的,同時(shí)結(jié)果與實(shí)驗(yàn)數(shù)據(jù)相符。