李燕玲 胡進 周錕

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

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

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

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

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

1 研究方法

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

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

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

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

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

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

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

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

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

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

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

此處，介紹一種簡單地生成光滑漸變?nèi)~片的方法。假設(shè)在葉片最大弦長處(靠近根部)采用翼型B，而在最小弦長處(葉片尖端)采用翼型C。翼型B和C可以從通用翼型庫中選擇(國家可再生能源實驗室官方網(wǎng)站上提供了幾十種二維翼型供選擇)。每種二維翼型，都可以從其前緣到后緣沿著弧長進行參數(shù)化，即以前緣為零點到前緣的表面弧長為坐標(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ù)雜的情形，即在葉片多處位置指定不同的二維翼型，而在任意兩段翼型之間進行同倫變形來生成光滑過渡翼型。

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

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

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

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

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

2 風(fēng)力機幾何模型及離散點示例

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

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

圖2 葉片模型尺寸以及表面離散點

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

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

3 模擬結(jié)果

3.1 二維機翼繞流

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

作為對比標(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é)果對比

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

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

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

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

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

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

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

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

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

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

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

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

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

4 結(jié)論

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

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

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

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

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