基于機(jī)器學(xué)習(xí)的彎管固液兩相流流動(dòng)特性研究1)

肖盛鵬 朱宏博 周 岱 包 艷

(上海交通大學(xué)船舶海洋與建筑工程學(xué)院,上海 200240)

引言

管道水力輸送作為典型的顆粒-流體的固-液兩相流,是采礦和石油等工業(yè)中一種常見運(yùn)輸方式,具有運(yùn)輸距離長(zhǎng)、安全性高、運(yùn)營(yíng)和維護(hù)成本低、環(huán)境友好的優(yōu)點(diǎn)[1],此外,輸送系統(tǒng)中的直管之間通過彎曲短管進(jìn)行連接,具有很好的布局靈活性.當(dāng)固液兩相流中的顆粒進(jìn)入彎曲段,與彎曲壁面首次碰撞后,會(huì)減速并以一定角度向彎管出口偏轉(zhuǎn),在顆粒從彎曲段出口流入直管段,加速并恢復(fù)到完全發(fā)展?fàn)顟B(tài)之前,會(huì)經(jīng)歷與壁面一次或多次的碰撞,整個(gè)過程導(dǎo)致了管道壓降與壁面磨損率的增加[2].

管道壓降是管道輸運(yùn)優(yōu)化設(shè)計(jì)的主要物理量之一[3-4].管道壓降過大,輸送兩相流所需要的能量則越多,導(dǎo)致資源浪費(fèi);壓降過小,則不能保證顆粒的順利輸送,甚至導(dǎo)致管道堵塞[5].周游[6]通過對(duì)垂直管道進(jìn)行數(shù)值模擬,發(fā)現(xiàn)壓降由顆粒重力、壁面摩擦等引起.Doron 等[7]和Song 等[8]對(duì)直管道壓降進(jìn)行了測(cè)試,發(fā)現(xiàn)管道壓降大小與流體速度、管道直徑、固體進(jìn)料濃度及其他參數(shù)有關(guān).Peretz 等[9]針對(duì)90°彎管,結(jié)合文獻(xiàn)數(shù)據(jù)和數(shù)值模擬,建立了包含幾何和流動(dòng)條件的彎管壓降經(jīng)驗(yàn)方程.

壁面磨損是管道輸運(yùn)優(yōu)化設(shè)計(jì)的另外一個(gè)主要的物理量.Vieira 等[10]通過理論與實(shí)驗(yàn)相結(jié)合,提出了適用于不同流率和粒徑的壁面磨損模型.Adedeji等[11]將數(shù)值模擬結(jié)果與前人的實(shí)驗(yàn)數(shù)據(jù)進(jìn)行對(duì)比,發(fā)現(xiàn)顆粒碰撞角、速度和顆粒質(zhì)量分?jǐn)?shù)是磨損程度的主要影響因素.Chen 等[12]研究了輸送細(xì)顆粒的45°,60°和90°彎管的壁面磨損,發(fā)現(xiàn)最大磨損率明顯不同.Pei 等[13]發(fā)現(xiàn),改變彎曲段的曲率半徑直接影響顆粒的流動(dòng)磨損特性.Zhang 等[14]在計(jì)算磨損率時(shí)考慮了顆粒碰撞速度和角度并討論了模型的適應(yīng)性.

隨著計(jì)算技術(shù)的發(fā)展,數(shù)值模擬已成為預(yù)測(cè)管道水力輸送流動(dòng)特性的有效工具[15-16].對(duì)固液兩相流進(jìn)行數(shù)值模擬的方法主要分為兩種:Euler-Euler方法和Euler-Lagrange 方法[17].與Euler-Euler 方法將固液兩相均視為連續(xù)介質(zhì)相比,Euler-Lagrange 方法將流體視為連續(xù)相而將顆粒視為離散相,更好地模擬了顆粒的碰撞[18].對(duì)于運(yùn)動(dòng)中的各向異性大顆粒(>1 mm),Euler-Lagrange 方法的計(jì)算結(jié)果更符合實(shí)際情況[19].離散單元法(DEM)是Euler-Lagrange框架中模擬顆粒之間、顆粒與壁之間相互作用的最常用方法之一.在該方法中,有限數(shù)量的離散顆粒通過接觸力和非接觸力相互作用,并具有由牛頓運(yùn)動(dòng)方程描述的平移和旋轉(zhuǎn)運(yùn)動(dòng).

近年來,機(jī)器學(xué)習(xí)方法在解決工程應(yīng)用和科學(xué)領(lǐng)域的復(fù)雜問題中表現(xiàn)出了優(yōu)異性,為管道水力輸送流動(dòng)特性的預(yù)測(cè)提供了一種新的思路,目前已被應(yīng)用于多相流流型識(shí)別、壓力預(yù)測(cè)和顆粒系數(shù)等的預(yù)測(cè)[20].Liu 等[21]采用卷積神經(jīng)網(wǎng)絡(luò)來計(jì)算不同體積分?jǐn)?shù)下固液兩相流的流速,精度比現(xiàn)有流速計(jì)算方法高21%.Li 等[22]使用K近鄰回歸模型對(duì)固液兩相流的湍流波動(dòng)進(jìn)行預(yù)測(cè),顯示出很好的一致性.顏建國(guó)等[23]研究了極限學(xué)習(xí)機(jī)等模型在螺旋流動(dòng)條件下對(duì)過冷沸騰兩相流的臨界熱流密度的預(yù)測(cè)效果,比傳統(tǒng)經(jīng)驗(yàn)公式有大幅提升.賀登輝等[24]基于隨機(jī)森林模型,實(shí)現(xiàn)了對(duì)離心泵氣液兩相流的壓升問題的準(zhǔn)確預(yù)測(cè).但針對(duì)顆粒-液體兩相流的壓降和磨損問題,相關(guān)機(jī)器學(xué)習(xí)模型還鮮有應(yīng)用.

另一方面,與對(duì)管道水力輸送系統(tǒng)中的垂直、水平和傾斜直管道的廣泛研究相比[7,18,25],對(duì)彎管的研究要少得多.然而,在管道系統(tǒng)的所有部件中,彎管是最有可能出現(xiàn)問題的.與水力輸送相類似的顆粒輸送方式——?dú)饬斔鸵炎C明了這一情況,在此系統(tǒng)中,由顆粒在彎曲段碰撞后再恢復(fù)完全發(fā)展?fàn)顟B(tài)所引起的壓降增加較為明顯,在一些系統(tǒng)中占總壓降的90%[26],同時(shí),彎曲段的磨損率可能比直管的磨損率高50 倍[27].Wang 等[28]基于氣固流動(dòng)數(shù)據(jù),比較了不同機(jī)器學(xué)習(xí)模型對(duì)彎管最大磨損率的預(yù)測(cè)效果并與實(shí)驗(yàn)結(jié)果基本一致.然而,對(duì)于大顆粒水力輸送,這些問題研究較少,如果處理不好,會(huì)導(dǎo)致輸送過程中的流動(dòng)不穩(wěn)定、能耗增加和彎曲磨損破壞,甚至導(dǎo)致整個(gè)輸送過程的失敗.因此,有必要評(píng)估和預(yù)測(cè)彎曲段的壓力損失、管道壁面磨損來確保輸送穩(wěn)定性,以便更好地設(shè)計(jì)和規(guī)劃輸送系統(tǒng).

本文采用CFD-DEM 方法,計(jì)算了彎曲段水力提升管道在不同的幾何、流場(chǎng)和顆粒參數(shù)條件下的壓降和磨損性質(zhì),分析了各參數(shù)的影響機(jī)理,并結(jié)合機(jī)器學(xué)習(xí)方法,基于CFD-DEM 方法計(jì)算的數(shù)據(jù),研究了不同機(jī)器學(xué)習(xí)模型對(duì)彎曲提升管的壓降和磨損率的預(yù)測(cè)能力,以期為管道系統(tǒng)的設(shè)計(jì)和優(yōu)化提供參考.

1 兩相流模型

1.1 流體控制方程

將流體視為連續(xù)相,其控制方程以N-S 方程為基礎(chǔ),在動(dòng)量方程中引入了由于顆粒存在而引起的流體擾動(dòng)項(xiàng).對(duì)于不可壓縮黏性流體,方程為

其中,ρf是流體密度,uf是流體速度,pf是流體壓力,g是重力加速度,fsf是顆粒對(duì)流體的體積力,Rf是應(yīng)力張量,包括黏性應(yīng)力和雷諾應(yīng)力,εf為流體體積分?jǐn)?shù),表示為

其中Vi,part是顆?？傮w積,Vcell是網(wǎng)格體積,n是與網(wǎng)格體積重疊的顆粒數(shù).

1.2 顆粒控制方程

在DEM 模型中,顆粒的平移和旋轉(zhuǎn)運(yùn)動(dòng)分別通過動(dòng)量和角動(dòng)量方程來考慮

其中mp是顆粒質(zhì)量,up是顆粒速度,ffs是流體對(duì)顆粒的作用力

其中fd是曳力[29],fl是升力[30],fp是壓力梯度力,fam是附加質(zhì)量力.fc表示顆粒之間的接觸力,采用Hertz-Mindlin 非滑動(dòng)接觸理論模型.Ip為質(zhì)點(diǎn)轉(zhuǎn)動(dòng)慣量,ωp為顆粒的角速度,Mct為顆粒與其他所有顆粒的接觸力矩,Mfs為顆粒受流體作用的力矩.

1.3 磨損率方程

采用Zhang 等[31]設(shè)計(jì)的E/CRC 模型

式中,ER為磨損率,定義為單位時(shí)間和面積下管壁因顆粒碰撞而損失的質(zhì)量;Aface為管道系統(tǒng)內(nèi)壁總面積;BH為壁面材料的布氏硬度;Fs為顆粒形狀系數(shù);Vp為顆粒撞擊速度;θ是以弧度為單位的撞擊角;n與C為經(jīng)驗(yàn)常數(shù),分別為2.17×10-7和2.41.

1.4 CFD-DEM 耦合過程

首先計(jì)算得到流場(chǎng)參數(shù),如流場(chǎng)壓力和速度,用于DEM 計(jì)算,計(jì)算顆粒所受的流場(chǎng)力和顆粒對(duì)流場(chǎng)的反作用力,之后將顆粒對(duì)流場(chǎng)的作用力代入流體控制方程,得到相應(yīng)的流場(chǎng)信息.最后,將所有信息作為下一次迭代的參數(shù),重復(fù)迭代直至收斂,計(jì)算完成.

2 數(shù)值模擬設(shè)計(jì)與驗(yàn)證

針對(duì)由水平直管道、彎曲管道和提升直管道組成的組合管道系統(tǒng),研究其在不同的管道彎曲角度Ab、管道彎曲半徑Rb、輸入速度U、顆粒濃度Cp、顆粒直徑dp條件下的壓降ΔP和磨損率ER變化規(guī)律.假設(shè)管道變形極小,并簡(jiǎn)化為剛性管道.水平直管道和提升直管道的設(shè)置是為使進(jìn)出彎曲管道的兩相流發(fā)展到穩(wěn)定狀態(tài),以便更完整準(zhǔn)確地模擬兩相流流經(jīng)彎管的過程.

2.1 邊界條件和參數(shù)設(shè)置

如圖1 所示,水平段和提升段均長(zhǎng)1 m(獨(dú)立性驗(yàn)證見后),管徑D為3.06 cm,管道兩端設(shè)置速度入口和壓力出口邊界條件,壁面采用無滑移邊界條件.管道劃分成六面體網(wǎng)格,通過O 型網(wǎng)格拓?fù)浣Y(jié)構(gòu),可以細(xì)化壁面附近的網(wǎng)格來更好地預(yù)測(cè)壁面附近的液體流動(dòng).邊界層共有9 層,伸長(zhǎng)率為1.2,第一層厚度為5.0×10-4m,網(wǎng)格數(shù)量為1.25×105左右.流體的密度與動(dòng)力黏度分別為998 kg/m3和0.001 Pa·s,顆粒的楊氏模量、泊松比和密度分別為1.0×107Pa,0.3 和2450 kg/m3.顆粒彼此間碰撞恢復(fù)系數(shù)以及顆粒與壁面間的碰撞恢復(fù)系數(shù)分別為0.85 和0.95;顆粒間的滑動(dòng)摩擦系數(shù)和靜摩擦系數(shù)分別為0.01 和0.1;顆粒與壁面間的滑動(dòng)摩擦系數(shù)和靜摩擦系數(shù)分別為0.01 和0.2.

圖1 模擬彎管的幾何和網(wǎng)格示意Fig.1 Geometry and mesh representation of the simulated bend

2.2 數(shù)值格式

控制方程采用有限體積法離散.瞬態(tài)問題時(shí)間離散采用二階隱式格式,空間離散采用二階迎風(fēng)格式.液相模擬采用可實(shí)現(xiàn)的k-ε湍流模型.在CFDDEM 模擬中,顆粒相的時(shí)間步長(zhǎng)受Rayleigh 時(shí)間?tp的限制,Li 等[32]計(jì)算的Rayleigh 時(shí)間為

其中,μp,Ep分別為顆粒的泊松比和楊氏模量.考慮準(zhǔn)確率和時(shí)間代價(jià),顆粒時(shí)間步長(zhǎng)設(shè)定為1.0×10-5s.根據(jù)Tsuji 等[33]的工作,流體時(shí)間步長(zhǎng)可以是固體時(shí)間步長(zhǎng)的10～100 倍.為減少模擬時(shí)間,將液相時(shí)間步長(zhǎng)定為1.0×10-4s.

2.3 模型驗(yàn)證

Toda 等[34]針對(duì)此類組合管進(jìn)行了純流體與流體-顆粒兩相流的實(shí)驗(yàn).實(shí)驗(yàn)的顆粒密度為2.5 g/cm3,彎曲半徑為12 cm,管道直徑為3.02 cm,輸入速度為2.18 m/s.實(shí)驗(yàn)取水平管中距離彎管入口0.69,0.48,0.3 和0.09 m 的橫截面以及提升管中距離彎管出口0.12,0.33,0.51 和0.84 m 的橫截面,8 個(gè)截面之間的管道依次為1～7 段,第4 段主要為彎管,對(duì)每一段進(jìn)行了壓降測(cè)量.模擬結(jié)果與實(shí)驗(yàn)結(jié)果相似,每段的壓降見圖2,證明了壓降模擬的準(zhǔn)確性,也說明了彎曲段的壓降顯著大于直管段,在管道系統(tǒng)總壓降中占比較大.

圖2 壓降的模擬結(jié)果與實(shí)驗(yàn)結(jié)果驗(yàn)證Fig.2 Verification of simulated (SL) and experimental (EXP) results of pressure drop

利用E/CRC 磨損模型,Zhang 等[31]進(jìn)行了水-硅粉流的壁面磨損模擬,Chen 等[12]進(jìn)行了管道直徑為4 cm、彎曲半徑為6 cm、輸送速度為3 m/s 的水-砂流管道磨損模擬,均與實(shí)驗(yàn)測(cè)量的結(jié)果基本一致,證明了在CFD-DEM 框架中應(yīng)用E/CRC 磨損模型的計(jì)算準(zhǔn)確性,本文不再驗(yàn)證.

2.4 獨(dú)立性驗(yàn)證

網(wǎng)格獨(dú)立性依次對(duì)細(xì)、中和粗網(wǎng)格進(jìn)行驗(yàn)證,模擬結(jié)果相似,見圖3(a).考慮到計(jì)算精度和計(jì)算時(shí)間,采用了中網(wǎng)格.計(jì)算時(shí)間獨(dú)立性取每0.5 s 結(jié)果的均值,見圖3(b),從第2 s 以后,每0.5 s 的均值基本相同,認(rèn)為計(jì)算收斂,因此計(jì)算時(shí)間取3.5 s.

圖3 計(jì)算結(jié)果隨(a)網(wǎng)格數(shù)量和(b)計(jì)算時(shí)間的變化曲線Fig.3 Curve of the change in calculation results with (a) the number of grids and (b) computing time

管道長(zhǎng)度獨(dú)立性驗(yàn)證部分,對(duì)水平段和提升段均為2 m 的同類管道系統(tǒng)進(jìn)行模擬,得到管道各橫截面的顆粒徑向體積分?jǐn)?shù)分布,如圖4,可以發(fā)現(xiàn)在水平段和提升段內(nèi)經(jīng)過約0.75 m 后兩相流均發(fā)展到穩(wěn)定狀態(tài).因此,水平段和提升段總計(jì)算域均取1 m.水平段的前0.75 m 作為兩相流發(fā)展段不計(jì)入組合管道的壓降和磨損率測(cè)量,后0.25 m 作為組合管道的水平直管部分計(jì)入測(cè)量,提升直管部分則全部測(cè)量.

圖4 長(zhǎng)度獨(dú)立性驗(yàn)證Fig.4 Length independence verification

3 數(shù)值模擬結(jié)果分析

管道壓降與水力輸送的能量消耗有關(guān),由顆粒與顆粒、顆粒與流體、顆粒與壁面、流體與壁面的相互作用以及顆粒重力等的耗能導(dǎo)致[6],管道磨損率則主要由于顆粒碰撞和沖擊導(dǎo)致,兩者是工程實(shí)際應(yīng)用中的主要關(guān)注點(diǎn).下面展開管道彎曲角度、管道彎曲半徑、輸入速度、顆粒濃度和顆粒直徑這5 個(gè)因素對(duì)壓降和磨損率的影響分析.

3.1 輸入速度U

管道系統(tǒng)的彎曲角度為90°,彎曲半徑為0.06 m,粒徑0.003 m,顆粒濃度4%,輸送速度分別為2,3.2,4.4 和5.6 m/s.圖5 顯示了顆粒的瞬時(shí)流動(dòng)速度,可劃分為3 個(gè)過程:入口直管段充分發(fā)展?fàn)顟B(tài),彎曲段的碰撞減速狀態(tài)以及出口直管段加速至恢復(fù)充分發(fā)展?fàn)顟B(tài).如圖6 所示,隨著輸入速度增加,管道系統(tǒng)的壓降和磨損率增大.

圖5 不同輸入速度的顆粒流速Fig.5 Particle velocity at different input velocity

圖6 壓降、磨損率隨輸入速度的變化曲線Fig.6 Curve of pressure drop and erosion rate with input velocty

當(dāng)輸入速度增加,顆粒以更大的速度與彎曲壁面碰撞,對(duì)壁面的磨損程度增加,與1.3 節(jié)磨損率模型中磨損率與顆粒撞擊速度的關(guān)系相符合,而且顆粒與壁面的碰撞耗能增大;同時(shí),從圖5 可以看出碰撞后顆粒的速度減幅變大,動(dòng)能損失增大,且隨著輸入速度增大,顆粒間的相互碰撞更加頻繁和劇烈,碰撞后的顆粒一部分向管道壁面運(yùn)動(dòng),導(dǎo)致顆粒與壁面的碰撞次數(shù)也增多,顆粒與壁面的每次碰撞都是一次減速再加速的耗能過程,導(dǎo)致壁面磨損與輸送耗能增大;由于輸送速度的增加,管道每秒輸送的顆粒增多,導(dǎo)致顆?？傊亓Φ脑黾?流體需提供更多能量給顆粒,增大耗能;隨速度的增大,液體與壁面的相互作用也增加.管道系統(tǒng)總能耗增大,因此壓降增加,且磨損率增大.

3.2 顆粒濃度Cp

管道系統(tǒng)的彎曲角度為90°,彎曲半徑為0.06 m,粒徑0.003 mm,顆粒濃度分別為1%,2%,4%和8%,輸送速度為3.2 m/s,顆粒瞬時(shí)流動(dòng)速度見圖7.由圖8 可知,隨著輸送濃度增大,管道系統(tǒng)的壓降和磨損率增大.

圖7 不同輸送濃度的顆粒流速Fig.7 Particle velocity at different particle concentration

圖8 壓降、磨損率隨顆粒濃度的變化曲線Fig.8 Curve of pressure drop and erosion rate with particle concentration

當(dāng)輸送濃度增大,輸送的顆粒數(shù)量增多,流體需要轉(zhuǎn)遞更多的能量到顆粒以克服顆粒重力,同時(shí),管道內(nèi)顆粒數(shù)量的增多導(dǎo)致顆粒與顆粒以及顆粒與壁面的碰撞次數(shù)增加,顆粒產(chǎn)生的變速運(yùn)動(dòng)以及旋轉(zhuǎn)同時(shí)也對(duì)流場(chǎng)產(chǎn)生極大擾動(dòng),顆粒與流體間相互作用增大.以上因素導(dǎo)致了耗能與磨損,因此壓降增加,且管壁的磨損率增加.

3.3 顆粒直徑dp

管道系統(tǒng)的彎曲角度為90°,彎曲半徑為0.06 m,粒徑分別為0.0015,0.003,0.0045 和0.006 m,濃度為4%,輸入速度為3.2 m/s,顆粒瞬時(shí)流動(dòng)速度見圖9.圖10 表明,隨著粒徑增大,管道系統(tǒng)的壓降和磨損率增大.

圖9 不同顆粒直徑的瞬時(shí)顆粒流速Fig.9 Particle velocity at different particle diameter

圖10 壓降、磨損率隨顆粒直徑的變化曲線Fig.10 Curve of pressure drop and erosion rate with particle diameter

相同濃度下,隨著粒徑增大,顆粒數(shù)量急劇減少,輸送過程中顆粒間的相互碰撞減少,因此顆粒的動(dòng)量損失更少,會(huì)以更大的速度與彎曲壁面相碰撞;且大顆粒與小顆粒相比,貼近彎曲壁面的一層顆粒之間距離更近,具有一定的“屏障”作用,減少了非最外層顆粒與彎曲壁面的碰撞,如圖9(a)和圖9(b)中的彎曲段較內(nèi)側(cè)的顆粒比貼近壁面?zhèn)阮w粒的速度減幅更小,因此大顆粒比小顆粒對(duì)壁面的磨損率要大.對(duì)于大顆粒,在顆粒重力不變的情況下,雖然顆粒間的相互碰撞減少導(dǎo)致耗能減少,但由于顆粒與壁面碰撞的速度更大次數(shù)更多,導(dǎo)致總能耗增大,壓降增大.

3.4 彎曲半徑Rb

管道系統(tǒng)的彎曲角度為90°,彎曲半徑分別為0.03,0.06,0.09 和0.12 m,粒徑0.003 mm,濃度為4%,輸入速度為3.2 m/s,顆粒瞬時(shí)流動(dòng)速度見圖11.

圖11 不同彎曲半徑的瞬時(shí)顆粒流速Fig.11 Particle velocity at different bending radius

由圖11 可以看出,在顆粒流經(jīng)彎曲段的過程中,由于離心力的作用,顆粒呈現(xiàn)貼近外側(cè)壁面的流動(dòng)狀態(tài)并沿外側(cè)壁面分散.因碰撞外側(cè)壁面而減速形成低速流,且速度減小幅度隨彎曲半徑的增大而減小,即彎曲段半徑越大越平緩,其中的顆粒流速越接近完全發(fā)展?fàn)顟B(tài)的速度,顆粒以更大的速度與彎曲壁面發(fā)生碰撞,且彎曲段距離隨彎曲半徑的增大而增大導(dǎo)致顆粒貼近外側(cè)壁面的流動(dòng)距離增大,碰撞范圍更大,碰撞和摩擦次數(shù)增加,導(dǎo)致磨損率增加.隨彎曲半徑的增大,顆粒在彎曲段的速度減幅變小,動(dòng)能損失減小,與顆粒和壁面間的相互作用增大導(dǎo)致的耗能的增量基本相抵消,因此不同彎曲半徑的管道系統(tǒng)壓降變化不大.

由圖12 可以看出,隨彎曲半徑的增大,壓降基本不變,磨損率增大.

圖12 壓降、磨損率隨彎曲半徑的變化曲線Fig.12 Curve of pressure drop and erosion rate with bending radius

3.5 彎曲角度Ab

管道系統(tǒng)的彎曲角度分別為30°,60°,90°,120°和150°,彎曲半徑為0.06 m,粒徑0.003 mm,濃度4%,輸入速度3.2 m/s,顆粒瞬時(shí)流動(dòng)速度見圖13.

圖13 不同彎曲角度的瞬時(shí)顆粒流速Fig.13 Particle velocity at different bending angle

如圖13 所示,隨著彎曲角度從150°減小到90°,彎折程度趨于平緩,彎曲段的長(zhǎng)度減小,顆粒在彎曲段貼近壁面流動(dòng)的距離減小,與壁面的碰撞和摩擦減少,磨損率降低,同時(shí)造成了耗能降低,導(dǎo)致壓降的下降.隨著彎曲角度從90°減小到30°,彎折程度更加平緩,接近于水平直管,彎曲形狀對(duì)顆粒的影響減小,顆粒貼近彎曲壁面的程度降低而更加均勻地分散在管道中,顆粒間的相互作用減少,耗能減少,壓降變小;顆粒速度減幅有所降低,以相對(duì)大一些的速度與壁面碰撞,導(dǎo)致顆粒與壁面的磨損率小幅上升.

如圖14 所示,隨著彎曲角度從150°減小到90°,壓降和磨損率下降;隨著彎曲角度從90°減小到30°,壓降下降,磨損率小幅增加.

圖14 壓降、磨損率隨彎曲角度的變化曲線Fig.14 Curve of pressure drop and erosion rate with bending angle

4 機(jī)器學(xué)習(xí)數(shù)據(jù)集及模型

由于固液兩相流不同條件下流動(dòng)狀態(tài)的復(fù)雜性,以及彎管幾何特征的多樣性,導(dǎo)致多參數(shù)條件下彎管的壓降和磨損率預(yù)測(cè)難度較大,基于CFDDEM 的方法將耗費(fèi)大量的計(jì)算時(shí)間和資源.因此,采用機(jī)器學(xué)習(xí)方法,通過對(duì)已有數(shù)據(jù)的學(xué)習(xí),來實(shí)現(xiàn)良好預(yù)測(cè)效果.

4.1 機(jī)器學(xué)習(xí)數(shù)據(jù)集

表1 列出了前文彎曲角度、彎曲半徑、顆粒輸送濃度、顆粒直徑和輸入速度5 個(gè)參數(shù)的常見取值.將參數(shù)組合為工況時(shí),若各參數(shù)一一組合,工況數(shù)量將十分龐大,故采用Pairwise 配對(duì)法,它是對(duì)正交分析方法優(yōu)化后得到的方法,在保證較大測(cè)試覆蓋度的情況下,縮減測(cè)試用例,降低算例耗時(shí).本文在完成數(shù)值模擬驗(yàn)證的基礎(chǔ)上,對(duì)5 個(gè)參數(shù)的不同取值進(jìn)行組合,共計(jì)255 個(gè)工況.經(jīng)CFD-DEM 數(shù)值模擬計(jì)算后,檢查并剔除計(jì)算異常的數(shù)據(jù),剩下234 個(gè)有效數(shù)據(jù),組成了最終的數(shù)據(jù)集.

表1 工況的各參數(shù)取值Table 1 Parameter values for cases

使用Python 中sklearn 庫將此數(shù)據(jù)集按照80%和20%的比例隨機(jī)劃分為訓(xùn)練集和測(cè)試集,并以訓(xùn)練集為標(biāo)準(zhǔn)對(duì)訓(xùn)練集和測(cè)試集進(jìn)行歸一化,以提高數(shù)據(jù)的表現(xiàn)

其中xmax和xmin數(shù)據(jù)所在特征列的最大值和最小值,經(jīng)歸一化后x的范圍為[0,1].利用訓(xùn)練集來訓(xùn)練機(jī)器學(xué)習(xí)模型,并通過測(cè)試集來對(duì)模型的準(zhǔn)確性進(jìn)行評(píng)估.

4.2 機(jī)器學(xué)習(xí)模型

基于4.1 節(jié)的數(shù)據(jù)集,總共6 個(gè)模型被開發(fā)和評(píng)估.接下來對(duì)每個(gè)模型進(jìn)行簡(jiǎn)要介紹.

4.2.1 線性回歸(LR)

線性回歸是基于輸入特征和輸出是線性關(guān)系的假設(shè),進(jìn)而求解線性方程.當(dāng)數(shù)據(jù)量較小時(shí),能達(dá)到與復(fù)雜模型相當(dāng)?shù)木?模型如下

4.2.2K近鄰回歸(KNN)

K近鄰回歸是回歸分析中經(jīng)典的非參數(shù)方法.在給定輸入變量x的情況下,K近鄰算法使用輸入空間中最接近x的K個(gè)點(diǎn)的均值來確定輸出值.模型的預(yù)測(cè)效果受K值的影響.

4.2.3 支持向量回歸(SVM)

支持向量回歸基于從訓(xùn)練數(shù)據(jù)確定的一維或多維超平面對(duì)數(shù)據(jù)進(jìn)行擬合,使靠超平面最遠(yuǎn)的樣本點(diǎn)之間的間隔最大,并限制超平面與真實(shí)值的偏差必須小于等于ε.同時(shí),由于無法保證所有點(diǎn)都在ε范圍內(nèi),故對(duì)每個(gè)樣本設(shè)置超出范圍上下界的松弛變量ξ+和ξ-,并乘正則化參數(shù)C.模型如下

SVM 的預(yù)測(cè)效果主要受到正則化參數(shù)C和范圍ε這兩個(gè)超參數(shù)的影響.

4.2.4 人工神經(jīng)網(wǎng)絡(luò)(ANN)

人工神經(jīng)網(wǎng)絡(luò)由多個(gè)人工神經(jīng)元組成,這些神經(jīng)元以一定權(quán)重相連組成神經(jīng)網(wǎng)絡(luò),并分為輸入層、隱藏層和輸出層.整個(gè)過程為:開始時(shí)為連接輸入層、隱藏層和輸出層中神經(jīng)元的權(quán)重wij分配隨機(jī)值,輸入值被饋送到輸入層中的神經(jīng)元中,乘以權(quán)重wij并加偏置bj,再經(jīng)非線性激活函數(shù)的處理,輸出到隱藏層的神經(jīng)元中,并經(jīng)過相同過程繼續(xù)輸出和傳遞,最終傳遞到輸出層輸出.每個(gè)神經(jīng)元的計(jì)算過程為

其中f為非線性激活函數(shù).將輸出層輸出的預(yù)測(cè)值與真實(shí)值進(jìn)行比較,兩者之間的差值用于通過反向傳播算法調(diào)整權(quán)重.重復(fù)前饋和反向傳播算法以迭代地調(diào)整權(quán)重,直到真實(shí)值和預(yù)測(cè)值之間的誤差達(dá)到可接受的范圍.模型的預(yù)測(cè)效果受到網(wǎng)絡(luò)結(jié)構(gòu)、學(xué)習(xí)率等較多超參數(shù)的影響.

4.2.5 隨機(jī)森林(RF)

隨機(jī)森林是由許多決策樹組成的集成模型,即采用不同的隨機(jī)選擇的子集和特征屬性建立多個(gè)獨(dú)立的決策樹,最后將多個(gè)決策樹的結(jié)果進(jìn)行平均或加權(quán)平均,合并成一個(gè)效果更佳的集成學(xué)習(xí)的模型

式中,Ti代表每顆決策樹的結(jié)果,xinput為決策每棵樹采用的子集數(shù)據(jù),ntree為決策樹的數(shù)量.模型的預(yù)測(cè)能力主要與決策樹的最大深度和數(shù)量有關(guān).

4.2.6 XGBoost

XGBoost 屬于集成學(xué)習(xí)中的boosting 框架里的算法,采用多個(gè)基學(xué)習(xí)器,即不斷生成新的基學(xué)習(xí)器,每個(gè)基學(xué)習(xí)器都是基于前面基學(xué)習(xí)器和目標(biāo)值的差值來學(xué)習(xí),從而不斷降低模型的偏差,最終模型的預(yù)測(cè)結(jié)果是所有基學(xué)習(xí)器預(yù)測(cè)結(jié)果的加和

4.2.7 模型評(píng)價(jià)指標(biāo)

為了量化機(jī)器學(xué)習(xí)模型的預(yù)測(cè)準(zhǔn)確性,使用了3 個(gè)評(píng)價(jià)指標(biāo),即決定系數(shù)(R2)(范圍在0～1 之間)、均方根誤差(RMSE)和平均絕對(duì)誤差(MAE)

其中yi表示真實(shí)值,表示預(yù)測(cè)值,表示各真實(shí)值的平均值,n是數(shù)據(jù)的數(shù)量.

R2表示模型能夠解釋數(shù)據(jù)方差的比例,通常用于比較不同模型的表現(xiàn),MAE計(jì)算了預(yù)測(cè)值和真實(shí)值之間差異的確切幅度,RMSE計(jì)算了誤差的平均大小,并且對(duì)異常值敏感.越接近1 的R2值和越低的RMSE和MAE值表明模型準(zhǔn)確性越高.

5 機(jī)器學(xué)習(xí)模型預(yù)測(cè)效果及討論

5.1 特征相關(guān)性分析

圖15 給出了各特征之間的相關(guān)性,可以初步看出,5 個(gè)輸入特征相互之間的相關(guān)性均在0.1 以下,相關(guān)性低,獨(dú)立性高,較好地表示管道固液兩相流的不同特征.輸出特征壓降與彎曲角度、輸入速度、顆粒濃度和顆粒直徑呈正相關(guān),其中與輸送速度正相關(guān)性程度最高,與彎曲半徑呈弱負(fù)相關(guān);磨損率與彎曲角度呈弱負(fù)相關(guān),與彎曲半徑、輸入速度、顆粒濃度和顆粒直徑呈正相關(guān);兩個(gè)輸出特征壓降與磨損率之間也存在一定的相關(guān)性,即磨損率的增加代表著顆粒與壁面的相互作用增大,產(chǎn)生的能量損耗增大,一定程度導(dǎo)致壓降增大.各特征的相關(guān)性與前文第3 節(jié)分析的規(guī)律大致吻合.為了便于后續(xù)定量計(jì)算分析輸入特征對(duì)于壓降和磨損率各自的影響程度大小,對(duì)這兩個(gè)物理量分別建立機(jī)器學(xué)習(xí)模型并預(yù)測(cè).

圖15 特征相關(guān)性Fig.15 Feature correlation

5.2 各機(jī)器學(xué)習(xí)模型對(duì)壓降的預(yù)測(cè)效果

通過Python 的sklearn 庫以及torch 庫建立機(jī)器學(xué)習(xí)模型,并對(duì)模型中由用戶自行設(shè)置大小的超參數(shù)采用隨機(jī)搜索方法進(jìn)行自動(dòng)尋優(yōu),以得到最好的模型預(yù)測(cè)效果.

預(yù)測(cè)壓降的各模型超參數(shù)設(shè)置和在測(cè)試集的預(yù)測(cè)結(jié)果見表2,以R2為準(zhǔn)確率的主要評(píng)估指標(biāo)進(jìn)行排序,可以發(fā)現(xiàn),神經(jīng)網(wǎng)絡(luò)、隨機(jī)森林和XGBoost 的準(zhǔn)確率均為0.96 左右,其中以神經(jīng)網(wǎng)絡(luò)的0.9654 為最高,同時(shí)3 個(gè)模型的MAE和RMSE也較為相近,以神經(jīng)網(wǎng)絡(luò)MAE=0.0325 和RMSE=0.0472 為最低,3 個(gè)模型均有較好的預(yù)測(cè)效果,在測(cè)試集的預(yù)測(cè)表現(xiàn)見圖16.預(yù)測(cè)效果最差的為K近鄰模型,R2僅為0.8212.

表2 預(yù)測(cè)壓降的各模型超參數(shù)及預(yù)測(cè)效果Table 2 Hyper parameters and prediction effects of various models for pressure drop

圖16 各模型對(duì)壓降的預(yù)測(cè)效果Fig.16 Prediction effect of various models on pressure drop

5.3 各機(jī)器學(xué)習(xí)模型對(duì)磨損率的預(yù)測(cè)效果

預(yù)測(cè)磨損率的各模型超參數(shù)設(shè)置和在測(cè)試集的預(yù)測(cè)結(jié)果見表3,可以發(fā)現(xiàn),神經(jīng)網(wǎng)絡(luò)、隨機(jī)森林和XGBoost 的準(zhǔn)確率同樣較為相近,其中以神經(jīng)網(wǎng)絡(luò)的0.9907 為最高,同時(shí)3 個(gè)模型的MAE和RMSE也較為相近,以XGBoost 的MAE=0.0086 和神經(jīng)網(wǎng)絡(luò)RMSE=0.0126 為最低,3 個(gè)模型仍均有良好的預(yù)測(cè)效果,在測(cè)試集的預(yù)測(cè)表現(xiàn)見圖17.預(yù)測(cè)效果最差的為線性回歸模型,R2僅為0.6184.

5.4 數(shù)據(jù)量大小與預(yù)測(cè)準(zhǔn)確率的關(guān)系

采用5.2 節(jié)和5.3 節(jié)中對(duì)壓降和磨損率預(yù)測(cè)準(zhǔn)確率最高的神經(jīng)網(wǎng)絡(luò)模型,分別用數(shù)據(jù)量為50,100,150,200 和234 的數(shù)據(jù)集進(jìn)行訓(xùn)練和評(píng)估.由表4可知,對(duì)于壓降,當(dāng)數(shù)據(jù)集從50 增長(zhǎng)到150,預(yù)測(cè)準(zhǔn)確率總體呈增長(zhǎng)的趨勢(shì),從150 增長(zhǎng)到234 后,預(yù)測(cè)準(zhǔn)確率基本持平,R2在0.96 左右;對(duì)于磨損率,則更為明顯,當(dāng)數(shù)據(jù)集從50 增長(zhǎng)到150,預(yù)測(cè)準(zhǔn)確率顯著增長(zhǎng),從150 增長(zhǎng)到234 后,預(yù)測(cè)準(zhǔn)確率基本持平,R2在0.987 左右.說明數(shù)據(jù)量已經(jīng)基本達(dá)到機(jī)器學(xué)習(xí)模型的要求,從計(jì)算成本的角度出發(fā),234 個(gè)數(shù)據(jù)的量可以實(shí)現(xiàn)模型的較高準(zhǔn)確性.

表4 數(shù)據(jù)集大小與預(yù)測(cè)準(zhǔn)確性的關(guān)系Table 4 The relationship between data size and prediction accuracy

5.5 特征相對(duì)重要性分析

機(jī)器學(xué)習(xí)模型可以計(jì)算每個(gè)輸入特征的相對(duì)重要性,以衡量每個(gè)特征對(duì)構(gòu)建預(yù)測(cè)模型的相關(guān)性和貢獻(xiàn)程度.選取對(duì)壓降和磨損率預(yù)測(cè)R2分別最高的神經(jīng)網(wǎng)絡(luò)模型,計(jì)算各特征的相對(duì)重要性,結(jié)果見表5.可以看出,輸入速度對(duì)于壓降的預(yù)測(cè)影響最大,彎曲半徑對(duì)壓降的預(yù)測(cè)影響最小;顆粒濃度對(duì)磨損率的預(yù)測(cè)影響最大,彎曲角度對(duì)磨損率的預(yù)測(cè)影響最小,為彎管固液兩相流的進(jìn)一步研究指明了方向.

表5 各特征對(duì)于壓降和磨損率預(yù)測(cè)的相對(duì)重要性Table 5 Relative importance of each feature for prediction of pressure drop and erosion rate

6 結(jié)論

通過CFD-DEM 耦合方法和機(jī)器學(xué)習(xí)方法,針對(duì)彎管內(nèi)顆粒水力輸送的壓降和磨損率問題,探討了彎曲角度、彎曲半徑、輸入速度、顆粒濃度和顆粒直徑這5 個(gè)特征因素的影響并對(duì)結(jié)果進(jìn)行預(yù)測(cè),主要結(jié)論如下.

(1)基于CFD-DEM 數(shù)值模擬計(jì)算的結(jié)果,發(fā)現(xiàn)彎管的壓降隨著輸入速度、顆粒濃度、顆粒直徑和彎曲角度的增大而增大,受彎曲半徑的影響較小;彎管的磨損率隨輸入速度、顆粒濃度、顆粒直徑和彎曲半徑的增大而增大,隨彎曲角度的增大,在90°前先略有下降,在90°后增大.

(2)根據(jù)Pairwise 配對(duì)法確定了工況并進(jìn)行數(shù)值模擬計(jì)算,最終建立了包含數(shù)百個(gè)數(shù)據(jù)的數(shù)據(jù)集.基于此數(shù)據(jù)集訓(xùn)練了6 個(gè)機(jī)器學(xué)習(xí)模型,分別對(duì)彎管的壓降和磨損率進(jìn)行預(yù)測(cè).其中,神經(jīng)網(wǎng)絡(luò)、隨機(jī)森林和XGBoost 3 個(gè)模型的預(yù)測(cè)效果最為良好且相近,對(duì)壓降和磨損率的R2表現(xiàn)分別在0.96 和0.99 左右.并發(fā)現(xiàn)輸入速度和顆粒濃度分別是對(duì)壓降和磨損率的預(yù)測(cè)影響程度最大的因素.

相比于對(duì)磨損率0.99 左右的準(zhǔn)確率,機(jī)器學(xué)習(xí)模型對(duì)壓降的預(yù)測(cè)準(zhǔn)確率仍有提升空間,因此下一步的一項(xiàng)工作是將目前效果良好且相近的多個(gè)模型進(jìn)行組合,形成集成學(xué)習(xí)模型,以進(jìn)一步提高機(jī)器學(xué)習(xí)模型的準(zhǔn)確性;另一方面,針對(duì)更加完整和復(fù)雜的輸流管路,增加對(duì)管道彎曲的垂直或水平朝向、管道直徑等因素的考慮,計(jì)算更多工況以訓(xùn)練適應(yīng)度更高的機(jī)器學(xué)習(xí)模型.