風(fēng)力機(jī)氣動特性的浸入邊界法模擬1)

李燕玲 胡進(jìn) 周錕

(武漢科技大學(xué)省部共建耐火材料與冶金國家重點(diǎn)實(shí)驗(yàn)室，武漢430081)

據(jù)統(tǒng)計(jì)[1]，全球可再生能源(風(fēng)能，太陽能，以及生物質(zhì)能等)于2019年占據(jù)全部一次能源的5%，并且實(shí)現(xiàn)了12.1%年度增長。其中，風(fēng)能占全部可再生能源的51%，并且以30%的年度速度增長。而中國是可再生能源增長速度最快的國家。風(fēng)力機(jī)是將自然界的風(fēng)能轉(zhuǎn)換成機(jī)械能或電能的裝置。其重要工作部件是葉片，決定了風(fēng)能轉(zhuǎn)換的效率。因此，提高葉片的氣動性能是研發(fā)風(fēng)力機(jī)的基礎(chǔ)與關(guān)鍵，也是難點(diǎn)之一。

風(fēng)力機(jī)也隸屬于葉輪機(jī)械。葉輪機(jī)葉片氣動性能研究主要包括理論研究、數(shù)值模擬和風(fēng)洞實(shí)驗(yàn)。理論研究發(fā)展較早[2-3]，如常用的Betz理論、葉素理論、動量理論以及葉素?動量理論等發(fā)展相對成熟。然而理論研究的計(jì)算模型、求解條件都比較復(fù)雜，并且實(shí)際工作環(huán)境中存在許多不可預(yù)測因素，通常需要對模型簡化以得到相關(guān)解[4-5]。風(fēng)洞試驗(yàn)雖然能夠準(zhǔn)確地控制實(shí)驗(yàn)條件，受氣候條件影響小，但需要花費(fèi)大量的時(shí)間和資金。隨著研究的深入和計(jì)算機(jī)技術(shù)的發(fā)展，越來越多的研究人員傾向用數(shù)值模擬方法對旋轉(zhuǎn)葉片流場特性進(jìn)行研究分析。而當(dāng)前的計(jì)算流體力學(xué)(computational fluid dynamics,CFD)方法，可分為以下幾類。(1)多重旋轉(zhuǎn)坐標(biāo)系方法[6-7]：在葉片周圍使用旋轉(zhuǎn)坐標(biāo)系，設(shè)置旋轉(zhuǎn)域和葉片轉(zhuǎn)速相同，模擬過程中旋轉(zhuǎn)域和葉片之間不發(fā)生網(wǎng)格的相對運(yùn)動，將非定常問題轉(zhuǎn)化為定常問題，可大幅節(jié)省模擬計(jì)算的時(shí)間，但計(jì)算精度不是非常高。(2)滑移網(wǎng)格法[8]：需要劃分包含風(fēng)機(jī)葉輪的旋轉(zhuǎn)域和剩余的固定域，在旋轉(zhuǎn)的過程中兩域在交界面處相對滑動，但此方法需要對交界面進(jìn)行特殊設(shè)置。(3)重疊網(wǎng)格法[9-10]：將復(fù)雜的流動區(qū)域分成幾何邊界比較簡單的子區(qū)域，各子區(qū)域中的計(jì)算網(wǎng)格獨(dú)立生成，彼此存在著重疊關(guān)系，但在計(jì)算域運(yùn)動的過程中，需要不斷檢測網(wǎng)格的重合區(qū)域，流場信息通過插值在重疊區(qū)邊界進(jìn)行匹配和耦合。(4)動網(wǎng)格法[11]：動網(wǎng)格法被運(yùn)用在瞬態(tài)模擬中，風(fēng)機(jī)葉片結(jié)構(gòu)在旋轉(zhuǎn)過程中不再是一個(gè)剛體，而會由于受力導(dǎo)致變形，因此每一時(shí)間步中均需對計(jì)算域中的網(wǎng)格進(jìn)行變形與重構(gòu)，確保網(wǎng)格質(zhì)量滿足模擬計(jì)算要求。

基于浸入邊界法[12]，本文發(fā)展了一種簡單的建模，網(wǎng)格離散，與數(shù)值模擬統(tǒng)一方法框架，來模擬風(fēng)力機(jī)的氣動問題。在此框架下，整個(gè)風(fēng)力機(jī)(包括葉片，輪轂，機(jī)艙，以及塔架)與周圍的空氣流場所占據(jù)的空間都用統(tǒng)一的均勻網(wǎng)格來離散劃分。而風(fēng)力機(jī)與空氣的界面，采用散點(diǎn)來離散。在這些離散點(diǎn)上，通過對空氣流場施加合適的作用力來模擬流固耦合作用。對于葉片構(gòu)型，本文提出一種簡單的同倫變形，來光滑融合不同的二維翼型從而生成更加復(fù)雜多樣的三維葉片。同時(shí)，采用仿射變換，來處理葉片的扭轉(zhuǎn)與漸縮問題。對于所有的風(fēng)力機(jī)部件，本文都給出很簡單的離散點(diǎn)生成方法。相較于基于貼體網(wǎng)格的各類模擬方法而言，本文采用的方法在模型離散時(shí)更加強(qiáng)健，可以非常容易地處理模型各部件之間的狹縫，區(qū)域重疊等常見問題，且數(shù)值穩(wěn)定性強(qiáng)。此方法框架的另一個(gè)顯著特點(diǎn)，是將建模，網(wǎng)格離散，數(shù)值模擬融為一體，這使得有可能在此框架下發(fā)展設(shè)計(jì)選型(關(guān)于翼型，尺寸，扭轉(zhuǎn)角等參數(shù))自動優(yōu)化算法--而這是本課題組的長期研究方向之一。

為了檢驗(yàn)方法與程序的正確性，本文模擬了二維典型翼型在不同攻角下的升阻力問題，并且與高精度的譜元方法模擬結(jié)果進(jìn)行了對比。針對浸入邊界法中阻力過估明顯的問題，提出了一種簡單有效的理查森線性外推誤差修正方法。同時(shí)，模擬研究了拱曲度以及厚度對二維翼型升阻力的影響。隨后，模擬研究了單風(fēng)力機(jī)(包含塔架)在不同尖速比下的功率系數(shù)，并對塔架與葉片間的相互氣動作用進(jìn)行了初步分析。最后，模擬研究了雙風(fēng)力機(jī)在風(fēng)場中不同前后間隔距離下的氣動干涉問題。

1 研究方法

1.1 流場模擬方法——浸入邊界法

浸入邊界法是一種很常見的處理流固耦合問題的方法。在此方法中，包含固體與流體的全部空間，都采用歐拉坐標(biāo)下的納維?斯托克斯方程來描述

式中，u，p，ρf和v分別是流體的速度、壓力、密度和運(yùn)動黏度，?為哈密頓算子，fIBM表示固體對流體的作用力。

另一方面，對于流體中固體的運(yùn)動，則是采用拉格朗日坐標(biāo)系下的牛頓運(yùn)動方程來描述。比如固體的平動，其方程為

式中，ρp與up分別表示固體的密度與平動速度，而FIBM則表示作用于固體上的作用力。當(dāng)固體進(jìn)行旋轉(zhuǎn)時(shí)，還要建構(gòu)類似的轉(zhuǎn)動方程。

理想狀態(tài)下，固體對流體的作用力?FIBM與流體對固體的作用力FIBM應(yīng)滿足牛頓第三定律，即大小相等，方向相反。因?yàn)楣腆w對流體的作用力?FIBM與流體對固體的作用力FIBM是分別定義在兩種不同坐標(biāo)下(即歐拉坐標(biāo)與拉格朗日坐標(biāo))，一般情況下，兩者空間位置并不重合，因而需要采用不同坐標(biāo)間的相互插值來進(jìn)行交互。最常用的則是采用離散δ函數(shù)。

浸入邊界法中最后一步，是構(gòu)建流固間作用力模型。對于剛體，最流行的方式是直接作用力方法[12]

式中，?U是在固體邊界位置處流體的插值速度，?t是時(shí)間步長，即流固間的相互作用力等于流體與固體在界面處速度的差異除以時(shí)間步長。換言之，直接作用力是通過施加一個(gè)正比于流固速度差的回復(fù)力，而使得流固間的不可滑移邊界條件得到滿足。

本研究中采用的程序，是基于經(jīng)典譜方法槽道流動求解器，并在其中引入額外的作用力項(xiàng)來模擬流固耦合作用。關(guān)于浸入邊界法的技術(shù)細(xì)節(jié)以及程序?qū)崿F(xiàn)，可參照文獻(xiàn)[12-14]。

1.2 風(fēng)力機(jī)葉片建模及網(wǎng)格離散

風(fēng)力機(jī)葉片是收集風(fēng)能的核心部件，其作用原理與機(jī)翼非常類似?？諝饬鬟^葉片時(shí)，在葉片上同時(shí)產(chǎn)生升力與阻力。風(fēng)力機(jī)葉片多采用高升阻比翼型，亦即更大比例地利用升力作用來驅(qū)動葉片轉(zhuǎn)動，從而獲得更高的風(fēng)能捕獲效率。出于對效率和結(jié)構(gòu)強(qiáng)度等多方面的考慮，葉片形狀非常多樣。一般來說，為了提高結(jié)構(gòu)強(qiáng)度，在葉片根部采用更厚重的翼型；而在尖部采用更輕薄的翼型。一方面，在改變翼型的相對厚度的同時(shí)，朝著尖端方向葉片的弦長也不斷減小，從而使得整個(gè)葉片呈現(xiàn)尖錐狀。另一方面，沿著翼展方向葉片也呈現(xiàn)不同角度的扭轉(zhuǎn)。這是因?yàn)槿~片的旋轉(zhuǎn)線速度沿翼展方向不斷增大，從而使得空氣來流與葉片的相對速度沿翼展方向也不斷改變。為了保證相對來流的攻角在整個(gè)翼展方向都處于理想范圍，需要沿翼展方向適當(dāng)扭轉(zhuǎn)葉片。

此處，介紹一種簡單地生成光滑漸變?nèi)~片的方法。假設(shè)在葉片最大弦長處(靠近根部)采用翼型B，而在最小弦長處(葉片尖端)采用翼型C。翼型B和C可以從通用翼型庫中選擇(國家可再生能源實(shí)驗(yàn)室官方網(wǎng)站上提供了幾十種二維翼型供選擇)。每種二維翼型，都可以從其前緣到后緣沿著弧長進(jìn)行參數(shù)化，即以前緣為零點(diǎn)到前緣的表面弧長為坐標(biāo)。這些坐標(biāo)，都采用總弧長作為歸一化參數(shù)，從而使得坐標(biāo)從0到1之間變化(本文將此坐標(biāo)記為l)。在翼型B和C之間，采用同倫變形來生成中間的翼型，使得翼型從B光滑變化到C。此同倫變化數(shù)學(xué)上可表示為

式中，D表示位于B與C之間s處的翼型，而s是沿翼展方向距離B的歸一化距離(即s=0表示B的位置，而s=1則是C的位置)。上面的同倫變形，給出了在任一位置s，弧長坐標(biāo)l處的翼型坐標(biāo)，從而唯一確定了B與C之間所有的翼型。這些翼型從B沿翼展方向光滑變化到C。這種構(gòu)建復(fù)雜翼型的方法非常簡單，且可直接推廣到更復(fù)雜的情形，即在葉片多處位置指定不同的二維翼型，而在任意兩段翼型之間進(jìn)行同倫變形來生成光滑過渡翼型。

此處介紹的同倫變形翼型構(gòu)建方法，還可以非常便捷地生成葉片的表面離散網(wǎng)格。作者關(guān)于浸入邊界法的前期研究表明[14]，當(dāng)固體的表面離散網(wǎng)格點(diǎn)呈現(xiàn)均勻分布，與整個(gè)空間的均勻歐拉網(wǎng)格相仿時(shí)，將具有更好的數(shù)值精度。所以本文研究模擬的葉片(包括后面討論的風(fēng)力機(jī)其他組件)，都采用盡可能均勻的表面離散點(diǎn)。生成葉片表面離散點(diǎn)時(shí)，首先給定離散點(diǎn)之間的需要間隔dx，此間隔應(yīng)該與歐拉網(wǎng)格的尺寸相仿。然后沿著翼展方向，以此間隔dx，垂直于翼展方向截取葉片形成一系列的截?cái)嗝妗τ诿總€(gè)截?cái)嗝?，先?jì)算出其總周長，再以總周長除以間隔dx來確定所需的離散點(diǎn)數(shù)目，最后沿著周向等間距生成此數(shù)目的離散點(diǎn)。

對于真實(shí)葉片，還有兩個(gè)參數(shù)需要考慮：錐度與扭轉(zhuǎn)角。這兩個(gè)參數(shù)都可以采用仿射變換來統(tǒng)一處理。仿射是包含平移、旋轉(zhuǎn)、鏡像和放縮等一系列幾何相似變形的組合。在數(shù)學(xué)上可以通過將原始幾何坐標(biāo)左乘一個(gè)變換矩陣來實(shí)現(xiàn)。對于葉片而言，所需要的變換是在沿著翼展方向的每個(gè)截?cái)嗝嫔?，對每個(gè)翼型周圍的離散點(diǎn)進(jìn)行適當(dāng)?shù)姆趴s和旋轉(zhuǎn)，從而獲得所需的錐度和扭轉(zhuǎn)角。這種仿射變換表達(dá)式為(不失一般性，假設(shè)二維離散點(diǎn)處于x?y平面上)

式中，S為由離散點(diǎn)x，y坐標(biāo)以及數(shù)字1填充而成的3×N矩陣(N是離散點(diǎn)數(shù)目)，S′則是對應(yīng)的變形之后的離散點(diǎn)坐標(biāo)矩陣。方程右邊三個(gè)矩陣分別表示離散點(diǎn)需要的平移量xt和yt，旋轉(zhuǎn)角度θ，以及放縮系數(shù)W。平移量直接由葉片最終的空間位置確定，旋轉(zhuǎn)角度就是所需的扭轉(zhuǎn)角，而放縮系數(shù)則由錐度確定。一般來說，葉片的弦長沿翼展方向大多呈線性變化(固定錐度)，此時(shí)放縮系數(shù)沿整個(gè)翼型方向都是固定值；但旋轉(zhuǎn)角θ沿翼展方向通常是變化的。值得特別說明的是，此處采用了擴(kuò)展的變換矩陣來統(tǒng)一處理平衡、旋轉(zhuǎn)與放縮問題。在最后結(jié)果S′矩陣中的第三行數(shù)據(jù)是冗余的，只有第一與第二行數(shù)據(jù)分別包含所需的變換后離散點(diǎn)x與y坐標(biāo)。

1.3 風(fēng)力機(jī)輪轂，機(jī)艙，塔架建模及網(wǎng)格離散

風(fēng)力機(jī)輪轂，機(jī)艙以及塔架等的氣動力直接作用遠(yuǎn)小于葉片。在本文中，輪轂、機(jī)艙和塔架模型分別采用簡單的半橢球體、長方體和細(xì)長圓臺。對于半橢球及圓臺這兩種中心對稱的幾何體，其離散方法與前面的葉片沿翼展方向的離散完全相同，而對于長方體表面則可直接采用正交均勻網(wǎng)格離散。

2 風(fēng)力機(jī)幾何模型及離散點(diǎn)示例

本文將風(fēng)力機(jī)置于經(jīng)典槽道流動中來模擬風(fēng)力機(jī)的氣動力問題。槽道橫向、流向以及垂直方向分別記為x，y和z。槽道的底部設(shè)置為不可滑移的零速度邊界條件來模擬地面。頂部采用給定的速度邊界來模擬穩(wěn)定的空氣來流。在其他方向都采用周期性邊界(因?yàn)椴捎玫母道锶~譜方法)。對于葉片，選擇NREL S825和S805A兩種翼型(分別對應(yīng)圖1中B和C)，按照前文介紹的同倫變形方法來生成(其三維投影見圖2)。圖1中O對應(yīng)于30%弦長位置，是翼型扭轉(zhuǎn)中心，同時(shí)也是葉片的裝配中心。兩種翼型的方向?qū)?yīng)于各自的扭轉(zhuǎn)角(葉片旋轉(zhuǎn)平面為x?z)。翼型S825具有更大的相對厚度，使其位于葉片近根部具有最大弦長的位置(圖2中B位置)。而相對厚度較小的翼型S805A則位于葉片尖端(圖2中C位置)。翼型S805A的弦長設(shè)置為最大弦長的一半，使得整個(gè)葉片形成錐狀。沿著翼展方向葉片不斷扭轉(zhuǎn)，其扭轉(zhuǎn)角(翼弦與葉片轉(zhuǎn)動平面的夾角)采用文獻(xiàn)[15]中給出的非線性分布。在葉片根部扭轉(zhuǎn)角最大為30°；而在尖端扭轉(zhuǎn)角較小為?1.9°。除了B至C段的漸變翼型，葉片模型的圓形的裝配底座(圖2中A位置)以及從A到B的過渡段的離散點(diǎn)同樣采用同倫變形方法來生成。整個(gè)葉片，從轉(zhuǎn)軸到翼尖為5(本文所有長度都采用無量綱化數(shù)值，以葉片的最大弦長為標(biāo)準(zhǔn)長度)。以旋轉(zhuǎn)軸為零點(diǎn)，裝配底座沿翼展方向位于0.5處，最大弦長位置為1.25。圖2給出了整個(gè)葉片的真實(shí)比例結(jié)構(gòu)圖和離散點(diǎn)間隔dx=0.1條件下的表面散點(diǎn)圖。在此條件下單個(gè)葉片大約有860個(gè)表面離散點(diǎn)。

圖1 葉片最大弦長及最小弦長處的二維翼型

圖2 葉片模型尺寸以及表面離散點(diǎn)

整個(gè)風(fēng)機(jī)(輪轂，機(jī)艙以及塔架)的裝配模型以及離散點(diǎn)如圖3所示(全部風(fēng)力機(jī)表面離散點(diǎn)約為6100個(gè))。不同于貼體網(wǎng)格類型的CFD方法，浸入邊界法對幾何模型的表面離散要求非常低，并不要求相臨離散單元之間滿足任何拓?fù)湎嗲㈥P(guān)系。換言之，不同部件離散單元之間重合和穿透等，都不會給數(shù)值方法本身造成任何困難。對于風(fēng)力機(jī)這樣復(fù)雜的結(jié)構(gòu)，可以通過本文介紹的方法，在不依賴任何復(fù)雜的幾何建模和網(wǎng)格生成軟件的前提下，進(jìn)行成功模擬。

圖3 風(fēng)力機(jī)整機(jī)裝配模型以及表面離散點(diǎn)

3 模擬結(jié)果

3.1 二維機(jī)翼繞流

為了驗(yàn)證數(shù)值模擬結(jié)果的精度，本文模擬了不同網(wǎng)格精度、不同來流攻角α條件下二維翼型NREL S805A的繞流問題。分析攻角分別為0°，5°，10°，15°以及20°，雷諾數(shù)Re(雷諾數(shù)定義為來流速度乘以弦長并除以運(yùn)動黏性)為10和100的升阻力結(jié)果。值得指出的是，此處模擬的雷諾數(shù)非常低，遠(yuǎn)低于風(fēng)力機(jī)工程實(shí)際應(yīng)用中的雷諾數(shù)。因?yàn)楸狙芯克薪Y(jié)果都是采用直接數(shù)值模擬，而眾所周知，直接數(shù)值模擬方法的計(jì)算量近似正比于雷諾數(shù)的三次方[16]，使得無法真正模擬實(shí)際工況。此算例主要是用來便捷地檢驗(yàn)算法的網(wǎng)格依賴性與精度。在后文模擬的三維風(fēng)力機(jī)，雷諾數(shù)更加接近工程實(shí)際，并將雷諾數(shù)對結(jié)果的影響進(jìn)行了探討。

作為對比標(biāo)準(zhǔn)，本文還采用了基于貼體網(wǎng)格的高階精度譜元方法Nek5000[17]，模擬了相同工況下翼型的升阻力。圖4給出了Re=10條件下的阻力系數(shù)(Cd=2Fd/(ρfu2∞c)，其中Fd是阻力，u∞是來流的速度，c是弦長)與升力系數(shù)(Cl=2Fl/(ρfu2∞c)，其中Fl是升力)隨攻角的變化。在每一幅子圖中，都包含5組數(shù)據(jù)，分別對應(yīng)于4種不同的網(wǎng)格精度與譜元方法的結(jié)果。這4種網(wǎng)格精度分別為dx=0.05，0.025，0.012 5，0.006 25。對于升力系數(shù)，浸入邊界法與譜元方法升力結(jié)果吻合很好。這說明浸入邊界法對于升力(也就是翼型上下表面的壓差)模擬非常準(zhǔn)確。但是對于阻力，則與譜元方法結(jié)果相差明顯，這種差距明顯依賴于網(wǎng)格大小。

圖4 不同攻角α，不同網(wǎng)格尺寸dx條件下二維翼型S805A阻力與升力與譜元方法Nek5000結(jié)果對比

為了進(jìn)一步研究阻力與網(wǎng)格尺寸之間的關(guān)系，圖5給出了不同網(wǎng)格大小下浸入邊界法阻力與譜元方法結(jié)果的相對誤差表示譜元方法阻力)。圖中還包括Re=10條件下的結(jié)果。從圖中可以清楚看出，不論是在何種攻角或者雷諾數(shù)條件下，相對誤差都幾乎正比于網(wǎng)格尺寸dx。即浸入邊界法對于阻力的模擬，相對網(wǎng)格尺寸具有較嚴(yán)格的一階精度。這與作者之前的平板黏性阻力理論研究結(jié)果[14]完全相符，并進(jìn)一步表明之前的理論結(jié)果對于翼型這種復(fù)雜的構(gòu)型也適用。一方面，因?yàn)轱L(fēng)力機(jī)屬于大升阻比類機(jī)械，主要是依賴升力驅(qū)動，而升力的精確模擬表明浸入邊界法可以較精確地模擬風(fēng)力機(jī)的功率。另一方面，可以利用不同網(wǎng)格尺寸下的阻力數(shù)值，通過線性外推法來預(yù)估理想條件下(即網(wǎng)格尺寸趨近于零)的阻力。具體來說，可以選定兩種不同尺寸網(wǎng)格重復(fù)模擬同一問題。以圖5中Re=10結(jié)果為例，可以使用dx=0.2，0.1這兩組最粗網(wǎng)格下的模擬結(jié)果來構(gòu)建線性的誤差曲線。而此曲線與縱軸的相交點(diǎn)，即為使用IBM方法在極限條件dx=0時(shí)的誤差。采用此線性外推方法，得到的結(jié)果誤差不到1%(圖6中數(shù)據(jù)線性外推，與坐標(biāo)原點(diǎn)非常接近)。對于Re=100的結(jié)果，采用dx=0.05，0.025這兩組數(shù)據(jù)，線性外推的結(jié)果誤差約為3%；如若采用dx=0.025，0.012 5這兩組網(wǎng)格數(shù)據(jù)線性外推，誤差則下降到1%以下?？偟膩碚f，雷諾數(shù)增大時(shí)需要更精細(xì)的網(wǎng)格來保證精度；而采用線性外推法，可以使用相對稀疏的兩種不同尺寸網(wǎng)格，得到精度遠(yuǎn)高于采用單一精細(xì)網(wǎng)格進(jìn)行模擬所得的結(jié)果。因?yàn)閷τ诙S問題，網(wǎng)格加倍時(shí)計(jì)算量上升4倍；而對于三維問題，計(jì)算量上升8倍。就圖5中的二維問題，如果采用極其細(xì)密網(wǎng)格來得到與線性外推法同樣精度的結(jié)果(誤差1%以下)，所需要的計(jì)算量則上升2到3個(gè)數(shù)量級。由此可以看出，這種簡單的線性外推法對于黏性阻力的模擬誤差糾正，是極其有效的。而這種方法能成功的理論根源，則在作者前期工作中進(jìn)行了深入闡述[14]。此處的模擬結(jié)果，則進(jìn)一步驗(yàn)證了前期理論的普適正確性。值得指出的是，這里的翼型模擬結(jié)果仍然限制于邊界層產(chǎn)生分離之前的流動情況。在更大雷諾數(shù)下，邊界層產(chǎn)生局部分離時(shí)，前述的線性外推法是否依然精確高效，將是進(jìn)一步要研究的工作。

圖5 不同網(wǎng)格尺寸dx以及攻角條件下，二維翼型阻力相對誤差

圖6 不同葉尖速比下力矩系數(shù)與功率系數(shù)

二維翼型的拱曲度(即拱起來的程度)與厚度(即最大內(nèi)切圓的直徑)是影響其升阻力的兩個(gè)主要形狀參數(shù)。此處，模擬了三個(gè)不同拱曲度的翼型，即S825，S805A，以及S809(三者拱曲度從大到小)。模擬結(jié)果表明，翼型的拱曲度越大，升力系數(shù)越大，這與文獻(xiàn)[20-21]中的研究結(jié)果一致。比較三種翼型的厚度，S809厚度最大，其次是翼型S825，最小的是翼型S805A，且翼型S809的厚度明顯大于S825和S805A，S825翼型厚度略大于S805A。翼型S805A屬于薄翼型，而S825和S809屬于厚翼型。三種不同厚度翼型在不同雷諾數(shù)及來流攻角下的模擬結(jié)果表明，當(dāng)雷諾數(shù)越小時(shí)，厚度較大的翼型升阻比越大。雷諾數(shù)變大時(shí)，在小攻角下，翼型的厚度對升阻比的影響不敏感，但隨著攻角的增大，厚度較大的翼型首先出現(xiàn)失速現(xiàn)象。

3.2 單風(fēng)力機(jī)流場模擬

基于前文介紹的方法與模型，模擬了三葉片風(fēng)力機(jī)在不同葉尖速比λ(翼尖轉(zhuǎn)動線速度與來流風(fēng)速之比)下的流動。本研究中，基于來流與葉片長度的雷諾數(shù)為2.5×104，這比工程應(yīng)用中的實(shí)際雷諾數(shù)要低兩個(gè)數(shù)量級左右。這是由于本研究中采用直接數(shù)值模擬，無法實(shí)現(xiàn)對應(yīng)于工程應(yīng)用中的高雷諾數(shù)。盡管如此，也有理由相信此處關(guān)于風(fēng)力機(jī)力矩與功率的數(shù)值模擬結(jié)果與實(shí)際工況比較吻合。這主要是因?yàn)轱L(fēng)力機(jī)升力為主型機(jī)械，而在本文模擬的雷諾數(shù)下，葉片繞流的整體結(jié)構(gòu)與工程實(shí)際非常類似，從而使得葉片上模擬的壓力分布具有較高的可信度(盡管阻力誤差較大)，最終保證葉片上的升力模擬所具有的精度。事實(shí)上，關(guān)于類似翼型的風(fēng)洞試驗(yàn)結(jié)果表明[22]，在湍流繞流條件下，不同攻角下翼型升力系數(shù)在整個(gè)實(shí)驗(yàn)范圍內(nèi)(雷諾數(shù)從6.0×104至4.6×105)幾乎完全相同，而只有阻力系數(shù)才對雷諾數(shù)呈現(xiàn)出較強(qiáng)的依賴性。

圖6給出了不同葉尖速比下葉片上的力矩系數(shù)(Cm=2M/(ρV2SR)，其中M，S，R分別是風(fēng)力機(jī)所受力矩、風(fēng)輪掃掠面積和風(fēng)輪半徑，ρ和V分別是流體的速度和來流風(fēng)速)與功率系數(shù)(CP=Cmλ，其中Cm和λ分別為力矩系數(shù)和葉尖速比)。雖然葉片阻力隨著葉尖速比增大近似線性增大(此處沒有顯示模擬結(jié)果)，但對力矩系數(shù)以及功率系數(shù)的影響則完全不同。在葉尖速比為1～3時(shí)，風(fēng)力機(jī)的功率系數(shù)加速增加，葉尖速比為3～5時(shí)，功率系數(shù)隨葉尖速比的增加變緩且在葉尖速比為5時(shí)達(dá)到最大值0.43，然后隨著葉尖速比的升高，風(fēng)力機(jī)發(fā)生失速現(xiàn)象導(dǎo)致功率系數(shù)下降，即此條件下風(fēng)力機(jī)的最佳葉尖速比為5，變化趨勢與文獻(xiàn)[18-19]的結(jié)果相同。

為了研究塔架對風(fēng)力機(jī)的影響，本文還模擬了沒有塔架的風(fēng)力機(jī)。對比兩者結(jié)果發(fā)現(xiàn)，就本文模型而言，塔架所受阻力約為總風(fēng)力機(jī)的30%，而葉片所受升力為7%。而塔架的作用，是使得葉片的上力矩減小了0.6%?？偟膩碚f，因?yàn)樗鼙旧磔^大的風(fēng)載荷，必需考慮這部分載荷的強(qiáng)度要求；而塔架對風(fēng)力機(jī)功率的影響則較小，可以忽略不計(jì)。

3.3 雙風(fēng)力機(jī)流場模擬

在風(fēng)力場中，多風(fēng)力機(jī)的氣動干涉問題，是決定風(fēng)力機(jī)布局的關(guān)鍵因素之一。一方面，希望風(fēng)力場占地盡可能小，要求風(fēng)力機(jī)布局盡可能緊密；但另一方面，為了保證單風(fēng)力機(jī)的風(fēng)能效率，又要求風(fēng)力機(jī)之間相互影響盡可能小。本文采用的浸入邊界法，可以很便捷地模擬多風(fēng)力機(jī)的氣動干涉問題。只需將前文介紹的風(fēng)力機(jī)表面離散點(diǎn)，復(fù)制放置于風(fēng)力場中不同位置，便可模擬多風(fēng)力機(jī)的流動問題。在此過程中，無需對均勻的背景歐拉網(wǎng)格進(jìn)行任何處理。

圖7給出了在葉塵速比為1，雷諾數(shù)為2.5×104條件下，兩風(fēng)力機(jī)前后間隔為3D時(shí)(D為葉片直徑)垂直截面瞬時(shí)速度云圖。從圖上可以清楚看出地面，塔架等對空氣來流的阻礙作用。為了分析前風(fēng)力機(jī)尾流的影響，圖8給出了雙風(fēng)力機(jī)在前后間隔分別為3D與5D時(shí)，沿葉片旋轉(zhuǎn)軸平均流動速度圖。第一臺風(fēng)力機(jī)后，軸線處的速度從零(近輪轂處)開始迅速增加，然后趨于穩(wěn)定；在接近第二個(gè)風(fēng)力機(jī)時(shí)快速減小到零，而在離開第二個(gè)風(fēng)力機(jī)后又快速增加。在風(fēng)力機(jī)間隔為3D條件下，第一個(gè)風(fēng)力機(jī)后尾流的最大速度恢復(fù)到空氣來流的63%。而在間隔為5D的條件下，最大速度恢復(fù)到來流的79%。因?yàn)轱L(fēng)力機(jī)的最大理想輸出功率正比于空氣來流速度的三次方，那么基于軸線上來流風(fēng)速進(jìn)行簡單估計(jì)，在風(fēng)力機(jī)間隔為3D時(shí)，第二個(gè)風(fēng)力機(jī)功率僅為第一個(gè)風(fēng)力機(jī)的25%；在間隔為5D時(shí)，第二個(gè)風(fēng)力機(jī)功率僅為第一個(gè)風(fēng)力機(jī)的一半。所以，為了保證風(fēng)力場的整體風(fēng)能效率，要盡可能避免風(fēng)力機(jī)位于其他風(fēng)力機(jī)的近尾流區(qū)。

圖7 雙風(fēng)力機(jī)速度云圖(垂直剖面)

圖8 前后雙風(fēng)力機(jī)沿軸向氣流平均速度

4 結(jié)論

浸入邊界法是一種非常常見的多相顆粒流直接數(shù)值模擬方法。本文將浸入邊界法應(yīng)用到風(fēng)力機(jī)的氣動力學(xué)數(shù)值模擬中，建立了從建模，網(wǎng)格離散，到數(shù)值模擬分析的統(tǒng)一框架模型。首創(chuàng)性地提出采用同倫變形和仿射變換的方法便捷地處理復(fù)雜葉片的建模問題。在模擬二維翼型的升阻力過程中，通過與高精度譜元方法結(jié)果對比，發(fā)現(xiàn)浸入邊界法相對空間離散具有較嚴(yán)格的一階精度。這一方面數(shù)值檢驗(yàn)并擴(kuò)展了作者的前期理論研究成果[14](即浸入邊界法在平板邊界層模擬中具有嚴(yán)格的一階精度)，證明在具有復(fù)雜壓力梯度的邊界層中，浸入邊界法仍然具有一階精度。另一方面，基于不同網(wǎng)格下的升阻力模擬結(jié)果，提出了一種簡單的理查森外推法來修正結(jié)果，從而得到遠(yuǎn)高于一階精度的模擬結(jié)果。

對不同拱曲度以及厚度的二維翼型模擬結(jié)果表明：翼型的拱曲度越大，升力系數(shù)越大；當(dāng)雷諾數(shù)越小時(shí)，厚度較大的翼型升阻比越大。雷諾數(shù)變大時(shí)，在小攻角下，翼型的厚度對升阻比的影響不敏感，但隨著攻角的增大，厚度較大的翼型首先出現(xiàn)失速現(xiàn)象。

對包含塔架的三葉片風(fēng)力機(jī)整機(jī)模擬表明，塔架本身能承受比較大比例的風(fēng)阻力載荷(此處為30%)，但塔架對風(fēng)力機(jī)葉片的轉(zhuǎn)運(yùn)力矩的氣動影響很小，可以忽略不計(jì)。而通過模擬不同葉尖速比下的風(fēng)力機(jī)運(yùn)行效率，發(fā)現(xiàn)在葉尖速比為5左右時(shí)，其有最大功率系數(shù)。

通過模擬雙風(fēng)力機(jī)在不同前后間隔下的空氣流場，發(fā)現(xiàn)在風(fēng)力機(jī)場中，前后風(fēng)力機(jī)間有很強(qiáng)的氣動干涉。間隔越短，干涉越強(qiáng)。在間隔為5倍葉片直徑條件下，空氣下游風(fēng)力機(jī)的功率只有上游風(fēng)力機(jī)功率的一半。因而，風(fēng)力機(jī)就盡可能避免處于上游風(fēng)力機(jī)的近尾流區(qū)。

本研究的最主要意義是提出一種用來模擬風(fēng)力機(jī)各種氣動作用的統(tǒng)一框架模型，可以簡單地處理葉片幾何構(gòu)型(拱曲度，厚度，扭轉(zhuǎn)角，錐度等)、部件(塔架，輪轂，機(jī)艙等)、以及多風(fēng)力機(jī)相互氣動干涉等眾多參數(shù)。在這種簡單的統(tǒng)一框架下，有可能實(shí)現(xiàn)關(guān)于風(fēng)力機(jī)氣動性能的自動優(yōu)化選型--這是正在進(jìn)行的深入研究工作。