劉士毅，岑海堂，劉光友

（內(nèi)蒙古工業(yè)大學(xué) 機械工程學(xué)院，呼和浩特 010051）

中國幅員遼闊，但其地貌復(fù)雜、各地區(qū)氣候存在很大差異。西北華北地區(qū)沙化嚴重，并形成中國八大沙漠[1-3]、東北冰雪多、東南沿海地區(qū)暴雨頻發(fā)，然而中國的風能資源又主要聚積于此，所以風力機的工作環(huán)境差，常年遭受風沙沖蝕、陽光曝曬以及極端溫差作用[4]，導(dǎo)致葉片涂層損傷嚴重、修復(fù)困難、維修成本高，嚴重影響葉片的氣動性能以及風力機組的輸出功率，極大地降低了風力機的經(jīng)濟性[5-7]。張永等[8]基于自制試驗臺，采用氣流挾沙噴射法研究了風沙作用下風力機葉片涂層沖蝕特性，進而得到涂層沖蝕磨損量與沖蝕速度的關(guān)系，并建立了涂層沖蝕磨損程度的評價公式。徐偉勝等[9]對航空發(fā)動機葉片表面硬質(zhì)涂層受沙粒高角度沖擊而出現(xiàn)裂紋的問題進行了研究，通過實驗?zāi)M硬質(zhì)顆粒以高角度重復(fù)沖擊TiN/Ti 硬質(zhì)涂層，并研究了調(diào)制比、層數(shù)結(jié)構(gòu)對硬質(zhì)涂層沖擊損傷的影響，結(jié)果表明硬質(zhì)顆粒重復(fù)沖擊作用引起涂層疲勞剝落和圓周疲勞裂紋的產(chǎn)生。戴麗萍等[10]針對含沙氣流對風力機葉片沖蝕作用問題進行了數(shù)值模擬仿真，并基于實驗數(shù)據(jù)建立了玻璃鋼材料的沙蝕沖蝕率模型，進而基于風力發(fā)電機葉輪三維流場及沙粒運動特點，計算得到了葉片表面在各工況下的沖蝕率，結(jié)果表明葉片中、外葉展前緣區(qū)域沖蝕損傷最為嚴重。本文基于已有試驗結(jié)果，利用EDEM 與Fluent 耦合進行風沙沖蝕進行數(shù)值模擬，得到隨機載荷數(shù)據(jù)，并對葉片涂層進行有限元分析，為涂層疲勞強度計算和后期維護提供理論基礎(chǔ)。

1 風機葉片涂層三維模型

本文參照3 MW 風力機結(jié)合Glauert 法確定風機主要參數(shù)，通過氣動分析軟件Profili 和三維建模軟件UG 建立涂層三維模型。

該風機葉片模型中，風輪直徑D= 90 m，取輪轂半徑rhub= 1.5 m，則葉片長度L= 43.5 m，風輪轉(zhuǎn)速n= 19.3 r/min，沿葉片展向選取6 個翼型截面，其占位分別為0.17、0.42、0.57、0.75、0.95 和1。本文結(jié)合風力機葉片專用翼型的相關(guān)特性選用FX77-343、S818、S830 及S831 這4 種翼型，并由氣動分析軟件Profili 計算得到各翼型在最大升阻比時的攻角和升力系數(shù)，并通過葛勞渥設(shè)計法的相關(guān)公式計算得到葉片各翼型的攻角αi、入流角φi、安裝角βi及弦長ci等建模參數(shù)，如表1 所示。

表1 葉片建模參數(shù)

通過Profili 生成上述各類型翼型輪廓數(shù)據(jù)，將數(shù)據(jù)導(dǎo)入UG 建立各翼型曲線串，建立的葉片三維模型，如圖1 所示；葉片表面建立的涂層三維模型，厚度為2 mm，如圖2 所示。

圖1 葉片三維模型

圖2 涂層三維模型

2 隨機風速模型建立

本文以內(nèi)蒙古地區(qū)為例，對大型風電場風力進行風速數(shù)據(jù)監(jiān)測如圖3 所示。得到實測風速樣本總均值和總方差分別為12.34 和5.24。

圖3 內(nèi)蒙古某大型風電場實測風速

自然風受大氣壓強的影響而瞬息萬變，導(dǎo)致風速的變化具有明顯的隨機性，工程實踐表明威布爾分布是目前最能描述自然風速的數(shù)學(xué)模型，其結(jié)果與風場實際監(jiān)測數(shù)據(jù)較為吻合，被廣泛應(yīng)用于風能相關(guān)領(lǐng)域的計算中。

求解方程為：

式中：V為風速；U為[0,1]之間的隨機數(shù)。

對于給定樣本的均值和方差，求解威布爾分布的形狀參數(shù)k和尺度參數(shù)c[11]，即

最后可得到基于威布爾分布模型的風速數(shù)學(xué)為

通過MATLAB 自帶隨機函數(shù)可得威布爾風速模型概率密度如圖4 所示，威布爾風速模型的分布如圖5所示，威布爾風速模型的隨機風速時間序列如圖6 所示。

圖4 威布爾風速模型的概率密度圖

圖5 威布爾風速模型的分布圖

圖6 威布爾風速模型的時間序列圖

3 隨機風沙沖蝕仿真

由于自然風速不可避免的呈威布爾分布隨機變化，因此對于內(nèi)蒙古地區(qū)服務(wù)于沙塵環(huán)境下的風力機葉片涂層耐久性的研究，應(yīng)充分考慮變量（風速大小、沙粒直徑等）的隨機性，進而生成葉片涂層及流體域的三維網(wǎng)格模型、離散元分析軟件EDEM 計算沙粒數(shù)據(jù)以及流體分析軟件Fluent 計算流場數(shù)據(jù)，得到隨機風沙沖蝕葉片涂層的隨機載荷數(shù)據(jù)。

3.1 EDEM 和Fluent 耦合仿真的基本流程

Fluent 自帶的離散項也能進行數(shù)值模擬仿真，但是其限制離散項體積分數(shù)不超過10%，只能進行稀疏的低濃度兩項流的仿真，這為工程實例的仿真帶來極大的不便，因此本文通過Fluent 與EDEM 耦合仿真，以此實現(xiàn)對風沙兩項流沖蝕葉片涂層的模擬仿真?；诶窭嗜振詈戏抡娴幕舅悸房傻玫紼DEM 與Fluent 耦合仿真的基本流程，如圖7 所示。

圖7 EDEM 與Fluent 耦合仿真流程

3.2 網(wǎng)格劃分

創(chuàng)建計算模型流體域，設(shè)置流體域的入口邊界inlet、出口邊界outlet、壁面wall 以及流固耦合面FSI，設(shè)置旋轉(zhuǎn)的流體進行旋轉(zhuǎn)速度與方向風機葉片的旋轉(zhuǎn)方向為（0，0，1），旋轉(zhuǎn)速度為19.3 r/min。建立流體域四面體網(wǎng)格，其節(jié)點為270 136，單元數(shù)為1 154 386。網(wǎng)格劃分后的流體域模型如圖8 所示。

圖8 流體網(wǎng)格劃分

3.3 隨機風沙沖蝕仿真

本文以內(nèi)蒙古某地為例，其沙粒直徑分布及其含量分布如表2 所示。為了最大限度貼合實際工況，本次模擬仿真分析中采用的沙粒直徑數(shù)據(jù)直接參照表2 選取，結(jié)果如表3 所示。

表2 庫布齊沙漠沙粒粒徑分布及其含量分布

表3 沙粒粒徑分布及含量分布

由于沙粒中90%以上為二氧化硅，因此直接設(shè)置沙粒材料為二氧化硅物性參數(shù)如表4 所示，風力機葉片涂層材料為聚氨酯物性參數(shù)如表5 所示[12]。

表4 二氧化硅物性參數(shù)

表5 聚氨酯涂層物性參數(shù)

將涂層三維模型加載到EDEM，并按表5 設(shè)置涂層物性參數(shù)；由于不同形狀沙粒對葉片涂層最大磨損率的變化規(guī)律基本一致，將沙塵顆粒簡化為球形顆粒，忽略非球形顆粒的影響[13]。建立球形沙粒模型，通過顆粒粒徑自定義函數(shù)按表3 設(shè)置沙粒粒徑分布參數(shù)，設(shè)置接觸模型為Hertz-Mindlin，沙粒靜摩擦因數(shù)為0.5、動摩擦因數(shù)為0.15、恢復(fù)系數(shù)為0.5；沙粒靜摩擦因數(shù)為0.4、動摩擦因數(shù)為0.1、恢復(fù)系數(shù)為0.3。在流場入口建立四邊形顆粒工廠邊界，進而在邊界平面上建立動態(tài)顆粒工廠，并設(shè)置顆粒生成速率為1.2 kg/s，總質(zhì)量為5 kg。網(wǎng)格大小為最小顆粒半徑的20 倍，激活EDEM 耦合器，等待與Fluent 耦合仿真分析。

3.4 Fluent 隨機風蝕仿真分析

假定接觸線周圍流場流體不可壓縮，而且湍流流動發(fā)展足夠充分，選擇標準大渦模擬模型，并設(shè)置模型參數(shù)為缺省值。采用基于壓力的求解器，流體域求解算法為SIMPLE 算法，流體域離散方法為二階迎風離散格式[14]。

3.5 控制方程

3.5.1 連續(xù)性方程

所有流動問題都必須遵循質(zhì)量守恒定律，即：單位時間內(nèi)流體微元質(zhì)量的增量等于該時間間隔內(nèi)流入該微元體的凈質(zhì)量，由此可推出連續(xù)性方程，其微分形式為

式中： ρ為密度；t為時間；u、v、w分別為速度矢量u在x、y、z方向上的分量。

對于不可壓縮流體，密度ρ 為常數(shù)，則有

3.5.2 動量守恒方程

動量守恒定律是所有流動體系都必須遵循的基本定律，即：給定流體系統(tǒng)的動量對時間的變化率等于其所受外力的總和，根據(jù)該定律推導(dǎo)可得動量守恒方程，也稱為運動方程，或N-S 方程[15]：

式中：p為靜壓； τxx、 τyx、 τzx、 τxy、 τyy、 τzy、 τxz、 τyz及τzz為 應(yīng)力張量；Fx、Fy及Fz為體力，若體力為重力，且Z軸豎直向上，則Fx=0,Fy=0,Fz=-ρg。

3.6 流體域的離散

以上論述表明CFD 數(shù)學(xué)模型由偏微分方程組組成，因此很難得到其解析解。目前主要通過對流體域進行離散化，進而在離散域節(jié)點上建立代數(shù)方程，然后用有限差分法、有限元法及有限體積法等數(shù)值方法對所得的代數(shù)方程進行求解。

3.6.1 計算區(qū)域的離散

對于有限體積法，其計算域的離散實質(zhì)是將計算域分解為多個相互獨立的子區(qū)域，并確定各子區(qū)域的節(jié)點位置和該節(jié)點所代表的控制體積，計算域離散后可得到節(jié)點、控制體積、界面及網(wǎng)格線4 種幾何要素。

3.6.2 控制方程的離散

對于三維瞬態(tài)流場，其離散方程為

其中：

4 仿真結(jié)果分析

4.1 EDEM 和Fluent 耦合仿真結(jié)果分析

流場仿真開始后，F(xiàn)luent 會通過耦合程序關(guān)聯(lián)啟動EDEM 進行顆粒系統(tǒng)的仿真計算，由于隨機風沙沖蝕葉片涂層的模擬仿真為瞬態(tài)流場問題，所以要求流場在每個時間步都需迭代至收斂狀態(tài)（Fluent默認殘差下降到10-3為收斂）才會進行下一個時間步的求解。

整個仿真過程結(jié)束后可以得到入口速度、流固耦合面FSI 動態(tài)壓力的時間歷程數(shù)據(jù)以及沙粒沖擊葉片涂層的沖擊載荷，分別如圖9、圖10 以及圖11 所示。

圖9 入口速度

圖10 流固耦合面動壓

圖11 隨機沖擊載荷

由圖9 可知入口速度介于5~25 m/s 之間與風場實測風速數(shù)據(jù)較為吻合。

由圖10 可知流固耦合面FSI 的動壓介于20~140 Pa 之間。

由圖11 可知沙粒的沖擊載荷介于0.4~7.9 N 之間，三者均具有顯著的隨機性。

4.2 涂層仿真結(jié)果分析

通過上文對葉片表面的風沙沖蝕涂層的隨機載荷分析，得到外載荷主要為沙粒的沖擊載荷和風載荷。通過有限元軟件對隨機載荷作用到葉片涂層進行應(yīng)力響應(yīng)分析，得到葉片涂層在隨機載荷作用下應(yīng)力響應(yīng)的時間歷程數(shù)據(jù)以及葉片涂層應(yīng)力響應(yīng)最大時的應(yīng)力云圖，分別如圖12 和圖13 所示。

圖12 隨機載荷時間歷程曲線

圖13 0.86 s 時葉片涂層的應(yīng)力云圖

由圖12 可知：風沙環(huán)境下，葉片涂層的應(yīng)力響應(yīng)具有顯著的隨機性，最大應(yīng)力響應(yīng)介于0.4~2.79 MPa之間且主要集中在0.8~1.6 MPa 之間。0.86 s 時葉片涂層的應(yīng)力響應(yīng)在葉片尖端附近達到最大為2.79 MPa，如圖13 所示。

5 結(jié)論

基于EDEM 與Fluent 的耦合仿真分析得到了隨機風沙沖蝕作用下葉片涂層的載荷數(shù)據(jù)自然風速的隨機性以及沙粒直徑的不確定性使得風沙沖蝕葉片涂層的過程具有顯著的隨機性。

沖蝕過程中：風速V=13.9×[-ln(1-U)]1/2.54，葉片涂層所受的外載荷介于0.4~7.9 N 之間，葉片涂層的應(yīng)力響應(yīng)介于0.4~2.79 MPa之間且主要集中在0.8~1.6 MPa 之間。

葉片涂層的最大應(yīng)力響應(yīng)為2.79 MPa，位于葉片尖端前緣及其側(cè)面附近，由于最大應(yīng)力小于涂層材料的屈服強度，所以葉片涂層受載后未發(fā)生明顯的塑性變形，因此風沙環(huán)境下葉片涂層的失效形式表現(xiàn)為隨機脈動應(yīng)力下的高周疲勞破壞，而且葉片涂層的疲勞損傷主要集中在葉片尖端前緣及其側(cè)面附近。

根據(jù)計算結(jié)果，確定葉片涂層疲勞壽命，建立涂層精細化設(shè)計方法， 并提出預(yù)防性維修方案，提高葉片涂層防沙粒侵蝕性能。