基于離散單元法和人工神經網絡的近壁顆粒動力學特征研究1)

段總樣 趙云華 徐 璋

(浙江工業大學機械工程學院,杭州 310014)

引言

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

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

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

因此,本文擬采用DEM 方法對顆粒流進行直接模擬,獲得壁面附近顆粒微觀運動數據和壁面與顆粒直接作用力;在此基礎上統計分析顆粒的宏觀運動特征,為顆粒動力學理論推導的基本假設提供參考,并基于人工神經網絡模型挖掘顆粒旋轉運動變量與顆粒物性參數和平動運動變量之間的隱含關系.旨在為常規TFM 方法建立更加可靠的壁面邊界條件尋求可行的方法和基礎數據.

1 數學模型

在DEM 模擬中,每個顆粒的運動受牛頓第二定律支配,質量為mi的顆粒的平動和旋轉運動可由下式描述[20]

切向力的大小受到庫侖摩擦定律的限制.如果式(3)和式(4)計算結果滿足以下條件

式中,μ 為摩擦系數,則碰撞接觸點上會發生滑動,此時切向接觸力按下式計算

2 模擬結果與討論

2.1 模擬工況

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

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

表1 模擬參數Table 1 Simulation parameters

圖1 為瞬時顆粒分布,其中顏色表示速度大小.圖中顆粒速度在軸向具有較好的相似性,但在兩端面附近表層顆粒速度略有增加.Zhang 等[31]研究表明端壁摩擦對顆粒軸向擴散具有增強效應,但其研究的滾筒長徑比為1,而本模擬的滾筒長徑比為6.5,因此端壁摩擦導致的軸向非均勻性較弱,模擬結果在軸向上具有較好的對稱性.

圖1 滾筒內的顆粒分布Fig.1 Particle distribution in the drum

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

圖2 切向速度與實驗數據的比較Fig.2 Comparison of the tangential velocity with experiment data

在統計近壁處顆粒信息時,為保證結果的獨立性,選取貼壁網格應該盡量要小,同時又要保證一定的顆粒數以減小統計誤差.考慮到顆粒在滾筒軸向的對稱性,貼壁網格在軸向取滾筒全長,而在橫截面上則如圖3 陰影區域所示,由周向尺寸δθ和徑向尺寸Δ共同確定.其中徑向尺寸Δ決定了統計顆粒靠近壁面的程度,對近壁顆粒的統計特性更為重要.圖3顯示了貼壁網格的徑向尺寸Δ與平均顆粒變量(顆粒速度vt和 ωz以及顆粒溫度T和 Θ)的關系.其中顆粒溫度是顆粒速度脈動程度的度量,顆粒平動溫度和旋轉溫度的定義和統計公式如下所示

式中,和分別是顆粒平動和旋轉速度的平均值,運算符 〈〉 表示系綜平均,N為貼壁網格內的顆粒樣本總數.從圖3 中可以看出,當貼壁網格周向尺寸δθ=10°,徑向尺寸Δ達到1.8 mm 時平均顆粒速度和溫度都趨于穩定,因此選取貼壁網格大小δθ=10°和Δ=1.8 mm 進行后續分析研究.

圖3 網格大小與平均顆粒變量的關系Fig.3 Relationship between mesh size and averaged particle variables

2.2 顆粒微觀運動

2.2.1 平動速度分布

在歐拉方法中,通常由顆粒動力學理論提供描述顆粒運動特性的本構關系.顆粒動力學借鑒分子運動論,假設顆粒運動速度近似滿足正態分布,顆粒壁面邊界條件也是基于這樣的假設[13].對顆粒運動信息進行統計分析,發現近壁處顆粒的平動速度的確較好地符合正態分布,如圖4 所示.圖中統計參數取自滾筒最低點偏右θ=30°貼壁網格(網格見圖3)內的顆粒.由于受到壁面剪切作用,切向的顆粒平動速度均值和標準差SD最大;而在壁面阻礙和軸向對稱性作用下,徑向和軸向的顆粒平動速度均值趨近零,標準差SD也相對較小.顆粒速度分布的標準差反映了顆粒速度脈動的強弱,顯然,在壁面作用下顆粒速度脈動呈現出各向異性.

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

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

表2 平動速度分布的標準差和決定系數Table 2 Standard deviation and determination coefficient of translational velocity distribution

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

2.2.2 旋轉速度分布

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

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

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

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

表3 為不同摩擦系數下的顆粒旋轉速度統計分析結果.從表3 中可以看出,摩擦系數為0.3 時,旋轉速度分布的擬合決定系數接近1;但隨著摩擦系數增加,軸向旋轉速度分布的擬合決定系數明顯減小,摩擦系數在0.5,0.7 和0.9 時,擬合決定系數小于0.9,表明其嚴重偏離正態分布,這與圖6 旋轉速度的直觀分布是一致的.同時,不同摩擦系數下,軸向旋轉速度分布的標準差SD均明顯大于其他兩個方向,表現出較強的各向異性.總的來說,切向和徑向的旋轉速度分布仍然較好的滿足正態分布,但軸向的旋轉速度分布隨著摩擦系數的增加越偏離正態分布.因此,與顆粒平動速度分布不同,在動力學理論中將壁面附近顆粒旋轉速度假設為正態分布需謹慎.

表3 旋轉速度分布的標準差和決定系數Table 3 Standard deviation and determination coefficient of rotational velocity distribution

2.3 無因次旋轉溫度

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

圖7 為本文采用的BP 神經網絡結構示意圖.BP 神經網絡是一種多層的前饋神經網絡,其主要的特點是:信號是前向傳播的,而誤差是反向傳播的,傳播過程主要分為兩個階段,第一階段是信號的前向傳播,從輸入層經過隱含層,最后到達輸出層;第二階段是誤差的反向傳播,從輸出層到隱含層,最后到輸入層,依次調節隱含層到輸出層的權重和偏置,輸入層到隱含層的權重和偏置.每層的神經元節點通過權值與偏置與下一層的節點相連接,輸入信號通過激活函數轉換成輸出結果,最后網絡輸出如下

圖7 本文使用的神經網絡結構Fig.7 Structure of neural network used in present work

式中,yp為輸出值,對于隱含層下標p為神經元序號,xpq則為第p個神經元的q個輸入,wpq和bp為對應的權值和偏置,Q為該層輸入值總數.f是對應層的激活函數,本文隱含層使用Sigmoid 函數,輸出層使用Linear 函數.

當神經網絡輸出值與目標值y? 不等時,存在輸出誤差E,定義如下

為了減小誤差E,BP 神經網絡使用梯度下降反向傳播算法對輸入與隱含層間以及隱含層與輸出層間的權值和偏置進行調整.

本研究輸入參數為平動溫度T和滑移速度Vslip,預期輸出值為無因次旋轉溫度 λ .為獲取大范圍滑移速度樣本,選取角度θ為0°,10°,20°,30°,40°,50°和60°的貼壁網格.為加快收斂速度,采用最大-最小標準化方法處理訓練樣本.訓練時,以隨機數對神經網絡權值和偏置進行賦值,然后根據誤差E反向傳播對權值與偏置進行優化調整.本文BP 神經網絡模型的具體細節見表4 所示,訓練完的代數模型數學表達式如下

表4 人工神經網絡模型的細節Table 4 Details of the artificial neural network model

圖8 為模型預測結果與DEM 統計結果的對比,其中紅色圓點為BP 神經網絡預測的無因次旋轉溫度λ,黑色的正方形為DEM 模擬結果,δ為統計的平均相對誤差.由圖可見,相同摩擦系數下,無因次旋轉溫度隨無因次滑移速度增大呈現出增長趨勢;且相同無因次滑移速度下,隨著摩擦系數的增加,無因次旋轉溫度也表現為增長趨勢.此外,即使樣本數據點較為離散,BP 神經網絡預測結果的平均相對誤差依然較小,表明采用神經網絡模型建立顆粒旋轉溫度的代數型封閉關聯式是可行的.

圖8 無因次旋轉溫度的模型預測與DEM 統計結果對比Fig.8 Comparison of predicted dimensionless rotational temperature with DEM results

圖8 無因次旋轉溫度的模型預測與DEM 統計結果對比(續)Fig.8 Comparison of predicted dimensionless rotational temperature with DEM results (continued)

2.4 壁面處顆粒應力分析

歐拉雙流體模型中,顆粒壁面邊界條件主要是通過壁面上顆粒相法向應力來關聯切向應力.因此,通過DEM 模擬并統計分析壁面所受顆粒的切向和法向應力,可以為雙流體模型中顆粒壁面邊界條件的構建提供數據參考.

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

圖9 不同摩擦系數下壁面切向與法向應力比Fig.9 Wall stress ratio under different friction

Johnson 和Jackson 壁面邊界條件采用簡單常數來描述應力比與無因次滑移速度的關系,這與實際數據是不相符的,只能作為一定范圍內的近似或平均.Zhao 等[34]對近壁顆粒平動和旋轉速度采用各向同性的正態分布,推導了應力比與滑移速度理論公式,對比Johnson 和Jackson 壁面邊界條件即可將光滑因子表示如下[15]

將神經網絡模型預測的無因次旋轉溫度代入上述光滑因子,再結合Johnson 和Jackson 壁面邊界條件,可以計算出壁面上應力比如圖9 中實線所示.當摩擦系數為0.3 時,由表2 和表3 可知顆粒速度分布的確較好的符合正態分布,因此公式預測結果與DEM 統計結果相吻合;但當摩擦系數提高后,由于顆粒旋轉速度偏離正態分布,公式預測和DEM 統計結果偏差增加.

3 結論

采用DEM 方法數值模擬了滾筒內顆粒流動狀況.通過統計和ANN 算法研究了壁面附近微觀顆粒的運動特征,得出了以下結論.

(1) 顆粒平動速度都基本符合正態分布,小摩擦系數下顆粒旋轉速度也滿足正態分布,但隨著摩擦系數的提高,壁面剪切產生的軸向旋轉速度分布偏離正態分布;同時顆粒平動和旋轉速度脈動呈現較強的各向異性.

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

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