劉士毅，岑海堂，劉光友

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

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

1 風(fēng)機(jī)葉片涂層三維模型

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

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

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

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

圖1 葉片三維模型

圖2 涂層三維模型

2 隨機(jī)風(fēng)速模型建立

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

圖3 內(nèi)蒙古某大型風(fēng)電場(chǎng)實(shí)測(cè)風(fēng)速

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

求解方程為：

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

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

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

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

圖4 威布爾風(fēng)速模型的概率密度圖

圖5 威布爾風(fēng)速模型的分布圖

圖6 威布爾風(fēng)速模型的時(shí)間序列圖

3 隨機(jī)風(fēng)沙沖蝕仿真

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

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

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

圖7 EDEM 與Fluent 耦合仿真流程

3.2 網(wǎng)格劃分

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

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

3.3 隨機(jī)風(fēng)沙沖蝕仿真

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

表2 庫(kù)布齊沙漠沙粒粒徑分布及其含量分布

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

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

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

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

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

3.4 Fluent 隨機(jī)風(fēng)蝕仿真分析

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

3.5 控制方程

3.5.1 連續(xù)性方程

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

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

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

3.5.2 動(dòng)量守恒方程

動(dòng)量守恒定律是所有流動(dòng)體系都必須遵循的基本定律，即：給定流體系統(tǒng)的動(dòng)量對(duì)時(shí)間的變化率等于其所受外力的總和，根據(jù)該定律推導(dǎo)可得動(dòng)量守恒方程，也稱為運(yùn)動(dòng)方程，或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é)模型由偏微分方程組組成，因此很難得到其解析解。目前主要通過對(duì)流體域進(jìn)行離散化，進(jìn)而在離散域節(jié)點(diǎn)上建立代數(shù)方程，然后用有限差分法、有限元法及有限體積法等數(shù)值方法對(duì)所得的代數(shù)方程進(jìn)行求解。

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

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

3.6.2 控制方程的離散

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

其中：

4 仿真結(jié)果分析

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

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

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

圖9 入口速度

圖10 流固耦合面動(dòng)壓

圖11 隨機(jī)沖擊載荷

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

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

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

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

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

圖12 隨機(jī)載荷時(shí)間歷程曲線

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

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

5 結(jié)論

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

沖蝕過程中：風(fēng)速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)力小于涂層材料的屈服強(qiáng)度，所以葉片涂層受載后未發(fā)生明顯的塑性變形，因此風(fēng)沙環(huán)境下葉片涂層的失效形式表現(xiàn)為隨機(jī)脈動(dòng)應(yīng)力下的高周疲勞破壞，而且葉片涂層的疲勞損傷主要集中在葉片尖端前緣及其側(cè)面附近。

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