基于離散單元法和人工神經(jīng)網(wǎng)絡(luò)的近壁顆粒動(dòng)力學(xué)特征研究1)

段總樣 趙云華 徐 璋

(浙江工業(yè)大學(xué)機(jī)械工程學(xué)院,杭州 310014)

引言

密集顆粒流及氣粒兩相流大量應(yīng)用于食品和藥物加工、石油化工和能源轉(zhuǎn)化等行業(yè).由于顆粒與顆粒之間及氣體與顆粒之間存在強(qiáng)烈的非線(xiàn)性耗散作用,密集顆粒流及氣粒流動(dòng)中往往會(huì)出現(xiàn)復(fù)雜的非均勻多尺度流動(dòng)現(xiàn)象[1-5].近年來(lái),國(guó)內(nèi)外諸多學(xué)者致力于采用數(shù)值模擬方法來(lái)研究這種復(fù)雜流動(dòng),目前針對(duì)顆粒相的數(shù)值模擬方法主要有離散單元法(DEM)和雙流體方法(TFM).DEM 相當(dāng)于顆粒相的直接數(shù)值模擬,計(jì)算準(zhǔn)確性高,但計(jì)算量大;而TFM 方法中顆粒相都被視為完全互穿的連續(xù)體,由單獨(dú)的質(zhì)量、動(dòng)量和能量守恒方程描述.這種顆粒擬流體的連續(xù)性表示,使得計(jì)算量不直接取決于顆粒數(shù),計(jì)算量相對(duì)較小.但需要提供額外的顆粒相本構(gòu)模型以及氣體與顆粒相間作用模型,而這些模型的可靠性將直接影響到TFM 計(jì)算結(jié)果的準(zhǔn)確性.

顆粒動(dòng)力學(xué)理論[6-8](KTGF)被廣泛用于推導(dǎo)顆粒相本構(gòu)模型,其中包含了描述顆粒與顆粒作用的本構(gòu)模型和描述顆粒與固體壁面作用的壁面邊界條件.對(duì)于后者,理論研究相對(duì)較少,當(dāng)前應(yīng)用最廣泛的是Johnson 和Jackson 邊界條件[9],它包含彈性恢復(fù)系數(shù)e和光滑因子φ兩個(gè)輸入?yún)?shù),分別用以描述顆粒與壁面碰撞時(shí)法向和切向上的速度變化.光滑因子φ的直接實(shí)驗(yàn)測(cè)量是不可行的,通常都是通過(guò)調(diào)整數(shù)值來(lái)擬合實(shí)驗(yàn)數(shù)據(jù)而獲得,這非常耗時(shí)且適用范圍窄[10-11].為此,一些學(xué)者嘗試使用顆粒物性和流動(dòng)參數(shù)來(lái)建立光滑因子φ的函數(shù).Li 和Benyahia[12]提出了一個(gè)適用于低摩擦系數(shù)下的光滑因子函數(shù)關(guān)聯(lián)式,其中光滑因子被擬合為壁面上顆粒法向恢復(fù)系數(shù)、摩擦系數(shù)和無(wú)因次滑移速度的函數(shù).Jenkins[13]引入3 個(gè)可測(cè)物性參數(shù),庫(kù)倫摩擦系數(shù)μ、切向彈性恢復(fù)系數(shù) β 和法向彈性恢復(fù)系數(shù)e,來(lái)共同描述微觀(guān)顆粒與壁面的碰撞行為,并通過(guò)假設(shè)近壁顆粒速度滿(mǎn)足正態(tài)分布,建立了小摩擦和大摩擦兩種極限情況下的壁面邊界條件.Schneiderbauer等[14]進(jìn)一步推導(dǎo)了從小摩擦到大摩擦全域統(tǒng)一的邊界條件表達(dá)式.Zhao 等[15]和Yang 等[16]考慮顆粒旋轉(zhuǎn)運(yùn)動(dòng),建立了適用于粗糙顆粒的壁面邊界條件,并發(fā)現(xiàn)考慮顆粒旋轉(zhuǎn)效應(yīng)后,模型預(yù)測(cè)與Louge[17]模擬結(jié)果的一致性要優(yōu)于之前的模型.由此可見(jiàn),可靠的壁面邊界條件需要關(guān)聯(lián)各種顆粒物性參數(shù)和更詳盡的顆粒流動(dòng)狀態(tài)參數(shù),但由于壁面附近顆粒屬性真實(shí)分布函數(shù)的缺乏以及三維碰撞過(guò)程數(shù)學(xué)積分求解的復(fù)雜,想要通過(guò)理論精確推導(dǎo)出統(tǒng)一的壁面邊界條件較為困難.因而在前期研究中,通常假設(shè)顆粒速度滿(mǎn)足正態(tài)分布[12-16],并且忽略顆粒旋轉(zhuǎn)效應(yīng)[12-14];而考慮顆粒旋轉(zhuǎn)效應(yīng)時(shí)又會(huì)增加額外的旋轉(zhuǎn)變量,在TFM 中需要求解附加的輸運(yùn)方程進(jìn)行封閉[15-16].為此,本文嘗試在常規(guī)TFM 中對(duì)近壁顆粒旋轉(zhuǎn)變量進(jìn)行局部代數(shù)型封閉,從而避免在整個(gè)流場(chǎng)中求解旋轉(zhuǎn)變量的輸運(yùn)方程.

近年來(lái),隨著計(jì)算機(jī)技術(shù)不斷提高,使得機(jī)器學(xué)習(xí)等人工智能技術(shù)也快速發(fā)展,其中作為機(jī)器學(xué)習(xí)研究熱門(mén)之一的人工神經(jīng)網(wǎng)絡(luò)模型(ANN)已經(jīng)發(fā)展得較為完善.人工神經(jīng)網(wǎng)絡(luò)是對(duì)人腦組織和運(yùn)行機(jī)制的某種抽象、簡(jiǎn)化和模擬,根據(jù)已知的一系列訓(xùn)練集,利用黑箱式學(xué)習(xí)方法的一種高效的數(shù)據(jù)處理和預(yù)測(cè)方法,具有效率高、聯(lián)想記憶、預(yù)測(cè)效果好等優(yōu)點(diǎn),因此,人工神經(jīng)網(wǎng)絡(luò)模型逐漸應(yīng)用于不同學(xué)科的研究之中.如胡洲等[18]采用神經(jīng)網(wǎng)絡(luò)建立非球形散體顆粒的休止角模型,發(fā)現(xiàn)休止角隨顆粒形狀變量、摩擦因數(shù)的增加都呈現(xiàn)增大的趨勢(shì),與現(xiàn)有研究結(jié)果一致;閆盛楠等[19]采用人工神經(jīng)網(wǎng)絡(luò)模型,對(duì)非球形顆粒氣固曳力系數(shù)進(jìn)行了預(yù)測(cè)及分析,并將模擬結(jié)果同文獻(xiàn)中的實(shí)驗(yàn)數(shù)據(jù)進(jìn)行對(duì)比分析,結(jié)果表明,人工神經(jīng)網(wǎng)絡(luò)可用于非球形顆粒氣固曳力系數(shù)的預(yù)測(cè)研究.

因此,本文擬采用DEM 方法對(duì)顆粒流進(jìn)行直接模擬,獲得壁面附近顆粒微觀(guān)運(yùn)動(dòng)數(shù)據(jù)和壁面與顆粒直接作用力;在此基礎(chǔ)上統(tǒng)計(jì)分析顆粒的宏觀(guān)運(yùn)動(dòng)特征,為顆粒動(dòng)力學(xué)理論推導(dǎo)的基本假設(shè)提供參考,并基于人工神經(jīng)網(wǎng)絡(luò)模型挖掘顆粒旋轉(zhuǎn)運(yùn)動(dòng)變量與顆粒物性參數(shù)和平動(dòng)運(yùn)動(dòng)變量之間的隱含關(guān)系.旨在為常規(guī)TFM 方法建立更加可靠的壁面邊界條件尋求可行的方法和基礎(chǔ)數(shù)據(jù).

1 數(shù)學(xué)模型

在DEM 模擬中,每個(gè)顆粒的運(yùn)動(dòng)受牛頓第二定律支配,質(zhì)量為mi的顆粒的平動(dòng)和旋轉(zhuǎn)運(yùn)動(dòng)可由下式描述[20]

切向力的大小受到庫(kù)侖摩擦定律的限制.如果式(3)和式(4)計(jì)算結(jié)果滿(mǎn)足以下條件

式中,μ 為摩擦系數(shù),則碰撞接觸點(diǎn)上會(huì)發(fā)生滑動(dòng),此時(shí)切向接觸力按下式計(jì)算

2 模擬結(jié)果與討論

2.1 模擬工況

在研究顆粒與壁面作用關(guān)系時(shí),可以選擇相對(duì)簡(jiǎn)單的模擬對(duì)象,如滾筒、斜槽[23]或者庫(kù)埃特(Couette)流[17].本文以石油化工中常見(jiàn)的滾筒為模擬對(duì)象.事實(shí)上,已經(jīng)有許多學(xué)者對(duì)滾筒進(jìn)行各種實(shí)驗(yàn)和模擬研究,如顧叢匯等[24]通過(guò)實(shí)驗(yàn)和數(shù)值模擬研究絲狀散體顆粒在滾筒內(nèi)的停留時(shí)間;胡陳樞等[25]采用DEM 方法對(duì)滾筒內(nèi)二元顆粒在不同轉(zhuǎn)速下的運(yùn)動(dòng)進(jìn)行模擬;張立棟等[26]采用DEM 方法研究滾筒內(nèi)構(gòu)件對(duì)二元顆粒體系運(yùn)動(dòng)混合的影響,并分析其增混機(jī)理.但多數(shù)研究都是關(guān)注滾筒內(nèi)部顆粒宏觀(guān)運(yùn)動(dòng)特征,然而顆粒和壁面作用通常至關(guān)重要[27-29],對(duì)滾筒而言,顆粒和壁面間的相互作用驅(qū)動(dòng)著整個(gè)顆粒系統(tǒng)運(yùn)動(dòng).

本文模擬的滾筒為Parker 等[30]的實(shí)驗(yàn)裝置.為了最大限度的接近實(shí)驗(yàn),模擬采用與實(shí)驗(yàn)一致的全尺寸,顆粒則簡(jiǎn)化為均一直徑.滾筒壁和顆粒材料均為有機(jī)玻璃,相關(guān)結(jié)構(gòu)和物性參數(shù)如表1 所示.模擬運(yùn)行18 s,時(shí)間步長(zhǎng)為5.6 × 10-6s,為避免初始效應(yīng)的影響,用于統(tǒng)計(jì)的數(shù)據(jù)樣本取自8～ 18 s.

表1 模擬參數(shù)Table 1 Simulation parameters

圖1 為瞬時(shí)顆粒分布,其中顏色表示速度大小.圖中顆粒速度在軸向具有較好的相似性,但在兩端面附近表層顆粒速度略有增加.Zhang 等[31]研究表明端壁摩擦對(duì)顆粒軸向擴(kuò)散具有增強(qiáng)效應(yīng),但其研究的滾筒長(zhǎng)徑比為1,而本模擬的滾筒長(zhǎng)徑比為6.5,因此端壁摩擦導(dǎo)致的軸向非均勻性較弱,模擬結(jié)果在軸向上具有較好的對(duì)稱(chēng)性.

圖1 滾筒內(nèi)的顆粒分布Fig.1 Particle distribution in the drum

圖2 為各轉(zhuǎn)速下滾筒最低點(diǎn)偏右θ=30°位置上顆粒切向速度模擬和實(shí)驗(yàn)結(jié)果對(duì)比圖.從圖2(a)中可以看出,摩擦系數(shù)取為0.7 時(shí)得到的切向速度曲線(xiàn)與實(shí)驗(yàn)結(jié)果一致;圖2(b)中針對(duì)不同轉(zhuǎn)速下也獲得較好的預(yù)測(cè)結(jié)果,特別是在近壁區(qū)域,模擬預(yù)測(cè)與實(shí)驗(yàn)結(jié)果吻合較好.通過(guò)DEM 模擬獲得接近物理實(shí)際的結(jié)果,將保障后續(xù)數(shù)據(jù)統(tǒng)計(jì)分析結(jié)果的可靠性.

圖2 切向速度與實(shí)驗(yàn)數(shù)據(jù)的比較Fig.2 Comparison of the tangential velocity with experiment data

在統(tǒng)計(jì)近壁處顆粒信息時(shí),為保證結(jié)果的獨(dú)立性,選取貼壁網(wǎng)格應(yīng)該盡量要小,同時(shí)又要保證一定的顆粒數(shù)以減小統(tǒng)計(jì)誤差.考慮到顆粒在滾筒軸向的對(duì)稱(chēng)性,貼壁網(wǎng)格在軸向取滾筒全長(zhǎng),而在橫截面上則如圖3 陰影區(qū)域所示,由周向尺寸δθ和徑向尺寸Δ共同確定.其中徑向尺寸Δ決定了統(tǒng)計(jì)顆?？拷诿娴某潭?對(duì)近壁顆粒的統(tǒng)計(jì)特性更為重要.圖3顯示了貼壁網(wǎng)格的徑向尺寸Δ與平均顆粒變量(顆粒速度vt和 ωz以及顆粒溫度T和 Θ)的關(guān)系.其中顆粒溫度是顆粒速度脈動(dòng)程度的度量,顆粒平動(dòng)溫度和旋轉(zhuǎn)溫度的定義和統(tǒng)計(jì)公式如下所示

式中,和分別是顆粒平動(dòng)和旋轉(zhuǎn)速度的平均值,運(yùn)算符 〈〉 表示系綜平均,N為貼壁網(wǎng)格內(nèi)的顆粒樣本總數(shù).從圖3 中可以看出,當(dāng)貼壁網(wǎng)格周向尺寸δθ=10°,徑向尺寸Δ達(dá)到1.8 mm 時(shí)平均顆粒速度和溫度都趨于穩(wěn)定,因此選取貼壁網(wǎng)格大小δθ=10°和Δ=1.8 mm 進(jìn)行后續(xù)分析研究.

圖3 網(wǎng)格大小與平均顆粒變量的關(guān)系Fig.3 Relationship between mesh size and averaged particle variables

2.2 顆粒微觀(guān)運(yùn)動(dòng)

2.2.1 平動(dòng)速度分布

在歐拉方法中,通常由顆粒動(dòng)力學(xué)理論提供描述顆粒運(yùn)動(dòng)特性的本構(gòu)關(guān)系.顆粒動(dòng)力學(xué)借鑒分子運(yùn)動(dòng)論,假設(shè)顆粒運(yùn)動(dòng)速度近似滿(mǎn)足正態(tài)分布,顆粒壁面邊界條件也是基于這樣的假設(shè)[13].對(duì)顆粒運(yùn)動(dòng)信息進(jìn)行統(tǒng)計(jì)分析,發(fā)現(xiàn)近壁處顆粒的平動(dòng)速度的確較好地符合正態(tài)分布,如圖4 所示.圖中統(tǒng)計(jì)參數(shù)取自滾筒最低點(diǎn)偏右θ=30°貼壁網(wǎng)格(網(wǎng)格見(jiàn)圖3)內(nèi)的顆粒.由于受到壁面剪切作用,切向的顆粒平動(dòng)速度均值和標(biāo)準(zhǔn)差SD最大;而在壁面阻礙和軸向?qū)ΨQ(chēng)性作用下,徑向和軸向的顆粒平動(dòng)速度均值趨近零,標(biāo)準(zhǔn)差SD也相對(duì)較小.顆粒速度分布的標(biāo)準(zhǔn)差反映了顆粒速度脈動(dòng)的強(qiáng)弱,顯然,在壁面作用下顆粒速度脈動(dòng)呈現(xiàn)出各向異性.

圖4 近壁顆粒平動(dòng)速度分布Fig.4 Translational velocity distribution of near-wall particles

在研究壁面邊界條件時(shí),摩擦系數(shù)對(duì)壁面顆粒運(yùn)動(dòng)具有顯著影響[31].表2 給出了不同摩擦系數(shù)下,壁面附近顆粒平動(dòng)速度分布的統(tǒng)計(jì)結(jié)果.其中決定系數(shù)R2可以表征正態(tài)分布函數(shù)擬合的好壞,取值范圍為[0,1],越接近1 表明數(shù)據(jù)擬合地越好.由表2 可知,4 個(gè)摩擦系數(shù)下的決定系數(shù)均大于0.92,表明顆粒平動(dòng)速度都較好地符合正態(tài)分布;此外,不同摩擦系數(shù)下,切向的標(biāo)準(zhǔn)差SD總是最大的,表明顆粒平動(dòng)速度脈動(dòng)存在各向異性[32].

表2 平動(dòng)速度分布的標(biāo)準(zhǔn)差和決定系數(shù)Table 2 Standard deviation and determination coefficient of translational velocity distribution

隨著摩擦系數(shù)的減小,各個(gè)方向的標(biāo)準(zhǔn)差SD有減小的趨勢(shì);當(dāng)摩擦系數(shù)降為0.3 時(shí),由于顆粒與壁面的摩擦作用減小,顆粒受壁面剪切激發(fā)的程度也相對(duì)減弱,顆粒速度脈動(dòng)也因此減弱.這在顆粒動(dòng)力學(xué)理論中表現(xiàn)為顆粒從宏觀(guān)運(yùn)動(dòng)中獲得的能量減小,從而導(dǎo)致顆粒溫度下降.此外,摩擦系數(shù)減小后,除軸向速度外,切向和徑向速度的決定系數(shù)也明顯減小,即正態(tài)分布假設(shè)的可靠性有所減弱;這表明近壁顆粒速度分布的主要影響因素仍然是壁面的徑向阻礙作用,其次才是壁面切向的摩擦剪切作用.因此,在顆粒動(dòng)力學(xué)理論中考慮顆粒與壁面作用時(shí),引入顆粒平動(dòng)速度正態(tài)分布假設(shè)基本上是合理的,但更精細(xì)的模型還應(yīng)進(jìn)一步考慮顆粒平動(dòng)速度脈動(dòng)的各向異性.

2.2.2 旋轉(zhuǎn)速度分布

滾筒最低點(diǎn)偏右θ=30°貼壁網(wǎng)格內(nèi)顆粒旋轉(zhuǎn)速度分布如圖5 所示,其中 ωx和 ωy兩個(gè)旋轉(zhuǎn)速度可以通過(guò)坐標(biāo)旋轉(zhuǎn)得到切向和徑向旋轉(zhuǎn)速度,因此它們的分布規(guī)律共同體現(xiàn)了切向和徑向旋轉(zhuǎn)速度分布特征.由圖5 可知,x和y方向上旋轉(zhuǎn)速度較好的滿(mǎn)足正態(tài)分布,且均值接近零、標(biāo)準(zhǔn)差SD也大小相當(dāng);但軸向 ωz均值和標(biāo)準(zhǔn)差都相對(duì)較大,并且離正態(tài)分布有較大的偏離.因此,在動(dòng)力學(xué)理論中直接將顆粒旋轉(zhuǎn)速度假設(shè)為正態(tài)分布,將產(chǎn)生一定的誤差.根據(jù)Schneiderbauer 等[14]研究,顆粒在壁面上的旋轉(zhuǎn)速度ωz主要由顆粒與壁面的平動(dòng)滑移速度Vslip=-vwall驅(qū)動(dòng),如果顆粒速度滿(mǎn)足正態(tài)分布,則 ωz≈μVslip/d.根據(jù)圖4 中的平動(dòng)速度均值及滾筒轉(zhuǎn)速可以計(jì)算出ωz≈30 rad/s,與模擬統(tǒng)計(jì)均值約25.7 rad/s 相比誤差明顯.

圖5 近壁顆粒旋轉(zhuǎn)速度分布Fig.5 Rotational velocity distribution of near-wall particles

圖6 為摩擦系數(shù)為0.3,0.5 和0.9 時(shí),顆粒軸向旋轉(zhuǎn)速度分布.結(jié)合圖5(c),可以發(fā)現(xiàn)隨著摩擦系數(shù)的提高,顆粒軸向旋轉(zhuǎn)速度 ωz分布逐漸偏離正態(tài)分布;在摩擦系數(shù)為0.9 時(shí),ωz分布明顯趨近于雙峰分布.在Parker 等[30]的實(shí)驗(yàn)分析中,對(duì)Δ=6 mm 的貼壁網(wǎng)格也統(tǒng)計(jì)出雙峰分布的結(jié)果,并且雙峰效應(yīng)隨著滾筒轉(zhuǎn)速的增加而增強(qiáng);其解釋為貼壁網(wǎng)格徑向尺寸為6 mm 至少包含兩層直徑為3 mm 的顆粒,不同層顆粒因受壁面剪切影響程度不同而具有各自不同的旋轉(zhuǎn)特性,從而導(dǎo)致在Δ=6 mm 的貼壁網(wǎng)格中出現(xiàn)雙峰分布.而本文貼壁網(wǎng)格Δ=1.8 mm,所統(tǒng)計(jì)顆粒都貼近壁面,出現(xiàn)雙峰分布的主要原因仍是剪切效應(yīng);此時(shí),貼壁網(wǎng)格內(nèi)的顆粒受到下層壁面和上層顆粒的兩面剪切,不同摩擦系數(shù)下兩面剪切的影響程度不同,當(dāng)摩擦系數(shù)增大時(shí),兩面剪切差異增加從而產(chǎn)生雙峰結(jié)果.

圖6 不同摩擦系數(shù)下的近壁顆粒軸向旋轉(zhuǎn)速度分布Fig.6 Axial rotational velocity distribution of near-wall particles under different friction coefficients

表3 為不同摩擦系數(shù)下的顆粒旋轉(zhuǎn)速度統(tǒng)計(jì)分析結(jié)果.從表3 中可以看出,摩擦系數(shù)為0.3 時(shí),旋轉(zhuǎn)速度分布的擬合決定系數(shù)接近1;但隨著摩擦系數(shù)增加,軸向旋轉(zhuǎn)速度分布的擬合決定系數(shù)明顯減小,摩擦系數(shù)在0.5,0.7 和0.9 時(shí),擬合決定系數(shù)小于0.9,表明其嚴(yán)重偏離正態(tài)分布,這與圖6 旋轉(zhuǎn)速度的直觀(guān)分布是一致的.同時(shí),不同摩擦系數(shù)下,軸向旋轉(zhuǎn)速度分布的標(biāo)準(zhǔn)差SD均明顯大于其他兩個(gè)方向,表現(xiàn)出較強(qiáng)的各向異性.總的來(lái)說(shuō),切向和徑向的旋轉(zhuǎn)速度分布仍然較好的滿(mǎn)足正態(tài)分布,但軸向的旋轉(zhuǎn)速度分布隨著摩擦系數(shù)的增加越偏離正態(tài)分布.因此,與顆粒平動(dòng)速度分布不同,在動(dòng)力學(xué)理論中將壁面附近顆粒旋轉(zhuǎn)速度假設(shè)為正態(tài)分布需謹(jǐn)慎.

表3 旋轉(zhuǎn)速度分布的標(biāo)準(zhǔn)差和決定系數(shù)Table 3 Standard deviation and determination coefficient of rotational velocity distribution

2.3 無(wú)因次旋轉(zhuǎn)溫度

在顆粒動(dòng)力學(xué)理論中,顆粒溫度將單個(gè)顆粒的微觀(guān)小尺度脈動(dòng)行為與大量顆粒所表現(xiàn)的宏觀(guān)流動(dòng)性質(zhì)相關(guān)聯(lián),出現(xiàn)在顆粒黏度、壓力和擴(kuò)散系數(shù)等參數(shù)的本構(gòu)關(guān)系中.實(shí)際顆粒碰撞過(guò)程中,平動(dòng)能量與旋轉(zhuǎn)能量會(huì)相互轉(zhuǎn)化,因此顆粒平動(dòng)溫度與旋轉(zhuǎn)溫度應(yīng)相互耦合.顆粒與壁面作用時(shí),旋轉(zhuǎn)溫度的影響是不可忽略的[14],但常規(guī)歐拉雙流體模型中并未求解顆粒旋轉(zhuǎn)溫度,為此需建立顆粒旋轉(zhuǎn)溫度的代數(shù)型封閉關(guān)聯(lián)式.Jenkins 和Zhang[33]曾將球形顆粒無(wú)因次旋轉(zhuǎn)溫度 λ =2.5Θ/T關(guān)聯(lián)為顆粒物性參數(shù)的函數(shù);Zhao 等[15]曾將λ關(guān)聯(lián)為顆粒物性參數(shù)和無(wú)因次滑移速度的函數(shù).但這些封閉關(guān)聯(lián)式,主要是根據(jù)理論簡(jiǎn)化假設(shè)或簡(jiǎn)單數(shù)據(jù)擬合,可能忽略了某些起作用的因素.通過(guò)DEM 模擬,可以得到大量的近壁顆粒運(yùn)動(dòng)信息,傳統(tǒng)方法在處理這些數(shù)據(jù)時(shí)較為吃力;而數(shù)據(jù)挖掘方法不僅能夠處理大量數(shù)據(jù),同時(shí)能夠探索數(shù)據(jù)之間的隱藏規(guī)律,神經(jīng)網(wǎng)絡(luò)就是數(shù)據(jù)挖掘常用的一種方法.因此本文將DEM 模擬得到的顆粒運(yùn)動(dòng)信息包括顆粒速度和顆粒溫度等作為神經(jīng)網(wǎng)絡(luò)的輸入,采用人工神經(jīng)網(wǎng)絡(luò)進(jìn)行顆粒無(wú)因次旋轉(zhuǎn)溫度的預(yù)測(cè)學(xué)習(xí),嘗試建立顆粒無(wú)因次旋轉(zhuǎn)溫度的封閉關(guān)聯(lián)式.

圖7 為本文采用的BP 神經(jīng)網(wǎng)絡(luò)結(jié)構(gòu)示意圖.BP 神經(jīng)網(wǎng)絡(luò)是一種多層的前饋神經(jīng)網(wǎng)絡(luò),其主要的特點(diǎn)是:信號(hào)是前向傳播的,而誤差是反向傳播的,傳播過(guò)程主要分為兩個(gè)階段,第一階段是信號(hào)的前向傳播,從輸入層經(jīng)過(guò)隱含層,最后到達(dá)輸出層;第二階段是誤差的反向傳播,從輸出層到隱含層,最后到輸入層,依次調(diào)節(jié)隱含層到輸出層的權(quán)重和偏置,輸入層到隱含層的權(quán)重和偏置.每層的神經(jīng)元節(jié)點(diǎn)通過(guò)權(quán)值與偏置與下一層的節(jié)點(diǎn)相連接,輸入信號(hào)通過(guò)激活函數(shù)轉(zhuǎn)換成輸出結(jié)果,最后網(wǎng)絡(luò)輸出如下

圖7 本文使用的神經(jīng)網(wǎng)絡(luò)結(jié)構(gòu)Fig.7 Structure of neural network used in present work

式中,yp為輸出值,對(duì)于隱含層下標(biāo)p為神經(jīng)元序號(hào),xpq則為第p個(gè)神經(jīng)元的q個(gè)輸入,wpq和bp為對(duì)應(yīng)的權(quán)值和偏置,Q為該層輸入值總數(shù).f是對(duì)應(yīng)層的激活函數(shù),本文隱含層使用Sigmoid 函數(shù),輸出層使用Linear 函數(shù).

當(dāng)神經(jīng)網(wǎng)絡(luò)輸出值與目標(biāo)值y? 不等時(shí),存在輸出誤差E,定義如下

為了減小誤差E,BP 神經(jīng)網(wǎng)絡(luò)使用梯度下降反向傳播算法對(duì)輸入與隱含層間以及隱含層與輸出層間的權(quán)值和偏置進(jìn)行調(diào)整.

本研究輸入?yún)?shù)為平動(dòng)溫度T和滑移速度Vslip,預(yù)期輸出值為無(wú)因次旋轉(zhuǎn)溫度 λ .為獲取大范圍滑移速度樣本,選取角度θ為0°,10°,20°,30°,40°,50°和60°的貼壁網(wǎng)格.為加快收斂速度,采用最大-最小標(biāo)準(zhǔn)化方法處理訓(xùn)練樣本.訓(xùn)練時(shí),以隨機(jī)數(shù)對(duì)神經(jīng)網(wǎng)絡(luò)權(quán)值和偏置進(jìn)行賦值,然后根據(jù)誤差E反向傳播對(duì)權(quán)值與偏置進(jìn)行優(yōu)化調(diào)整.本文BP 神經(jīng)網(wǎng)絡(luò)模型的具體細(xì)節(jié)見(jiàn)表4 所示,訓(xùn)練完的代數(shù)模型數(shù)學(xué)表達(dá)式如下

表4 人工神經(jīng)網(wǎng)絡(luò)模型的細(xì)節(jié)Table 4 Details of the artificial neural network model

圖8 為模型預(yù)測(cè)結(jié)果與DEM 統(tǒng)計(jì)結(jié)果的對(duì)比,其中紅色圓點(diǎn)為BP 神經(jīng)網(wǎng)絡(luò)預(yù)測(cè)的無(wú)因次旋轉(zhuǎn)溫度λ,黑色的正方形為DEM 模擬結(jié)果,δ為統(tǒng)計(jì)的平均相對(duì)誤差.由圖可見(jiàn),相同摩擦系數(shù)下,無(wú)因次旋轉(zhuǎn)溫度隨無(wú)因次滑移速度增大呈現(xiàn)出增長(zhǎng)趨勢(shì);且相同無(wú)因次滑移速度下,隨著摩擦系數(shù)的增加,無(wú)因次旋轉(zhuǎn)溫度也表現(xiàn)為增長(zhǎng)趨勢(shì).此外,即使樣本數(shù)據(jù)點(diǎn)較為離散,BP 神經(jīng)網(wǎng)絡(luò)預(yù)測(cè)結(jié)果的平均相對(duì)誤差依然較小,表明采用神經(jīng)網(wǎng)絡(luò)模型建立顆粒旋轉(zhuǎn)溫度的代數(shù)型封閉關(guān)聯(lián)式是可行的.

圖8 無(wú)因次旋轉(zhuǎn)溫度的模型預(yù)測(cè)與DEM 統(tǒng)計(jì)結(jié)果對(duì)比Fig.8 Comparison of predicted dimensionless rotational temperature with DEM results

圖8 無(wú)因次旋轉(zhuǎn)溫度的模型預(yù)測(cè)與DEM 統(tǒng)計(jì)結(jié)果對(duì)比(續(xù))Fig.8 Comparison of predicted dimensionless rotational temperature with DEM results (continued)

2.4 壁面處顆粒應(yīng)力分析

歐拉雙流體模型中,顆粒壁面邊界條件主要是通過(guò)壁面上顆粒相法向應(yīng)力來(lái)關(guān)聯(lián)切向應(yīng)力.因此,通過(guò)DEM 模擬并統(tǒng)計(jì)分析壁面所受顆粒的切向和法向應(yīng)力,可以為雙流體模型中顆粒壁面邊界條件的構(gòu)建提供數(shù)據(jù)參考.

圖9 給出了不同摩擦系數(shù)下,壁面上切向和法向應(yīng)力比與顆粒無(wú)因次滑移速度之間的關(guān)系.由圖可知,當(dāng)摩擦系數(shù)為0.3 和0.5 時(shí),應(yīng)力比隨無(wú)因次滑移速度的增加而增大,并趨近于摩擦系數(shù)值;當(dāng)摩擦系數(shù)增大到0.7 和0.9 時(shí),無(wú)因次滑移速度的增加到一定程度后不再增加,應(yīng)力比出現(xiàn)非單值型分布,這和圖中Louge[17]針對(duì)摩擦系數(shù)0.75 和 ∞ 兩種狀態(tài)下給出的分布趨勢(shì)是一致的.本文統(tǒng)計(jì)結(jié)果和Louge 的數(shù)據(jù)分布規(guī)律相同,但數(shù)值上存在差異,并隨摩擦系數(shù)增加差異變大,其原因可能在于Louge的模擬中壁面上鑲嵌有半球體,增大了顆粒與壁面的切向作用.

圖9 不同摩擦系數(shù)下壁面切向與法向應(yīng)力比Fig.9 Wall stress ratio under different friction

Johnson 和Jackson 壁面邊界條件采用簡(jiǎn)單常數(shù)來(lái)描述應(yīng)力比與無(wú)因次滑移速度的關(guān)系,這與實(shí)際數(shù)據(jù)是不相符的,只能作為一定范圍內(nèi)的近似或平均.Zhao 等[34]對(duì)近壁顆粒平動(dòng)和旋轉(zhuǎn)速度采用各向同性的正態(tài)分布,推導(dǎo)了應(yīng)力比與滑移速度理論公式,對(duì)比Johnson 和Jackson 壁面邊界條件即可將光滑因子表示如下[15]

將神經(jīng)網(wǎng)絡(luò)模型預(yù)測(cè)的無(wú)因次旋轉(zhuǎn)溫度代入上述光滑因子,再結(jié)合Johnson 和Jackson 壁面邊界條件,可以計(jì)算出壁面上應(yīng)力比如圖9 中實(shí)線(xiàn)所示.當(dāng)摩擦系數(shù)為0.3 時(shí),由表2 和表3 可知顆粒速度分布的確較好的符合正態(tài)分布,因此公式預(yù)測(cè)結(jié)果與DEM 統(tǒng)計(jì)結(jié)果相吻合;但當(dāng)摩擦系數(shù)提高后,由于顆粒旋轉(zhuǎn)速度偏離正態(tài)分布,公式預(yù)測(cè)和DEM 統(tǒng)計(jì)結(jié)果偏差增加.

3 結(jié)論

采用DEM 方法數(shù)值模擬了滾筒內(nèi)顆粒流動(dòng)狀況.通過(guò)統(tǒng)計(jì)和ANN 算法研究了壁面附近微觀(guān)顆粒的運(yùn)動(dòng)特征,得出了以下結(jié)論.

(1) 顆粒平動(dòng)速度都基本符合正態(tài)分布,小摩擦系數(shù)下顆粒旋轉(zhuǎn)速度也滿(mǎn)足正態(tài)分布,但隨著摩擦系數(shù)的提高,壁面剪切產(chǎn)生的軸向旋轉(zhuǎn)速度分布偏離正態(tài)分布;同時(shí)顆粒平動(dòng)和旋轉(zhuǎn)速度脈動(dòng)呈現(xiàn)較強(qiáng)的各向異性.

(2) 無(wú)因次旋轉(zhuǎn)溫度和無(wú)因次滑移速度存在一定的正相關(guān),并受摩擦系數(shù)的影響.將滑移速度和顆粒平動(dòng)溫度作為輸入,采用人工神經(jīng)網(wǎng)絡(luò)構(gòu)建顆粒旋轉(zhuǎn)溫度模型,可以為常規(guī)歐拉雙流體模型框架下考慮顆粒旋轉(zhuǎn)溫度提供代數(shù)型封閉.

(3) 壁面與顆粒作用的切向與法向應(yīng)力比,受無(wú)因次滑移速度和摩擦系數(shù)的影響,摩擦系數(shù)增大后,其理論關(guān)系的建立需要考慮更真實(shí)的顆粒速度分布,如非正態(tài)分布、各向異性.