收稿日期：2023-02-15

基金項(xiàng)目：國家自然科學(xué)基金（51908430；52178476；52068019）；濰坊學(xué)院博士科研啟動基金（2019BS16）；海南省自然科學(xué)基金

（522RC605；520QN231）

通信作者：李正農(nóng)（1962—），男，博士、教授，主要從事風(fēng)能資源利用及結(jié)構(gòu)抗風(fēng)方面的研究。zhn88@263.net

DOI：10.19912/j.0254-0096.tynxb.2023-0154 文章編號：0254-0096（2023）12-0330-11

摘 要：采用無人機(jī)測風(fēng)手段實(shí)測獲得某沿海風(fēng)電場的風(fēng)場特征，將所測量的風(fēng)場參數(shù)應(yīng)用到后續(xù)數(shù)值模擬中。以某2 MW風(fēng)力機(jī)葉片全尺寸模型為研究對象，基于Workbench平臺建立流固耦合系統(tǒng)進(jìn)行計(jì)算，通過對比近尾流風(fēng)速分布的數(shù)值模擬結(jié)果與實(shí)測結(jié)果來驗(yàn)證數(shù)值計(jì)算的準(zhǔn)確性，進(jìn)而計(jì)算額定風(fēng)速下旋轉(zhuǎn)葉片的揮舞變形和等效應(yīng)力等響應(yīng)。結(jié)果表明：風(fēng)力機(jī)葉片的揮舞變形在0°方位角時最大，最大變形為1916.4 mm，接近靜力加載試驗(yàn)最大揮舞工況結(jié)果2299 mm；相同方位角下的揮舞變形沿著葉片展向呈非線性增大，0.60倍風(fēng)輪半徑處變形增長速度明顯提高。葉片等效應(yīng)力集中區(qū)位于大梁和前緣側(cè)交界處，最大值出現(xiàn)在60°方位角下0.78倍風(fēng)輪半徑處，同時在風(fēng)輪高度方向上呈非對稱分布。

關(guān)鍵詞：風(fēng)力機(jī)葉片；流固耦合；風(fēng)速分布；揮舞變形；等效應(yīng)力

中圖分類號：TK83 """" 文獻(xiàn)標(biāo)志碼：A

0 引 言

在“雙碳”時代背景下風(fēng)能領(lǐng)域?qū)Πl(fā)電效率的要求越來越高，葉片作為風(fēng)力機(jī)捕獲風(fēng)能的關(guān)鍵部件愈發(fā)趨向細(xì)長柔化演變。風(fēng)力機(jī)葉片是一個復(fù)雜的流-固耦合系統(tǒng)［1］，細(xì)長柔結(jié)構(gòu)葉片運(yùn)行過程中在重力、離心力和氣動力等交變載荷作用下更易發(fā)生非線性的變形和結(jié)構(gòu)振動［2］。同時，沿海風(fēng)力機(jī)多位于臺風(fēng)登陸區(qū)，自1998年至今每年7—9月份平均有5.4次臺風(fēng)從中國東南沿海地區(qū)登陸［3］。在極端風(fēng)環(huán)境的影響下，較高的風(fēng)速和湍流強(qiáng)度會引起風(fēng)力機(jī)葉片的無規(guī)則振動，失效現(xiàn)象頻發(fā)［4］，嚴(yán)重情況下可能發(fā)生葉片的撕裂和折斷［5］，而這些問題的解決都必須建立在對葉片動態(tài)響應(yīng)準(zhǔn)確計(jì)算的基礎(chǔ)上。因此，研究沿海風(fēng)力機(jī)葉片的流固耦合特性，準(zhǔn)確計(jì)算葉片的動態(tài)響應(yīng)，對風(fēng)力機(jī)葉片的設(shè)計(jì)、運(yùn)維都具有重要的工程意義。

目前，國內(nèi)外學(xué)者在求解風(fēng)力機(jī)葉片氣動載荷方面主要采用葉素動量理論、渦尾跡方法、致動模型和計(jì)算流體動力學(xué)（computational fluid dynamics， CFD）方法，在建立結(jié)構(gòu)模型方面主要運(yùn)用模態(tài)法、多體動力學(xué)方法和3D有限元模型，進(jìn)而進(jìn)行流固耦合計(jì)算葉片的動態(tài)響應(yīng)。陳佳慧等［6］結(jié)合葉素動量理論（blade element-momentum theory，BEM）和盒形梁理論建立風(fēng)力機(jī)葉片的動態(tài)響應(yīng)分析方法。趙鷹等［7］提出一種耦合BEM理論和復(fù)合材料鋪層結(jié)構(gòu)建模的風(fēng)力機(jī)葉片氣彈模型，計(jì)算氣彈耦合對風(fēng)力機(jī)功率、葉片結(jié)構(gòu)響應(yīng)等的影響。呂品等［8］采用分區(qū)弱耦合方法建立基于幾何精確梁理論和自由渦尾跡理論的風(fēng)力機(jī)葉片氣彈模型，計(jì)算了NREL 5 MW風(fēng)力機(jī)葉片的固有頻率及結(jié)構(gòu)響應(yīng)，并與FAST軟件的計(jì)算結(jié)果進(jìn)行對比驗(yàn)證。李德源等［2］將翼型振動速度和扭轉(zhuǎn)變形引入到風(fēng)力機(jī)翼型動態(tài)失速模型中，結(jié)合修正的BEM理論，建立葉片的多體氣彈耦合方程，得出氣彈耦合模型的葉根揮舞力矩與葉片位移均值均小于未考慮耦合的計(jì)算結(jié)果。楊斌等［9］提出一種基于自由渦尾跡方法和FAST結(jié)構(gòu)模塊耦合的風(fēng)力機(jī)葉片氣彈模型，并模擬了風(fēng)切變下的葉片動態(tài)響應(yīng)。Lee等［10］將Xfoil中的翼型壓力系數(shù)積分后得到氣動載荷，施加在葉片的3D有限元模型上，迭代計(jì)算了葉片結(jié)構(gòu)響應(yīng)。以上方法在氣動力計(jì)算方面多采用經(jīng)驗(yàn)?zāi)Ｐ?，?jì)算效率高但難以滿足復(fù)雜風(fēng)場的計(jì)算需要。

郭同慶等［11］建立CFD/CSD時域模型，耦合了N-S方程和結(jié)構(gòu)方程，結(jié)合一維梁模型進(jìn)行動態(tài)響應(yīng)計(jì)算，結(jié)果與渦尾跡法較為接近。戴麗萍等［12］采用非定常雷諾平均N-S方程和有限元松耦合的方法，開展Tj?reborg風(fēng)力機(jī)葉片模型的動態(tài)響應(yīng)計(jì)算，得到偏航角對風(fēng)力機(jī)功率、葉片變形和應(yīng)變的影響規(guī)律。Sayed等［13］提出一種基于流固耦合（fluid-structure interaction，F(xiàn)SI）的氣動彈性仿真方法，采用CFD方法求解葉片氣動載荷，利用CSD求解器Carat++計(jì)算了基于非線性Timoshenko梁單元的彈性變形，指出柔性葉片的的振動幅值顯著減小，在距離葉根90%的位置處最大可降低22%。以上研究中，葉片的結(jié)構(gòu)模型均采用簡化模型，無法準(zhǔn)確反映葉片的扭角、剛度中心等結(jié)構(gòu)特征及鋪層信息，從而降低了葉片動態(tài)響應(yīng)的計(jì)算精度。

李媛等［14］采用MpCCI平臺耦合CFD和CSD軟件進(jìn)行風(fēng)力機(jī)葉片流固耦合分析，得到3個風(fēng)速下葉片的變形及應(yīng)力分布，并指出高風(fēng)速下葉片流固耦合效應(yīng)不容忽視。劉海鋒等［15］指出考慮雙向流固耦合后，NREL 5 MW風(fēng)力機(jī)的近尾流區(qū)受葉片變形的影響速度虧損減小，遠(yuǎn)尾流區(qū)的速度恢復(fù)更快。趙元星等［16］考慮不同風(fēng)切變的影響，針對兩種翼型分析流固耦合作用下的應(yīng)力分布規(guī)律，并指出葉片在方位角等于30°時出現(xiàn)最大應(yīng)力值。宋力等［17］以縮尺的1.5 MW風(fēng)力機(jī)葉片為研究對象，基于Workbench平臺進(jìn)行流固耦合模擬計(jì)算，得出葉片位移變形與載荷、轉(zhuǎn)速的影響規(guī)律。以上CFD數(shù)值模擬中，較少考慮風(fēng)力機(jī)在實(shí)際風(fēng)場中的來流和尾流信息，且未經(jīng)過實(shí)測數(shù)據(jù)的驗(yàn)證。除了理論研究外，還可采用風(fēng)洞試驗(yàn)的方法對風(fēng)力機(jī)葉片動態(tài)響應(yīng)進(jìn)行研究?？率捞玫龋?8］以NREL 15 MW風(fēng)力機(jī)葉片為研究對象，開展縮尺模型的測振、測力風(fēng)洞試驗(yàn)，得到了葉片顫振臨界風(fēng)速后的響應(yīng)特性和能量耗散機(jī)制。周興［19］搭建一種風(fēng)力機(jī)葉片動力學(xué)分析實(shí)驗(yàn)臺，得到不同風(fēng)速下葉片的應(yīng)力應(yīng)變分布規(guī)律。但由于風(fēng)力機(jī)葉片沿展向的尺寸變化較大，尤其是葉尖處進(jìn)行幾何縮尺后與實(shí)際葉片結(jié)構(gòu)存在較大差別，如何在風(fēng)力機(jī)葉片的氣彈模型風(fēng)洞試驗(yàn)中選擇合適的相似比，仍是目前制約風(fēng)洞試驗(yàn)開展的關(guān)鍵問題。

綜上，數(shù)值模擬方法是開展風(fēng)力機(jī)葉片流固耦合研究的有效手段，但數(shù)值模擬的合理性仍缺少實(shí)測風(fēng)場特性的驗(yàn)證。為此，本文首先采用無人機(jī)測風(fēng)手段獲得某沿海風(fēng)電場的風(fēng)場數(shù)據(jù)，將其作為后續(xù)數(shù)值模擬的初始條件。以某2 MW風(fēng)力機(jī)葉片全尺寸模型為研究對象，采用Workbench中的Fluent模塊和瞬態(tài)結(jié)構(gòu)模塊建立流固耦合系統(tǒng)進(jìn)行數(shù)值模擬，并用實(shí)測擬合的風(fēng)力機(jī)下游近尾流風(fēng)速分布值驗(yàn)證數(shù)值計(jì)算的準(zhǔn)確性，進(jìn)而計(jì)算額定風(fēng)速下旋轉(zhuǎn)葉片的變形、應(yīng)力等動態(tài)響應(yīng)，以期為在役風(fēng)力機(jī)葉片的運(yùn)維提供參考。

1 風(fēng)場實(shí)測

1.1 實(shí)測概況

為獲得沿海風(fēng)電場的實(shí)際風(fēng)場特征，選擇浙江某沿海風(fēng)電場為實(shí)測場地，場內(nèi)布置20臺單機(jī)容量為2 MW的風(fēng)力機(jī)，位置及周圍環(huán)境如圖1所示。

實(shí)測設(shè)備為2臺搭載2D超聲風(fēng)速儀的無人機(jī)，前期經(jīng)過風(fēng)洞試驗(yàn)和現(xiàn)場實(shí)測驗(yàn)證了該實(shí)測手段的準(zhǔn)確性。通過設(shè)計(jì)合理的實(shí)測方案［20-21］，將2臺無人機(jī)同時置于風(fēng)力機(jī)前側(cè)，無人機(jī)A所在的輪轂高度處作為參考點(diǎn)，無人機(jī)B逐次測量10～120 m共12個高度處的數(shù)據(jù)，每個高度處采集時間為2 min，如圖2所示。

本研究在2020年9月和2021年6月兩次現(xiàn)場實(shí)測過程中共采集11條來流風(fēng)速剖面，采用平均風(fēng)速作為研究對象，以輪轂中心高度為參考高度計(jì)算擬合得到11條風(fēng)速剖面的剖面指數(shù)α，并與《建筑結(jié)構(gòu)荷載規(guī)范》（GB 50009—2012）［22］中規(guī)定的A類地貌指數(shù)值0.12進(jìn)行對比，如圖3所示，是實(shí)測數(shù)據(jù)擬合得到的風(fēng)速剖面指數(shù)與0.12之間的誤差棒圖，最終確定實(shí)測風(fēng)力機(jī)所處地貌接近于A類地貌。

需說明的是，在實(shí)測過程中通過實(shí)時數(shù)據(jù)傳輸來觀察來流風(fēng)向的變化，使2臺無人機(jī)始終處于相同來流風(fēng)向范圍內(nèi)。以實(shí)測獲得的1條風(fēng)速剖面為例進(jìn)行說明，圖4為該條來流風(fēng)速剖面2臺無人機(jī)采集的風(fēng)速、風(fēng)向分布統(tǒng)計(jì)圖，來流平均風(fēng)速集中在6.00～14.55 m/s，風(fēng)向?yàn)楸憋L(fēng)風(fēng)向，主要集中在342.34°～2.64°的風(fēng)向范圍，水平向平均風(fēng)向標(biāo)準(zhǔn)差平均為5.19°［20］，可認(rèn)為實(shí)測過程中的風(fēng)向較為穩(wěn)定。

1.2 實(shí)測結(jié)果

通過現(xiàn)場實(shí)測樣本數(shù)據(jù)分析獲得了風(fēng)電場的平均風(fēng)速剖面和平均湍流強(qiáng)度剖面等參數(shù)，如圖5所示。其中，來流平均風(fēng)速剖面公式與IEC 61400-1-2019中一致［20］，參考高度為輪轂中心高度。本次實(shí)測獲得的平均風(fēng)速擬合指數(shù)為0.12，接近于《建筑結(jié)構(gòu)荷載規(guī)范》（GB 50009—2012）［22］中規(guī)定的A類地貌指數(shù)值，符合實(shí)測風(fēng)電場所處的沿海灘涂地貌。在0～120 m高度Z范圍內(nèi)湍流強(qiáng)度分布在0.10～0.17范圍內(nèi)。采用編譯UDF自定義函數(shù)將上述實(shí)測風(fēng)場參數(shù)應(yīng)用到數(shù)值模擬中。

2 數(shù)值模擬計(jì)算

2.1 風(fēng)力機(jī)葉片三維幾何模型

本文選取的2 MW風(fēng)力機(jī)葉片截面翼型主要由DU系列和NACA-64618構(gòu)成，其中葉片展向長度9 m位置處的翼型截面弦長[c]最大，為3.5 m，扭角[β]為11.61°，如圖6所示為葉片的弦長和扭角。葉片其他基本參數(shù)詳見表1。

采用SpaceClaim建立風(fēng)力機(jī)風(fēng)輪的三維全尺寸模型，如圖7所示。將每個葉片沿展向分為12段，風(fēng)輪旋轉(zhuǎn)中心設(shè)為坐標(biāo)原點(diǎn)，[X]軸平行于風(fēng)輪平面，正方向表示葉片翼型截面由前緣指向后緣；[Y]軸平行于來流方向，正方向指向風(fēng)輪下游；[Z]軸為高度方向，向上為正方向，[Z=-80] m處為地面。

2.2 Fluent流體模型的建立

2.2.1 計(jì)算域的劃分

根據(jù)全尺寸葉片模型建立流體域，包括內(nèi)部旋轉(zhuǎn)域和外部靜止域，計(jì)算域設(shè)置為圓柱體類型。計(jì)算域的尺寸劃分參考文獻(xiàn)［23］，并考慮風(fēng)切變效應(yīng)、尾流效應(yīng)等因素對計(jì)算域進(jìn)行一定程度放大。如圖8所示，風(fēng)輪直徑[D=96]m，旋轉(zhuǎn)域直徑[1.2D]，厚度為[0.6D]，旋轉(zhuǎn)域中心為輪轂中心。風(fēng)輪中心距離靜止域入口[3.0D]，距離出口[10.0D]，靜止域高度為[4.2D]。

2.2.2 網(wǎng)格劃分

采用ICEM對外部靜止域和內(nèi)部旋轉(zhuǎn)域進(jìn)行結(jié)構(gòu)化網(wǎng)格劃分，通過O-block劃分靜止域和旋轉(zhuǎn)域的進(jìn)口、出口和壁面，并考慮邊界層的影響，局部加密6層網(wǎng)格，同時也對區(qū)域交界面進(jìn)行網(wǎng)格加密，如圖9a和圖9b所示。旋轉(zhuǎn)域內(nèi)葉片網(wǎng)格采用120°周期性網(wǎng)格劃分，分塊處理，如圖9c所示。整體單元數(shù)量為14663658（1466萬），整體節(jié)點(diǎn)數(shù)量為14964079（1496萬）。

此外，開展了網(wǎng)格無關(guān)性驗(yàn)證，將1496萬計(jì)算節(jié)點(diǎn)與1692萬和1826萬計(jì)算節(jié)點(diǎn)的網(wǎng)格進(jìn)行對比，如表2所示。由表2中的相對誤差（（設(shè)計(jì)功率-計(jì)算功率）/設(shè)計(jì)功率×100%）可見，當(dāng)節(jié)點(diǎn)總數(shù)為1496萬時，計(jì)算功率與設(shè)計(jì)功率的相對誤差為3.80%，與節(jié)點(diǎn)數(shù)更多的1692萬和1826萬網(wǎng)格相差較小。因此，綜合考慮網(wǎng)格的計(jì)算效率和計(jì)算精度，本文中采用的1496萬節(jié)點(diǎn)數(shù)網(wǎng)格是可行的。

2.2.3 Fluent中的相關(guān)設(shè)置

采用偽瞬態(tài)計(jì)算方法，選擇標(biāo)準(zhǔn)k-ε模型和標(biāo)準(zhǔn)壁面函數(shù)?？紤]葉片的旋轉(zhuǎn)效應(yīng)，利用運(yùn)動坐標(biāo)系設(shè)置風(fēng)輪繞[Y]軸順時針方向的旋轉(zhuǎn)速度[ω=1.78]rad/s。邊界條件采用速度入口和壓力出口，無滑移壁面，在旋轉(zhuǎn)域和外流場域之間設(shè)置網(wǎng)格交界面。速度入口基于實(shí)測獲得的來流平均風(fēng)速剖面采用UDF編譯加載，輪轂高度處風(fēng)速采用額定風(fēng)速10.5 m/s，風(fēng)切變指數(shù)為0.12，如圖10所示。

流體計(jì)算方法采用壓力-速度耦合的Simple算法，通量選用Rhie-Chow位移插值，壓力采用二階算法，湍流動能和時間項(xiàng)離散格式均采用二階迎風(fēng)格式。

2.3 結(jié)構(gòu)模型建立

風(fēng)力機(jī)葉片結(jié)構(gòu)建模除考慮葉片具有彎扭耦合特性外，還要將葉片的材料各向異性、鋪層設(shè)計(jì)等考慮在內(nèi)，才能獲得更精準(zhǔn)的結(jié)構(gòu)動力響應(yīng)特性［24］。風(fēng)力機(jī)葉片結(jié)構(gòu)主要包括大梁、前緣、后緣及腹板，其中主要的受力構(gòu)件是大梁，吸力面和壓力面的大梁由抗剪腹板連接，如圖11所示。由于風(fēng)力機(jī)葉片是復(fù)合材料，是由多種材料按照不同的方向鋪貼而成，為了簡化建模，鋪層厚度按照大梁、前后緣和腹板分區(qū)布置，葉片厚度從葉根51 mm漸變到葉尖5 mm，大梁從37 mm漸變到15 mm，抗剪腹板從53 mm漸變到35 mm，建模中葉片采用Shell99殼單元進(jìn)行模擬。復(fù)合材料的特性參數(shù)列于表3中。

結(jié)構(gòu)模型采用非結(jié)構(gòu)化網(wǎng)格，單元最小尺寸為0.3 m，共計(jì)25697個單元和11994個節(jié)點(diǎn)，如圖12所示。

2.4 流固耦合計(jì)算及驗(yàn)證

2.4.1 流固耦合計(jì)算設(shè)置

采用單向流固耦合的方法計(jì)算，將流體求解得到的表面壓力導(dǎo)入結(jié)構(gòu)模型中進(jìn)行求解，流體和結(jié)構(gòu)計(jì)算均采用瞬態(tài)計(jì)算，總計(jì)算時長為30 s，時間步長為0.01 s，對應(yīng)風(fēng)輪旋轉(zhuǎn)約1°，收斂殘差值達(dá)到1.0×10-6。結(jié)構(gòu)模型的計(jì)算同時考慮重力和離心力的作用，轉(zhuǎn)速為1.78 rad/s。

2.4.2 計(jì)算有效性驗(yàn)證

為驗(yàn)證數(shù)值模擬結(jié)果的準(zhǔn)確性，提取風(fēng)力機(jī)葉片的渦核旋轉(zhuǎn)區(qū)域示意圖和風(fēng)輪中心平面的風(fēng)速分布云圖進(jìn)行分析，如圖13所示。由圖13a可得，在風(fēng)輪后方近尾流區(qū)可明顯地

觀察到從葉尖和輪轂處有大量的渦脫落，且葉尖渦呈螺旋狀，這與文獻(xiàn)［25］提出的近尾流區(qū)葉尖渦的特征是一致的。如圖13b所示，風(fēng)切變及風(fēng)輪的旋轉(zhuǎn)效應(yīng)使得風(fēng)輪旋轉(zhuǎn)平面內(nèi)（0D）上、下半風(fēng)輪中心平面速度呈非對稱分布，在旋轉(zhuǎn)域區(qū)間存在顯著的速度虧損。

為進(jìn)一步驗(yàn)證數(shù)值模擬計(jì)算的準(zhǔn)確性，選取額定風(fēng)速10.5 m/s下流固耦合計(jì)算的尾流風(fēng)速剖面與實(shí)測尾流風(fēng)速剖面進(jìn)行對比。前期對無人機(jī)實(shí)測獲得的大量數(shù)據(jù)進(jìn)行擬合，提出預(yù)測風(fēng)力機(jī)近尾流風(fēng)速剖面的公式［20］。本文利用此公式得到額定風(fēng)速10.5 m/s風(fēng)力機(jī)下游2D處的風(fēng)速分布，并與數(shù)值模擬結(jié)果進(jìn)行對比，如圖14所示。

在圖14a中給出了垂直方向來流風(fēng)速數(shù)值模擬和實(shí)測結(jié)果的對比，在風(fēng)輪范圍內(nèi)風(fēng)速的相對差值（（實(shí)測結(jié)果-數(shù)值模擬）/實(shí)測結(jié)果×100%）在2.15%～14.75%以內(nèi)；圖14b中對比了水平方向數(shù)值模擬與實(shí)測結(jié)果的對比，最大相對誤差在輪轂中軸線為12.73%，其余均在10%之內(nèi)，且尾流分布規(guī)律較為接近。目前研究中大多利用轉(zhuǎn)矩、推力等驗(yàn)證CFD數(shù)值模擬的準(zhǔn)確性，朱瑩［26］對比CFD模型和FAST計(jì)算的轉(zhuǎn)矩相差5%之內(nèi)，推力相差10%之內(nèi)，由此認(rèn)為CFD模型可信。由于實(shí)測的復(fù)雜性，利用實(shí)測結(jié)果來驗(yàn)證數(shù)值模擬準(zhǔn)確性的研究較少，劉鑫等［27］的CFD數(shù)值模擬的結(jié)果雖然與激光雷達(dá)實(shí)測獲得的數(shù)據(jù)有差距但變化規(guī)律基本一致，可較好地反映尾流速度變化規(guī)律。綜上，本文的數(shù)值模擬與實(shí)測數(shù)據(jù)規(guī)律方向較為接近，且結(jié)果值相差較小，可用于后續(xù)的葉片動態(tài)響應(yīng)計(jì)算。

3 結(jié)果討論

3.1 葉片揮舞變形分析

風(fēng)力機(jī)在運(yùn)行中受到復(fù)雜載荷的作用導(dǎo)致葉片變形，其中葉片揮舞方向的變形較大。因此，為確保葉片與塔筒之間的凈空距離，需對流固耦合條件下風(fēng)力機(jī)葉片揮舞方向的變形進(jìn)行重點(diǎn)研究。定義當(dāng)葉片旋轉(zhuǎn)至豎直向上時方位角為0°，葉片順時針旋轉(zhuǎn)，旋轉(zhuǎn)方向與方位角的增加方向相同。為研究額定風(fēng)速10.5 m/s時旋轉(zhuǎn)葉片不同位置[r]處的揮舞方向變形，本文沿大梁與前緣的連接線在距離葉根6個展向位置處布置了監(jiān)測點(diǎn)，采用無量綱的表達(dá)形式，分別為[0.07R、0.20R、0.40R、0.60R、0.80R]和[0.97R]（[R]為風(fēng)輪半徑）。

本文提取葉片的揮舞變形分布云圖如圖15所示，可得風(fēng)力機(jī)葉片的最大揮舞變形出現(xiàn)在0°方位角的葉尖處，最大值1916.4 mm，占葉片展向長度的4.12%。

對于變速變槳風(fēng)力機(jī)，未達(dá)到額定風(fēng)速之前葉片的氣動載荷隨來流風(fēng)速的增加而增大，而當(dāng)來流風(fēng)速大于額定風(fēng)速時會進(jìn)行變槳操作，以限制氣動載荷的大小。因此，額定風(fēng)速工況下的氣動載荷最大，變形也最大。據(jù)此，本文求解的額定風(fēng)速下的葉片揮舞變形最大值1916.4 mm為該葉片最大揮舞變形。為驗(yàn)證這一結(jié)輪，將數(shù)值計(jì)算得到的0°方位角下葉片的揮舞變形與廠家的靜力加載試驗(yàn)結(jié)果（最大值為2299 mm）進(jìn)行了對比，如圖16所示?？傻茫瑪?shù)值計(jì)算結(jié)果與靜力加載試驗(yàn)中最大揮舞工況結(jié)果沿葉片展向的變化規(guī)律一致，兩者最大值相差約16%（（試驗(yàn)值-數(shù)值計(jì)算結(jié)果）/試驗(yàn)值×100%），滿足工程需求，同時也驗(yàn)證了結(jié)構(gòu)模塊計(jì)算結(jié)果的準(zhǔn)確性。進(jìn)一步分析圖16可得出，葉片的揮舞變形沿展向0.60R處開始呈現(xiàn)快速增大的趨勢，在葉尖處最大，說明葉片受到氣動力作用時，葉尖受到的擾動最大。以上結(jié)論可為葉片的預(yù)彎設(shè)計(jì)提供參考。

3.2 葉片等效應(yīng)力分析

強(qiáng)風(fēng)影響下，葉片在較大風(fēng)荷載、重力和離心力共同作用下會發(fā)生較大的變形，使得葉片表面產(chǎn)生應(yīng)力集中現(xiàn)象，嚴(yán)重時可導(dǎo)致葉片變形甚至?xí)a(chǎn)生結(jié)構(gòu)斷裂［28］。因此，研究在額定風(fēng)速下葉片的應(yīng)力分布特性對分析葉片損傷及運(yùn)行具有重要工程意義。

本文提取了葉片的等效應(yīng)力（Von-Mises）分布云圖，如圖17所示，可得葉片等效應(yīng)力（后文簡稱“應(yīng)力”）的最大值119.28 MPa出現(xiàn)在葉片展向約[0.78R]位置處，靠近葉尖，應(yīng)力最小值出現(xiàn)在葉根處。通過圖17b的局部視圖可發(fā)現(xiàn)，應(yīng)力最大值出現(xiàn)在葉片迎風(fēng)面大梁和前緣側(cè)的交界處。前期開展了葉片的有限元分析［29］，按照靜力加載試驗(yàn)進(jìn)行應(yīng)力計(jì)算，有限元分析結(jié)果如圖17c所示，應(yīng)力最大值也出現(xiàn)在前緣側(cè)與大梁連接的位置處，這進(jìn)一步驗(yàn)證了本文采用的流固耦合計(jì)算方法是合理可行的。

葉片在旋轉(zhuǎn)一周的過程中，應(yīng)力分布會隨著方位角發(fā)生變化。圖18給出了不同方位角下的葉片應(yīng)力分布，為了反映應(yīng)力最大值的所在位置，圖18中增加了展向位置[0.78R]。如圖18所示，給出了7個展向位置處應(yīng)力0°～360°的應(yīng)力分布，并用星形符號標(biāo)注了每一個展向位置處的最大值?？梢?，相同位置處，當(dāng)葉片旋轉(zhuǎn)至60°方位角時，葉片應(yīng)力最大，且風(fēng)輪上半部分和下半部分的應(yīng)力分布不對稱，上半部分大于下半部分。為了分析風(fēng)剪切的影響，圖18中將風(fēng)切變和均勻風(fēng)速入流兩種來流條件下的應(yīng)力分布進(jìn)行對比。結(jié)果表明，當(dāng)風(fēng)力機(jī)輪轂高度處的來流風(fēng)速一致時，風(fēng)切變?nèi)肓鲿r葉片應(yīng)力比均勻入流情況增大約1.32%（（風(fēng)切變?nèi)肓鳁l件119.28 MPa-均勻入流條件117.71 MPa）/風(fēng)切變?nèi)肓鳁l件），這與文獻(xiàn)［30］中指出的風(fēng)切變效應(yīng)使得應(yīng)力增大1.92%相接近。

為進(jìn)一步分析沿葉片展向不同位置處應(yīng)力應(yīng)變情況，本文以60°、90°、180°、270°和330°方位角為例提取7個位置處的應(yīng)力情況，兩種入流條件下的應(yīng)力展向分布如圖19所示?？傻?，切變和均勻風(fēng)速入流條件下的葉片應(yīng)力沿葉根向葉尖展向方向呈先增大后減小趨勢，葉片應(yīng)力集中區(qū)域主要出現(xiàn)在[0.60R～0.80R]范圍內(nèi)，在[0.78R]處出現(xiàn)最大值，這與文獻(xiàn)［31］中應(yīng)力集中區(qū)域位于葉片展向長度50%～85%范

圍內(nèi)的結(jié)論相近，這是由于葉片的揮舞變形在大于[0.60R]處（圖16）開始快速增大導(dǎo)致葉片應(yīng)力升高。葉片應(yīng)力集中位置容易形成疲勞損傷，因此在葉片設(shè)計(jì)時應(yīng)重點(diǎn)關(guān)注此區(qū)域。

同時，由圖19可得，節(jié)點(diǎn)應(yīng)力在葉片展向[0.40R]處有局部增大的趨勢，這是由于風(fēng)力機(jī)葉片在葉片展向22 m（[0.46R]）處開始產(chǎn)生預(yù)彎，該位置區(qū)域存在一定的應(yīng)力增加。如圖20所示為葉片預(yù)彎設(shè)計(jì)示意圖，預(yù)彎使得葉片展向沿著[Y]方向（圖7）產(chǎn)生[yi]的位移，并與風(fēng)輪旋轉(zhuǎn)平面中產(chǎn)生一個夾角[θ]。

4 結(jié) 論

本文以某2 MW風(fēng)力機(jī)葉片為研究對象，基于無人機(jī)現(xiàn)場實(shí)測，結(jié)合Workbench平臺中的Fluent和Mechanical模塊開展沿海地區(qū)風(fēng)力機(jī)葉片在流固耦合作用下的結(jié)構(gòu)響應(yīng)研究，得到葉片揮舞變形和應(yīng)力的變化規(guī)律，主要結(jié)論如下：

1）葉片揮舞變形最大值出現(xiàn)在0°方位角下的葉尖處，為1916.4 mm，約為葉片展向長度的4.12%，與靜力加載試驗(yàn)中最大揮舞工況的葉尖變形2299 mm相差約16%，滿足工程計(jì)算需要；相同方位角下?lián)]舞變形沿葉片展向長度呈非線性變化，從0.60倍風(fēng)輪半徑處變化曲線斜率明顯增大。

2）葉片應(yīng)力最大值119.28 MPa出現(xiàn)在葉片迎風(fēng)面0.78倍風(fēng)輪半徑處，主要位于前緣側(cè)與大梁的交界處，與靜力加載試驗(yàn)得到的應(yīng)力最大值出現(xiàn)位置一致。且受風(fēng)切變、風(fēng)輪旋轉(zhuǎn)及重力的共同影響，旋轉(zhuǎn)葉片應(yīng)力分布在風(fēng)輪上下部分呈非對稱特性，風(fēng)輪上部應(yīng)力大于風(fēng)輪下部。

3）葉片應(yīng)力受風(fēng)剪切的影響，比均勻風(fēng)速入流條件下約增大1.32%；受葉片預(yù)彎的影響，葉片應(yīng)力在0.40R處存在局部增大。

4）本文采用的風(fēng)力機(jī)葉片動態(tài)響應(yīng)的流固耦合計(jì)算方法經(jīng)過了實(shí)測尾流風(fēng)速分布數(shù)據(jù)的驗(yàn)證，揮舞變形和應(yīng)力計(jì)算結(jié)果也與靜力加載試驗(yàn)較為一致，該模型可用于后續(xù)的葉片損傷分析研究。

［參考文獻(xiàn)］

［1］"""" TANG D， BAO S Y， LUO L J， et al. Study on the aeroelastic responses of a wind turbine using a coupled multibody-FVW method［J］. Energy， 2017， 141： 2300-2313.

［2］""nbsp;" 李德源， 汪顯能， 莫文威， 等. 非定常條件下風(fēng)力機(jī)柔性葉片氣彈耦合分析［J］. 太陽能學(xué)報(bào)， 2017， 38（4）： 966-975.

LI D Y， WANG X N， MO W W， et al. Aeroelastic coupling analysis of flexible blades of wind turbine under unsteady"" conditions［J］." Acta"" energiae"" solaris"" sinica， 2017， 38（4）： 966-975.

［3］"""" 溫州臺風(fēng)網(wǎng)，http：//www.wztf121.com/analyse/wzs.

Wenzhou Typhoon Network，http：//www.wztf121.com/analyse/wzs.

［4］"""" 樓文娟， 余江， 潘小濤. 風(fēng)力機(jī)葉片揮舞擺振氣彈失穩(wěn)分析［J］. 工程力學(xué)， 2015， 32（11）： 236-242.

LOU W J， YU J， PAN X T. Calculating for aerodynamic stability response of wind turbine blade in flapwise and edgewise" direction［J］. Engineering mechanics， 2015， 32（11）： 236-242.

［5］"""" 王景全， 陳政清. 試析海上風(fēng)機(jī)在強(qiáng)臺風(fēng)下葉片受損風(fēng)險與對策： 考察紅海灣風(fēng)電場的啟示［J］. 中國工程科學(xué)， 2010， 12（11）： 32-34.

WANG J Q， CHEN Z Q. Analysis of risks and measures on the blade damage of offshore wind turbine during strong typhoons：" enlightenment"" from" Red" Bay" wind" farm［J］. Engineering sciences， 2010， 12（11）： 32-34.

［6］"""" 陳佳慧， 王同光. 考慮氣動彈性的風(fēng)力機(jī)葉片性能分析［J］. 空氣動力學(xué)學(xué)報(bào)， 2011， 29（3）： 396-400.

CHEN J H， WANG T G. Wind turbine performance analysis"" with"" aeroelastic" effect［J］." Acta" aerodynamica sinica， 2011， 29（3）： 396-400.

［7］"""" 趙鷹， 廖猜猜， 秦志文， 等. 基于高精度有限元模型的葉片氣彈耦合分析［J］. 工程熱物理學(xué)報(bào)， 2017， 38（4）： 748-753.

ZHAO Y， LIAO C C， QIN Z W， et al. Aeroelastic coupling analysis of blades based on high precision finite element model［J］. Journal of engineering thermophysics， 2017， 38（4）： 748-753.

［8］"""" 呂品， 廖明夫， 尹堯杰. 大型風(fēng)力機(jī)葉片非線性流固耦合模型［J］. 太陽能學(xué)報(bào)， 2017， 38（8）： 2126-2135.

LYU P， LIAO M F， YIN Y J. A structural-aerodynamic coupled method for nonlinear aeroelastic response of large-scaled hawt［J］. Acta energiae solaris sinica， 2017， 38（8）： 2126-2135.

［9］"""" 楊斌， 戴麗萍， 王曉東， 等. 大型風(fēng)力機(jī)葉片非線性氣彈模型［J］. 工程熱物理學(xué)報(bào)， 2021， 42（3）： 626-632.

YANG B， DAI L P， WANG X D， et al. Nonlinear aeroelastic model of large wind turbine blades［J］. Journal of engineering thermophysics， 2021， 42（3）： 626-632.

［10］""" LEE Y J， JHAN Y T， CHUNG C H. Fluid-structure interaction of FRP wind turbine blades under aerodynamic effect［J］. Composites part B： engineering， 2012， 43（5）： 2180-2191.

［11］""" GUO T Q， LU Z L， TANG D， et al. A CFD/CSD model for aeroelastic calculations of large-scale wind turbines［J］. Science China technological sciences， 2013， 56（1）： 205-211.

［12］""" DAI L P， ZHOU Q， ZHANG Y W， et al. Analysis of wind turbine blades aeroelastic performance under yaw conditions［J］. Journal of wind engineering and industrial aerodynamics， 2017， 171： 273-287.

［13］""" SAYED M， LUTZ T， KR?MER E， et al. Aeroelastic analysis of 10 MW wind turbine using CFD-CSD explicit FSI-coupling approach［J］. Journal of fluids and structures， 2019， 87： 354-377.

［14］""" 李媛， 康順， 王建錄， 等. 2.5 MW風(fēng)力機(jī)葉片流固耦合數(shù)值模擬［J］. 工程熱物理學(xué)報(bào)， 2013， 34（1）： 71-74.

LI Y， KANG S， WANG J L， et al. Numerical simulation of fluid-structure coupling of 2.5 MW wind turbine blades［J］. Journal of engineering thermophysics， 2013， 34（1）： 71-74.

［15］""" 劉海鋒， 孫凱， 胡丹梅. 大型風(fēng)力機(jī)尾跡雙向流固耦合特性分析［J］. 可再生能源， 2015， 33（11）： 1664-1673.

LIU H F， SUN K， HU D M. Characteristic analysis of large"" wind"" turbine"" wake"" on"" two-way"" fluid-structure interaction［J］." Renewable" energy" resources，" 2015，" 33（11）： 1664-1673.

［16］""" 趙元星， 汪建文， 白葉飛， 等. 切變來流下風(fēng)力機(jī)葉片應(yīng)力耦合性分析［J］. 工程熱物理學(xué)報(bào)， 2019， 40（11）： 2574-2581.

ZHAO Y X， WANG J W， BAI Y F， et al. Stress coupling analysis of wind turbine blades under wind shear flow［J］. Journal of engineering thermophysics， 2019， 40（11）： 2574-2581.

［17］""" 宋力， 黨永利， 謝曉鳳， 等. 基于流固耦合的風(fēng)力機(jī)葉片模擬分析［J］. 可再生能源， 2019， 37（10）： 1546-1550.

SONG L， DANG Y L， XIE X F， et al. Wind turbine blade simulation analysis based on fluid-solid coupling［J］. Renewable energy resources， 2019， 37（10）： 1546-1550.

［18］""" 柯世堂， 陸曼曼， 吳鴻鑫， 等. 基于風(fēng)洞試驗(yàn)15 MW風(fēng)力機(jī)葉片顫振后形態(tài)與能量圖譜研究［J］. 空氣動力學(xué)學(xué)報(bào)， 2022， 40（4）： 169-180， 168.

KE S T， LU M M， WU H X， et al. Experimental study on the post-flutter morphological chracteristics and energy dissipation"" of"" a" 15 MW" wind" turbine" blade［J］.nbsp;" Acta aerodynamica sinica， 2022， 40（4）： 169-180， 168.

［19］""" 周興. 大氣邊界層實(shí)驗(yàn)?zāi)M及水平軸風(fēng)力機(jī)葉片動力學(xué)特性實(shí)驗(yàn)研究［D］. 武漢： 華中科技大學(xué)， 2016.

ZHOU X. ABL experimental simulation and research on dynamic characteristics of horizontal axis wind turbine blade［D］. Wuhan： Huazhong University of Science and Technology， 2016.

［20］""" 潘月月， 李正農(nóng)， 徐海華， 等. 沿海灘涂風(fēng)力機(jī)近尾流風(fēng)速分布特性研究［J］. 太陽能學(xué)報(bào)， 2022， 43（9）： 193-201.

PAN Y Y， LI Z N， XU H H， et al. Near wake wind velocity distribution characteristics of wind turbine on coastal beach［J］. Acta energiae solaris sinica， 2022， 43（9）： 193-201.

［21］""" LI Z N， PU O， PAN Y Y， et al. A study on measuring wind turbine wake based on UAV anemometry system［J］. Sustainable energy technologies and assessments， 2022， 53： 102537.

［22］""" GB 50009—2012， 建筑結(jié)構(gòu)荷載規(guī)范［S］.

GB 50009—2012， Load code for the design of building structures ［S］.

［23］""" 胡丹梅， 張志超， 孫凱， 等. 風(fēng)力機(jī)葉片流固耦合計(jì)算分析［J］. 中國電機(jī)工程學(xué)報(bào)， 2013， 33（17）： 98-104， 18.

HU D M， ZHANG Z C， SUN K， et al. Computational analysis of wind turbine blades based on fluid-structure interaction［J］. Proceedings of the CSEE， 2013， 33（17）： 98-104， 18.

［24］""" 李德順， 朱瑩. 來流條件對風(fēng)力機(jī)葉片流固耦合特性影響研究［J］. 液壓氣動與密封， 2022， 42（5）： 17-24.

LI D S， ZHU Y. The effect of wind shear on wind turbine blades""""" fluid-structure""""" coupling"""" characteristics［J］. Hydraulics pneumatics amp; seals， 2022， 42（5）： 17-24.

［25］""" ZHANG W， MARKFORT C D， PORTé-AGEL F. Near-wake flow structure downwind of a wind turbine in a turbulent boundary layer［J］. Experiments in fluids， 2012， 52（5）： 1219-1235.

［26］""" 朱瑩. 來流條件對風(fēng)力機(jī)葉片流固耦合動力響應(yīng)特性的影響［D］. 蘭州： 蘭州理工大學(xué)， 2021.

ZHU Y. Fluid-structure interaction dynamic response of wind turbine with influence of inflow conditions［D］. Lanzhou： Lanzhou University of Technology， 2021.

［27］""" 劉鑫， 閆姝， 郭雨桐， 等. 復(fù)雜地形下風(fēng)力機(jī)尾流數(shù)值模擬研究和激光雷達(dá)實(shí)驗(yàn)對比［J］. 太陽能學(xué)報(bào)， 2020， 41（3）： 1-7.

LIU X， YAN S， GUO Y T， et al. Cfd simulation and lidar experimental study on wind turbines in complex terrain［J］. Acta energiae solaris sinica， 2020， 41（3）： 1-7.

［28］""" YE J J， CHU C C， CAI H， et al. A multi-scale model for studying failure mechanisms of composite wind turbine blades［J］. Composite structures， 2019， 212： 220-229.

［29］""" 李正農(nóng)， 朱勝兵， 潘月月. 考慮結(jié)冰質(zhì)量與剛度影響對風(fēng)力機(jī)葉片模態(tài)分析［J］. 玻璃鋼/復(fù)合材料， 2018（10）： 82-89， 19.

LI Z N， ZHU S B， PAN Y Y. Research about the effect of icing stiffness and mass on the mode of large wind turbine blade［J］. Fiber reinforced plastics/composites， 2018（10）： 82-89， 19.

［30］""" 趙元星. 氣動載荷影響下風(fēng)力機(jī)葉片動態(tài)應(yīng)力特性研究［D］. 呼和浩特： 內(nèi)蒙古工業(yè)大學(xué)， 2021.

ZHAO Y X. Research on dynamic stress characteristics of wind turbine blade under aerodynamic load［D］. Hohhot： Inner Mongolia University of Tehchnology， 2021.

［31］""" 姚世剛， 戴麗萍， 康順. 風(fēng)力機(jī)葉片氣動性能及流固耦合分析［J］. 工程熱物理學(xué)報(bào)， 2016， 37（5）： 988-992.

YAO S G， DAI L P， KANG S. Aerodynamic performance and fluid-structure coupling analysis of wind turbine blades［J］. Journal of engineering thermophysics， 2016， 37（5）： 988-992.

ANALYSIS OF FLUID-STRUCTURE COUPLING CHARACTERISTICS OF COASTAL WIND TURBINE BLADEDS BASED ON

MEASURED WIND FIELD

Pan Yueyue1，2，Li Zhengnong2，Zhang Yukun3，Huang Bin4

（1. College of Architectural and Civil Engineering， Weifang University， Weifang 261061， China；

2. College of Civil Engineering， Hunan University， Changsha 410082， China；

3. College of Civil Engineering and Architecture， Shandong University of Science and Technology， Qingdao 266590， China；

4. School of Civil Engineering and Architecture， Hainan University， Haikou 570228， China）

Abstract：In this paper， the wind field characteristics of a coastal wind farm are obtained by the UAVs （unmanned aerial vehicles） wind filed measurement method， and the measured wind field parameters are applied to the subsequent computational fluid dynamics simulation. Taking the full-scale model of a 2 MW wind turbine blade as the research object， a fluid-structure coupling system is established based on Workbench platform for calculation. The numerical simulation results of near-wake wind velocity distribution are compared with the measured results， and the rationality of the numerical calculation is verified. Then， the flapping deformation and equivalent stress responses of the rotating blades at rated wind velocity are calculated. Results show that the flapping deformation of wind turbine blade reaches its maximum value at the azimuth of 0°， and the maximum flapping deformation is 1916.4 mm， which is close to the results of the static loading test 2299 mm. The flapping deformation increases nonlinearly along the span-wise direction of the blade， and with a larger increasing slope at 0.60 times the radius of the wind wheel. The equivalent stress concentration area of the blade is located at the junction of the girder and the leading edge， and appears at azimuth of 60°， where the maximum value is 0.78 times of the radius of the wind wheel. Stress is asymmetrically distributed along the wind wheel height direction.

Keywords：wind turbine blades； fluid-structure interaction； wind velocity distribution； flapping deformation； equivalent stress