馮 永 原子然
(河南工業(yè)大學土木建筑學院,鄭州 450001)
自CUNDALL提出離散元法以來,PFC,DEM等基于離散元法的數(shù)值模擬軟件廣泛運用于筒倉領域的研究,并取得了一定的成果[1-3]。筒倉的相關研究中也經(jīng)常采用PFC中的ball單元來模擬筒倉中的糧食顆粒[4-6],但是ball單元并不適用于所有實際工程狀況,常常存在宏觀力學結(jié)果與試驗結(jié)果誤差較大、顆粒細觀力學參數(shù)不全面的情況[7]。產(chǎn)生誤差的一個重要原因就是自然條件下的糧食顆粒大多不以規(guī)則的形狀存在,而傳統(tǒng)的單一ball單元模擬不能反映小麥、稻谷等不規(guī)則糧食儲料的形態(tài)差異對卸糧過程中的影響,因此難以準確反映筒倉卸糧過程中的宏細觀力學機理,針對這個問題,國內(nèi)外學者提出clump模型模擬玉米、小麥、煤等顆粒[8-10],但是對于顆粒形狀及含塵率在筒倉卸糧方面對宏細觀參數(shù)模擬精確度的影響,至今鮮見定量的研究,對于筒倉卸糧實驗中含塵率的影響考慮不足[11]。
基于以上分析,本研究提出了一種由clump單元和小型ball單元組成的顆粒簇單元模型,該模型用clump單元模擬糧食顆粒,并添加小型ball單元模擬粉塵顆粒。在推導顆粒之間接觸、顆粒-倉壁接觸本構(gòu)關系的基礎上[12,13],以室內(nèi)卸糧物理模型試驗為基礎,采用顆粒簇單元模型進行數(shù)值模擬試驗,并和傳統(tǒng)采用單一ball單元的模擬結(jié)果進行對比分析。
本研究旨在針對小麥、稻谷等不規(guī)則糧食顆粒,建立一種顆粒簇單元模型,更加量化模擬糧食實際形態(tài)、含塵率、孔隙率等對卸糧宏細觀力學參數(shù)影響,客觀準確地揭示卸糧過程中的宏細觀力學機理,并對類似不規(guī)則顆粒的數(shù)值模擬也提供參考。
圖1 改進顆粒簇單元示意圖
本研究提出的顆粒簇單元由clump單元和傳統(tǒng)小型ball單元組成, 該模型用clump單元模擬不規(guī)則糧食顆粒,添加小型ball單元模擬粉塵顆粒。
clump是以三個剛性小球相互疊合成一個整體來模擬除剛性圓球外的其他形狀。針對本文研究不規(guī)則糧食顆粒(小麥、稻谷)的形態(tài),擬合橢圓形狀如圖1所示。
2 接觸本構(gòu)模型
結(jié)合現(xiàn)有研究現(xiàn)狀,PFC顆?;A模型多為線性接觸模型[14,15],根據(jù)數(shù)值模擬需要,建立相應的力-位移更新公式。
將clump單元簡化成一個規(guī)則的橢圓形,根據(jù)力-位移的定律,得到兩顆粒接觸時的單位矢量模型和顆粒接觸變形的情況。認為接觸點處的法向切向向量為:
(1)
(2)
兩個顆粒相互接觸如圖2所示。
圖2 顆粒單元之間接觸
(3)
(4)
則在接觸點處兩顆粒的相對速度為:
(5)
(6)
相對速度在法向和切向的分量如式(7)和式(8)所示。
(7)
(8)
(9)
(10)
(11)
相對速度在法向和切向的分量為:
(12)
(13)
圖3 clump單元和墻體之間接觸
對于顆粒和墻體,法向切向的速度分量都造成了相應方向的位移,根據(jù)上述方式所計算出的速度分量,可以得到一個時間段變化的位移分量為Δr,根據(jù)力-位移定律[16,17]可以得到單位時間內(nèi)力的變化量。
(14)
(15)
由此易得顆粒的運動方程:
(16)
(17)
式中:I為顆粒的慣性矩,∑F為顆粒體所受的合力。
在運算的每個微小時段,細觀參數(shù)都在發(fā)生改變。以顆粒A為例,力在一個微小時段的更新為:
(18)
(19)
速度在一個微小時段的更新為:
(20)
(21)
顆粒位置與旋轉(zhuǎn)角度在一個微小時段的更新為
(22)
(23)
按照實際常用筒倉合理尺寸20 ∶1的縮尺比例制作模型,試驗筒倉高1 000 mm,半徑250 mm,倉壁厚約5 mm,漏斗半頂角60°。在倉壁布置監(jiān)測點,每個監(jiān)測點都布置一個壓力感應器來檢測筒倉卸糧成拱過程中的側(cè)壓力,筒倉模型的具體尺寸和監(jiān)測點布置由圖4所示 。
圖4 筒倉模型尺寸和監(jiān)測點布置
小麥試樣選用河南產(chǎn)小麥,小麥長軸粒徑分布區(qū)間為5~6 mm,短軸粒徑分布區(qū)間為3~4 mm,試樣顆粒尺寸分布比較均勻。
所選用的壓力感應裝置為電阻式應變片,型號為B×120-2AA,采集儀器為DH3816N應變測試分析系統(tǒng)。
2018年9月15日進行裝糧,將試驗模型倉內(nèi)裝滿小麥,待裝料完成并靜置兩日后于9月18日撤掉料斗底板進行卸糧。倉壁壓力感應器將卸糧過程中倉壁的微小位移通過無氧銅導線輸入數(shù)據(jù)采集系統(tǒng)。
4 數(shù)值模擬
根據(jù)室內(nèi)試驗的筒倉尺寸,建立一個半徑250 mm,高1 000 mm的筒倉模型。用wall單元模擬墻體,用長軸直徑6 mm,短軸直徑4 mm的clump單元模擬橢圓形小麥顆粒,以10 ∶1的比例添加用直徑0.5 mm的ball單元模擬的粉塵顆粒。(原模型儲料為直徑5 mm的單一ball單元),單一ball單元模型中所采用的ball單元體積為19.63 mm3,改進顆粒簇模型中的clump單元體積為18.84 mm3,相比降低了4.02%。
在筒倉模型左右兩側(cè)遵循試驗中監(jiān)測點的布置分段建立墻體,監(jiān)控并記錄各wall單元的X軸方向所受到的壓力,與試驗中監(jiān)測點監(jiān)測記錄側(cè)壓力的作用相同。
根據(jù)實驗室測定以及糧食力學特性[18]選取參數(shù),并通過參數(shù)標定[19,20]方法,使靜態(tài)儲糧狀態(tài)下模擬值與實測值相符合,進一步調(diào)整選取參數(shù),最終確定選定的物理力學參數(shù)如表 1所示。wall,ball和clump分別表示墻體單元,圓顆粒單元和不規(guī)則顆粒單元,kn、ks則分別為它們的法向和切向剛度。
表1 模擬中顆粒和墻體力學參數(shù)
墻體和顆粒的線性接觸模型所選取的參數(shù)如表 2所示。
表2 模擬中主要物性參數(shù)
根據(jù)計算機運算能力和試驗條件進行參數(shù)取值后,使單一ball單元或改進顆粒簇單元(改進顆粒簇模型為clump單元和小型ball單元,數(shù)量比例大約為10 ∶1)以中心裝料的方法循環(huán)落入筒倉,設定線性接觸模型后,顆粒之間將按照第二章的接觸模型和更新方式進行運算。
由于如果顆粒內(nèi)摩擦力較大,則會發(fā)生裝料緩慢,計算復雜的狀況,所以在本研究中采用一種新的裝料方法,即在顆粒單元下落過程中不設置顆粒間的摩擦力,待分層裝料完成后,用ball property(clump property)命令賦予顆粒體摩擦系數(shù),再迭代一定時步達到穩(wěn)定。單一ball單元模型裝載完成后共含有23 000 個ball單元,堆料高度為0.994 m。改進顆粒簇模型裝載完成后含有24 676 個clump單元,2 382個小型ball單元堆料高度為1.013 m。兩模型規(guī)模相近。
選取左側(cè)點位將模擬所得靜態(tài)儲糧側(cè)壓力結(jié)果與室內(nèi)物理模型試驗的靜態(tài)儲糧側(cè)壓力結(jié)果對比(下圖5),結(jié)果表示數(shù)值相差不大,而整體趨勢比較契合,可以驗證此模擬試驗的真實性。
圖5 靜態(tài)裝載下倉壁側(cè)壓力值
考慮到上部監(jiān)測點側(cè)壓力波動幅度較小且在卸糧開始一段時間后就進入零壓力區(qū)而導致不利于觀測的問題,選定 1、3、5三個監(jiān)測點進行卸糧動態(tài)側(cè)壓力的數(shù)據(jù)分析。
改進后模型和原單一ball模型在自由卸糧過程中的倉壁動態(tài)側(cè)壓力如下圖6、圖7所示。
由圖可以發(fā)現(xiàn),在改進前后模型中,動態(tài)側(cè)壓力均遵循隨著深度的降低而減小的原則[21],改進后模型迭代1 080萬步后糧食卸空,對應物理時間為31.4 s,原模型迭代450萬步卸空,對應物理時間為29.1 s,改進模型的卸糧過程所迭代的步數(shù)較多說明卸糧經(jīng)歷時間比較長,出流相對比較滯澀。
圖6 ball單元模型動態(tài)側(cè)壓力模擬結(jié)果
圖7 改進顆粒簇模型的動態(tài)側(cè)壓力模擬結(jié)果
取1號監(jiān)測點為例進行分析,如圖8所示,從卸料開始到結(jié)束,原模型和改進顆粒簇模型以及室內(nèi)試驗結(jié)果的側(cè)壓力變化趨勢對比。經(jīng)過對比可以發(fā)現(xiàn),改進前后模型相較試驗結(jié)果雖然波動都比較劇烈,但是整體趨勢基本一致,且數(shù)值圍繞實驗值波動,其他監(jiān)測點也呈現(xiàn)類似規(guī)律,由此可以進一步驗證此模擬試驗的真實性。改進后模型的動態(tài)側(cè)壓力波動幅度比較小,波動趨勢較為平緩,這與筒倉卸料過程中改進模型相較于原模型摩擦更充分,顆粒之間的“自鎖現(xiàn)象”有關。
圖8 卸料過程中動態(tài)側(cè)壓力對比
模型誤差可以通過公式(24)計算。
(24)
注:pA為原模型或改進顆粒簇模型的模擬結(jié)果;pB為試驗結(jié)果。
計算可得,改進模型誤差為11.28%,原模型誤差為29.175%,證明改進模型可以有效地提升卸糧動態(tài)側(cè)壓力模擬的精確度。各點動態(tài)側(cè)壓力最大值如表 3所示,可以證明改進模型更符合試驗結(jié)果。
表3 動態(tài)側(cè)壓力最大值對比
以改進模型為例,在筒倉下部建立測量圓,對比改進顆粒簇單元和傳統(tǒng)的ball單元的細觀結(jié)構(gòu)參數(shù),分析兩者細觀力學行為。
利用PFC提供的measure命令生成測量圓[22],如圖9所示。
圖9 測量圓位置
5.1.1 孔隙率對比
根據(jù)孔隙率的定義[23],可以知道,孔隙率的確定方法為:
(25)
式中:Vp為筒倉模型中顆粒的體積;V為模型中筒倉區(qū)域的體積。可以得到圓球顆粒和橢圓顆粒計算體積的表達式為:
Vpb=nb×R2π
(26)
Vpc=nc×(Ra×Rb×π)
(27)
則原單一圓球顆粒模型的孔隙率表達式為:
(28)
改進后模型孔隙率表達式為:
(29)
筒倉在靜態(tài)儲糧狀態(tài)下,根據(jù)室內(nèi)試驗結(jié)果,小麥在糧倉儲存狀態(tài)下,孔隙率為40.04%,觀測模型中九個measure圓的平均孔隙率,本研究選用2D模型,由孔隙率轉(zhuǎn)換關系[24]:
ε3d=1-ξ(1-ε2d)
(30)
(31)
式中:Dr為糧食相對密實度;ε2d為糧食二維孔隙率;ε3d為糧食三維孔隙率。
計算可得轉(zhuǎn)換得到三維狀態(tài)下的改進顆粒簇模型孔隙率為39.13%,誤差2.27%,而單一ball單元模型孔隙率為42.20%,誤差5.39%。改進模型在儲糧孔隙率方面下降了3.07%,儲料更加密實,精確度提升了3.12%。
5.1.2 配位數(shù)對比分析
配位數(shù)也是重要的細觀結(jié)構(gòu)參數(shù)之一,與孔隙率一樣用以表征堆積顆粒的密實度,意為某一顆粒附近與顆粒相接觸的顆粒數(shù)量,配位數(shù)與孔隙率成反比[25]。原模型配位數(shù)為3.67,改進后模型配位數(shù)為4.66,證明改進模型可以有效增大顆粒集合的配位數(shù)。在細觀結(jié)構(gòu)參數(shù)方面,改進模型可以使糧食更密集。
通過對比兩種模型下的卸糧力鏈圖,比較分析其細觀力學動態(tài)變化,根據(jù)所得結(jié)果可以清晰地發(fā)現(xiàn),原模型接觸力數(shù)量為31 279,改進后模型的接觸力數(shù)目為61 403,可以得知改進后模型的力鏈更加密集,顆粒間的結(jié)構(gòu)網(wǎng)絡更加復雜。根據(jù)圖10所示的兩模型局部力鏈圖,改進后模型的力鏈更密集,接觸力分布更均勻。由于顆粒的自鎖現(xiàn)象[26],clump單元之間易形成堅實的顆粒鏈抵御外力作用。改進的顆粒簇單元模型中可以清晰地看到數(shù)條橫向的拱形力鏈,更加符合卸糧中瞬時拱此消彼長的實際情況[27]。
圖10 卸糧模擬力鏈圖
以室內(nèi)卸糧物理模型試驗為基礎,采用改進顆粒簇單元進行數(shù)值模擬試驗,并和傳統(tǒng)采用單一ball單元的模擬結(jié)果進行對比分析。研究表明顆粒形態(tài)不僅影響筒倉卸糧過程中的宏觀力學行為,對顆粒間的細觀結(jié)構(gòu),流動過程中顆粒間的力鏈動態(tài)演化過程都有明顯影響。
從宏觀角度來看,改進顆粒簇模型使得倉壁卸糧動態(tài)側(cè)壓力波動更平緩,卸糧過程持續(xù)時間更長。通過對比模擬與試驗卸糧過程中所得動態(tài)側(cè)壓力結(jié)果,改進模型誤差為11.28%,原模型誤差為29.175%,改進顆粒簇模型相比原模型降低了17.895%的誤差。
顆粒細觀結(jié)構(gòu)方面,筒倉在靜態(tài)儲糧狀態(tài)下,原模型孔隙率為42.20%,而改進模型孔隙率為39.13%,改進模型的顆粒堆積孔隙率下降了3.07%,精確度上升了3.12%,原模型配位數(shù)為3.67,改進后模型配位數(shù)為4.66,配位數(shù)增大了27%,說明改進顆粒簇模型可以有效地增大顆粒堆積的密實度,更符合實際糧食顆粒接觸情況。
糧食顆粒形態(tài)影響顆粒的細觀力學行為,模擬結(jié)果表明改進模型接觸力數(shù)量和傳統(tǒng)ball模型相比增加了30124,接觸力鏈分布更均勻密集,改進的顆粒簇單元模型中可以清晰地看到數(shù)條橫向的拱形力鏈,更能清晰反映出瞬時拱此消彼長的動態(tài)變化規(guī)律。
針對筒倉卸糧方面的改進顆粒簇模型,有效提高了模擬的精確度,不僅對日后糧倉的設計提供了參考,而且對于進一步模擬不規(guī)則糧食顆粒在倉內(nèi)的流動機理也具有重要借鑒價值。