隨機風(fēng)沙沖蝕葉片涂層的數(shù)值模擬研究

劉士毅，岑海堂，劉光友

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

中國幅員遼闊，但其地貌復(fù)雜、各地區(qū)氣候存在很大差異。西北華北地區(qū)沙化嚴(yán)重，并形成中國八大沙漠[1-3]、東北冰雪多、東南沿海地區(qū)暴雨頻發(fā)，然而中國的風(fēng)能資源又主要聚積于此，所以風(fēng)力機的工作環(huán)境差，常年遭受風(fēng)沙沖蝕、陽光曝曬以及極端溫差作用[4]，導(dǎo)致葉片涂層損傷嚴(yán)重、修復(fù)困難、維修成本高，嚴(yán)重影響葉片的氣動性能以及風(fēng)力機組的輸出功率，極大地降低了風(fēng)力機的經(jīng)濟性[5-7]。張永等[8]基于自制試驗臺，采用氣流挾沙噴射法研究了風(fēng)沙作用下風(fēng)力機葉片涂層沖蝕特性，進而得到涂層沖蝕磨損量與沖蝕速度的關(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]針對含沙氣流對風(fēng)力機葉片沖蝕作用問題進行了數(shù)值模擬仿真，并基于實驗數(shù)據(jù)建立了玻璃鋼材料的沙蝕沖蝕率模型，進而基于風(fēng)力發(fā)電機葉輪三維流場及沙粒運動特點，計算得到了葉片表面在各工況下的沖蝕率，結(jié)果表明葉片中、外葉展前緣區(qū)域沖蝕損傷最為嚴(yán)重。本文基于已有試驗結(jié)果，利用EDEM 與Fluent 耦合進行風(fēng)沙沖蝕進行數(shù)值模擬，得到隨機載荷數(shù)據(jù)，并對葉片涂層進行有限元分析，為涂層疲勞強度計算和后期維護提供理論基礎(chǔ)。

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

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

該風(fēng)機葉片模型中，風(fēng)輪直徑D= 90 m，取輪轂半徑rhub= 1.5 m，則葉片長度L= 43.5 m，風(fēng)輪轉(zhuǎn)速n= 19.3 r/min，沿葉片展向選取6 個翼型截面，其占位分別為0.17、0.42、0.57、0.75、0.95 和1。本文結(jié)合風(fēng)力機葉片專用翼型的相關(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 隨機風(fēng)速模型建立

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

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

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

求解方程為：

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

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

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

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

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

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

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

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

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

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

Fluent 自帶的離散項也能進行數(shù)值模擬仿真，但是其限制離散項體積分?jǐn)?shù)不超過10%，只能進行稀疏的低濃度兩項流的仿真，這為工程實例的仿真帶來極大的不便，因此本文通過Fluent 與EDEM 耦合仿真，以此實現(xiàn)對風(fēng)沙兩項流沖蝕葉片涂層的模擬仿真。基于拉格朗日耦合仿真的基本思路可得到EDEM 與Fluent 耦合仿真的基本流程，如圖7 所示。

圖7 EDEM 與Fluent 耦合仿真流程

3.2 網(wǎng)格劃分

創(chuàng)建計算模型流體域，設(shè)置流體域的入口邊界inlet、出口邊界outlet、壁面wall 以及流固耦合面FSI，設(shè)置旋轉(zhuǎn)的流體進行旋轉(zhuǎn)速度與方向風(fēng)機葉片的旋轉(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 隨機風(fēng)沙沖蝕仿真

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

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

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

由于沙粒中90%以上為二氧化硅，因此直接設(shè)置沙粒材料為二氧化硅物性參數(shù)如表4 所示，風(fēng)力機葉片涂層材料為聚氨酯物性參數(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ēng)蝕仿真分析

假定接觸線周圍流場流體不可壓縮，而且湍流流動發(fā)展足夠充分，選擇標(biāo)準(zhǔn)大渦模擬模型，并設(shè)置模型參數(shù)為缺省值。采用基于壓力的求解器，流體域求解算法為SIMPLE 算法，流體域離散方法為二階迎風(fēng)離散格式[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)的仿真計算，由于隨機風(fēng)沙沖蝕葉片涂層的模擬仿真為瞬態(tài)流場問題，所以要求流場在每個時間步都需迭代至收斂狀態(tài)（Fluent默認(rèn)殘差下降到10-3為收斂）才會進行下一個時間步的求解。

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

圖9 入口速度

圖10 流固耦合面動壓

圖11 隨機沖擊載荷

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

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

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

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

通過上文對葉片表面的風(fēng)沙沖蝕涂層的隨機載荷分析，得到外載荷主要為沙粒的沖擊載荷和風(fēng)載荷。通過有限元軟件對隨機載荷作用到葉片涂層進行應(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 可知：風(fēng)沙環(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 的耦合仿真分析得到了隨機風(fēng)沙沖蝕作用下葉片涂層的載荷數(shù)據(jù)自然風(fēng)速的隨機性以及沙粒直徑的不確定性使得風(fēng)沙沖蝕葉片涂層的過程具有顯著的隨機性。

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

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