考慮顆粒形狀的離散元法在巖土工程中的應用進展

張世民, 陳萬里, 劉嘉英, 孫宏磊

(1.浙江工業大學土木工程學院, 杭州 310014; 2.浙大城市學院工程學院, 杭州 310015; 3.城市基礎設施智能化浙江省工程研究中心, 杭州 310015)

巖土工程中往往涉及大量的散粒體材料,散粒體物理力學性質復雜,以往常用的宏觀連續體力學理論體系分析其應力應變特征,但這種傳統方法忽略了散粒體的離散性,對散粒體內部作用機理認識不足,較難描述應變軟化、應力路徑相關性等特征,致理論與實際相離。為探究巖土類散粒體力學特征的內部作用機理,節約研究成本,提高研究效率,離散元法(discrete element method, DEM)[1]由Cundall等[2]于1979年提出,并逐漸在巖土數值分析中得到了應用。

離散元法的基本思想為將整個介質材料看作一系列離散的單元,每一單元顆粒均可以獨立運動并遵從牛頓運動定律,對每一單元求解其運動方程以進一步迭代計算。通常而言,構成散粒體的顆粒,如巖土工程中的粗粒土、砂土等,具有不規則的幾何形狀,因此在離散元模擬中,也應考慮單元顆粒的幾何形狀特性。

早在20世紀30年代,就有學者認識到顆粒形狀會對土體力學特性產生一定程度的影響[3-4],土體顆粒由于形狀所帶來的摩擦咬合,表面凹凸所產生的嵌鎖作用都是提供土體抗剪強度不可忽略的重要因素。真實工況中,由于介質材料的幾何形狀的不規則性,現場試驗較難精確還原材料形狀參數,許多學者提出了顆粒形狀的離散元模擬方法:①以簡單形狀或形狀參數為分類依據的典型形狀模擬法[5-7];②基于三維掃描重現顆粒輪廓的實物掃描法[8-9]。采用這些方法進行離散元試驗,試驗結果均表明了顆粒形狀對散粒體材料力學特性產生影響[10-12]。然而大多數的研究都受限于模擬形狀方法的局限性,還有待進一步完善。

基于目前巖土工程計算模擬領域離散元法應用的現狀,現細致介紹離散元法的模擬方式及模擬成果,分析各模擬顆粒形狀的方法的優劣并對各方法的實際應用性進行評價,最后對考慮顆粒形狀離散元模擬方法進行展望。

1 離散元方法應用現狀

1.1 細觀機制分析

采用離散元法進行散粒體物理力學性質模擬,很大一部分是用來分析顆粒材料復雜力學響應的細觀機制。目前常對配位數、微觀摩擦角、組構張量等進行分析,進而建立宏微觀聯系。對散粒體材料力學性質的機制分析目前還處于定性階段,部分學者提出的定量關系的確定方法具有一定的參考意義,為本構模型的建立提供了思路。楊升等[13]模擬直剪試驗分析了顆粒初始配位數與內摩擦角的關系以及不同顆粒試樣的體積變化規律,從細觀上描述了不同形狀顆粒的空間活躍程度。張玲等[14]對碎石和黏土顆粒進行數值三軸試驗分析了筋材開孔率、筋材抗拉剛度等因素對界面摩擦特性的影響。徐小敏等[15]模擬三軸試驗探究了初始楊氏模量和初始剪切模量與顆粒剛度比之間的關系,顆粒剛度比和初始泊松比之間的關系。劉勇等[16]根據粗粒土室內試驗結果和數值試驗對參數進行標定,研究了顆粒間摩擦系數、細觀連接強度以及孔隙率對宏觀力學特性的影響。王怡舒等[17]進行了不同滑動摩擦系數和滾動摩擦系數組合下的常規三軸模擬,研究了顆粒接觸摩擦對顆粒宏細觀力學特性的影響規律。王舒永等[18]利用三維掃描技術建立了不同含石量的土石混合體離散元模型,從微觀角度解釋了土石混合土剪切過程中一系列的變形并建立了宏細觀參數之間的線性方程。路德春等[19]確定了組構張量與應變張量的關系并用隱式方程的方式表達了出來,其中定量關系的確定方法為建立大尺度的本構提供了一種潛在思路。張家偉等[20]通過離散元法生成不同數量和角度的裂縫巖,定量分析了裂隙數量和巖石損傷變量之間的關系。閆洪超等[21]基于Hertz接觸理論模擬砂土直剪試驗推導出具有普適性的砂土宏細觀力學關系的模型,為定量構建土的本構關系奠定了基礎。

1.2 復雜應力路徑模擬

散體類巖土材料在實際中常經歷復雜應力路徑作用,常規物理試驗無法還原許多客觀存在的三維復雜應力路徑,離散元數值試驗方法能夠簡便地實現對復雜應力路徑的模擬,揭示復雜應力路徑下散粒體力學性質的演化規律。薛龍等[22]通過對3個主應力大小和方向的任意控制實現復雜應力路徑,顯示了數值試驗方法在深入研究三維復雜應力路徑下粒狀介質力學響應方面所擁有的能力和優勢。劉嘉英等[23-24]進行真三軸應力路徑模擬試驗分析了加載過程中宏微觀力學參數的演化過程以及組構張量與應力張量的聯系。李濤等[25]為探究結構性土粒間膠結的關鍵特征,建立了結構性土離散元模型開展多種應力路徑下的離散元數值分析,再現了結構性土體的主要宏觀力學性質。郝曉平[26]研究了砂巖試樣在不同應力路徑下宏觀應力-應變響應及細觀接觸力的分布規律發現試樣破壞時的峰值應力及體積應變與圍壓有強相關性。

1.3 多尺度耦合分析

為研究環境及地質災害、第三方施工行為等復雜載荷擾動下基坑、管道、樁體等宏觀變形與散粒體顆粒的細觀力學行為,需要考慮基坑、管道、水體等與周圍土體的耦合作用,鑒于離散元法在土體分層、壓縮、運移、流變等大變形領域時有其獨特的優勢,可采用離散元與其他數值方法耦合這樣不僅提高了計算效率,還提升了模擬精度,進而能夠更好地應用于巖土工程之中。徐濤龍等[27]利用SPH-FEM算法和PFC-FLAC耦合分析技術對特定管土耦合場景成功實踐,豐富了土中多物理場耦合領域。高濤等[28]采用離散元-有限元差分法的跨尺度耦合進行樁和周圍土體接觸過程中的穩定性分析,詳細解釋了下沉過程中周圍土的細觀變形、應力分布和樁體自身的變化及原因。禹海濤等[29]采用有限差分法(finite difference method,DEM-FDM)方法研究隧道開挖誘發近斷層錯動的物理機制。楊江坤等[30]基于DEM-FDM耦合方法構建了三維分離式霍普金森壓桿(split hopkinson pressure bar,SHPB)沖擊數值模型為細觀角度理解巖石動態沖擊破壞提供了新的思路。李偉一等[31]利用CFD-DEM研究間斷級配砂土的滲流現象對土體抗剪強度的影響,為從本源上理解砂土滲透侵蝕做出了重要貢獻。熊書春等[32]基于CFD-DEM提出了優化算法使用粘結顆粒模型提升了計算效率還解決了三維情況下非球形大顆粒在流場中的運動模擬難題。徐汪豪等[33]使用剛體動力學(rigid body dynamics, RBD)-離散元(DEM)耦合的數值模擬方法,實現了在數值模擬中的滾刀動態被動轉動,進一步貼近工程實際。

2 考慮顆粒形狀的離散元方法

自然界中的顆粒形狀并非絕對的圓球形,因此在顆粒材料離散元數值分析中,已有關于顆粒形狀的研究[5-9]。從目前已有的模擬顆粒形狀方法來看大致可以分為兩大類:一類是模擬顆粒的典型形狀,一類是模擬顆粒的真實形狀。

2.1 模擬典型顆粒形狀

模擬典型顆粒形狀是目前的主流模擬方法,其方法為設置特殊形狀顆粒替換圓形顆粒以達到模擬典型形狀的目的,起初只能較為籠統地表達顆粒的形狀,后發展為利用形狀參數對顆粒形狀進行量化,對顆粒形狀量化再與圓形顆粒作對照能簡便、高效地探究顆粒形狀與宏細觀參數之間的內部機理,適用于機制分析。John等[5]首先嘗試采用數值模擬方法模擬不同的顆粒形狀但僅局限于籠統得表達顆粒形狀如較圓、不規則形狀等,隨后Wadell[3]提出以球形度n(表征橢圓顆粒、多面體顆粒的形貌參數,形貌上越接近球的顆粒,其球形度越接近于1)。來表示顆粒形狀與圓球形的偏離程度。在此基礎上劉清秉等[34]提出以球形度為關鍵形狀特征參數生成試樣。定量分析了球形度對土體力學特性的影響。

宋二祥[7]參考實際工程中堆石料的幾何特性將堆石料顆粒形狀大致分為4類,利用粘合法將半徑相同的球顆粒粘結成密排六方結構并按照單顆粒的精細程度及計算效率分別生成不同球度的試樣。分別采用13、23、33顆粒簇作為球形度0.878、0.617和0.446的顆粒模型并保證各組試樣的圓度和表面粗糙度相同。顆粒之間采用平行黏結模型無需修正密度。在參考已有研究的基礎上采用最大與最小顆粒粒徑路徑之比為1.25,試樣邊長與顆粒直徑之比在5～9。對不同顆粒形狀的試樣進行常規三軸試驗分析了峰值摩擦角、平均接觸力等微觀變量的變化規律。

周光軍等[35]研究礫石含量以及礫石長寬比對砂-礫石混合物宏觀力學特性的影響。顆粒形狀的量化方法經歷了一系列的發展逐漸變得成熟,但量化參數卻缺乏工程統一的應用標準。有不少學者采用長細比[36]、磨圓度[37]、軸向系數[38]、扁平度[39]等參數對顆粒形狀進行表征取得了不少研究成果但不同形狀量化參數適用的工況缺乏統一的標準,導致同一工況下不同學者研究不同的顆粒形狀參數,彼此之間的聯系與參考價值有所影響長此以往不利于此領域的發展。

王家全等[40]在保證生成的Clump顆粒體積不變的前提下,在三個大小相同的圓形顆粒中間塞入一個可改變直徑D的圓形顆粒控制生成的Cump顆粒的球形度n。圖1為Clump顆粒生成示意圖,圖2為球形度示意圖。林呈祥等[41]運用類似方法對月壤顆粒進行量化統計后研究顆粒凹凸度對月壤抗剪強度指標的影響。魏婕等[42]通過Clump命令預生成幾種形狀的顆粒在0～360°方位角內隨機分布來模擬不同形狀的顆粒和圓形顆粒之間的區別。楊升等[13]通過Clump命令生成橢圓形顆粒和啞鈴型顆粒來分析不同顆粒的配位數和體積變化及空間轉動特性也取得了較好的結果。利用Clump方法模擬顆粒典型形狀的方法在定量生成所需量化參數具有便捷性,在研究顆粒形狀對土體力學特性的影響上直觀得出結論,便于機理研究,但實際工程中參數標定困難、計算量大、誤差不可控等原因都是實際工程中應用所必須克服的難題。

圖1 Clump顆粒生成示意圖[40]Fig.1 The schematic diagram of clump generation[40]

圖2 Clump顆粒示意圖[40]Fig.2 The schematic diagram of particle shape[40]

2.2 模擬顆粒真實形狀

顆粒真實形狀模擬能夠真實反映散粒體實際情況,模擬顆粒真實形狀首先需基于掃描設備對顆粒外輪廓進行掃描,其次利用圓形顆粒填充外輪廓和利用數學方程再現顆粒形貌來達到模擬真實形狀的目的,這是目前模擬顆粒真實形狀的主流方法。

2.2.1 顆粒真實形狀掃描

純利用數學方程或顆粒變形對顆粒真實形狀建模十分困難,所以常借助成像設備來對顆粒真實形狀重構。二維真實形狀模擬常用到顯微鏡或高清照相機并通過圖像處理、數學分析方法計算顆粒真實形狀二維形狀參數。陳海洋等[43]利用光學顯微鏡獲取鈣質砂的二維幾何形狀成功對鈣質砂的幾何形態進行描述,研究發現鈣質砂顆粒形狀具有明顯分形特征,與顆粒粒徑密切相關。為進一步提高重構形狀的準確性,近年來X射線計算機斷層掃描技術(computed tomography, CT)成為研究顆粒三維形態特征的主要手段。當X射線穿透物體后,其強度會隨物質成分、密度、光源強度等變化,而CT掃描過程則是基于X射線的變化來生成物體的三維影像,將試件置于X射線下旋轉360°完成CT掃描過程[44]。早期的三維掃描技術主要用于研究三維剪切破壞過程中顆粒的運動以及接觸行為[42]。由于所獲得的CT影像體元像素具有鋸齒狀的特性,近期高精度CT掃描技術也逐漸應用于掃描顆粒形態特征和微觀結構[8-9],通過后期圖像處理技術可精確計算體積、表面積等三位尺寸[45]。但是邊界搜索算法具有較大的誤差,進而使得砂粒表面積和表面曲率等形態參數的精確測定成為一個巨大的挑戰[46]。崔博等[47]運用三維激光掃描技術提出礫石形狀的評定標準,通過三軸壓縮模擬試驗分析礫石土的力學性質。徐文杰等[48]根據現場試驗反演得到土石混合體參數,為離散元法在實際巖土工程中的應用提供了新的思路。

2.2.2 真實顆粒形狀建模

Matsushim等[49]和Katagiri等[50]在一定區域內設置數個大小不同的圓形顆粒,讓顆粒之間互相運動變形,最后對穩定下來的各個圓形顆粒黏結成為一個整體,這樣就完成了對實際形狀的模擬。但是由于起始圓的位置具有隨機性,所以試驗會有一定誤差,同時無法量化生成的顆粒形狀,總體而言此模擬方法效率較低,工程應用困難。Wang等[51]在此方法的基礎上在基本圓相鄰較多的區域用一個更大的基本圓來代替相鄰的基本圓,提高了模擬效率,此種方法不僅可以模擬平面形狀還可以對空間立面進行模擬。薛明華等[52]采用Clump模塊生成真實天然砂樣模型進行直剪試驗模擬揭示顆粒形狀對砂性土宏細觀力學性質的影響。為了優化計算Ferellec等[53]設置了一系列參數來優化計算效率,但是由于不同環境和材料的優化參數差異較大導致此種方法較難推廣應用。Shamsi[54]采用此方法來研究晶粒形狀對于顆粒土力學行為的影響。雖然保證了模擬精度但參數選取較為繁瑣。李皋等[55]運用填充法對礫巖真實形狀建模并開展單軸抗壓模擬試驗,分析了顆粒形態和顆粒粒徑大小對礫巖力學性質的影響,更好地反映了礫巖的力學性質。袁斌等[56]基于等效橢球體改進棱度指數的計算公式可定量生成量化后的顆粒形狀,雖參數調整較為繁瑣但已經實現了實際中的應用,應用性大幅提升,在應用碎石鐵道方面有較大潛力。

徐文杰等[48]基于三維掃描設備利用無嵌套方式將相鄰球顆粒互相嵌套在一起來,在塊體三維幾何空間范圍內均勻填充一定比例的直徑在1D～1.5D(D為標準圓球顆粒的直徑)的顆粒,以塊體三維幾何網絡作為固定邊界,采用Yade進行計算并不斷調整球體顆粒半徑直至所有顆粒的接觸力在規定范圍之內。基于原位土體直剪試驗結果反演得到離散元細觀參數,通過數值模擬與現場試驗的剪切應力關系曲線驗證了此建模方法的可行性和合理性。所構建的塊石三維形態數據庫極大地優化了土石混合體模擬效率,也為在真實工況下模擬土石混合體顆粒真實形狀打下了結實的基礎。

邊學成等[57]建立計算機成像裝置利用三個互相垂直的高清攝像機獲取顆粒外輪廓,進行道砟顆粒材料直剪試驗。由于依賴三個互相垂直的高清攝像機獲取顆粒外輪廓,成像設備的精度會對模擬精度的產生的較大影響,此外僅靠三個互相垂直的成像設備會丟失顆粒表面特性,但誤差在允許范圍內。Reid等[58]運用LS-DEM提供的技術模擬三軸試驗,這種模擬試驗的預測水平較高,不僅能夠預測剪切帶的位置,還能預測剪切帶的方向、厚度和傾角。

部分學者運用數字方程進行真實顆粒形狀建模,此種方法模擬精度十分靈活,計算效率的優化較為成熟,采用數學方法來模擬顆粒形狀無疑為顆粒形狀的模擬提供了一種新的思路。Bowman等[59]首次使用數學傅里葉級數描述符來表征二維空間中顆粒的外輪廓,然而,真實顆粒的微觀形態遠比二維顆粒的外輪廓要復雜得多。Mollon等[60]將二維傅里葉擴展至三維利用三個真實的二維正交平面來構建顆粒的三維表面。然而,由于人工選擇原始剖面的三個正交截面,會丟失顆粒的表面特性尤其是十分不規則的顆粒。

Shi等[61]開發的球諧函數分析程序,用以表征顆粒的微觀形態,并以離散的方式確定介觀效應。此方法可根據所需的精細程度來選擇球諧函數的在再現階數。蘇棟等[62]獲取顆粒表面的三維幾何坐標,再通過球諧函數進行三維重構實現了對骨料顆粒的精確重構和表征分析。Cui等[63]通過球諧函數生成不規則塊石,用以研究土石混合體的力學性質。進行土石混合體的崩滑流試驗模擬結果與實驗測試結果吻合較好。利用球諧函數重構可對粗、細骨料進行模擬在土石混合體上的應用效果較好,計算精度靈活,計算效率優化較為成熟。李景哲等[64]基于高斯差分小波提出了重構骨料形狀的方法與球諧重構法相比可有效避免振鈴效應導致的重構誤差增大的現象。王嗣強等[65]發展了基于水平集方法的任意形態接觸算法用以準確計算單元間的接觸方向和重疊量,最大程度保留顆粒外輪廓幾何特性。

利用初始圓球顆粒運動和變形的方法由于初始圓顆粒的位置和大小具有隨機性,需反復調試,變形參數難以參考所以導致此種方法較難推廣應用。由于掃描設備成像具有鋸齒狀特性所以對精度要求較高并且模擬的過程十分耗時,雖有學者提出優化模擬效率的參數,但對于不同的環境和材料還需另外試驗選取新的優化參數導致應用性不強。從另一個角度,由于掃描設備的問題這反過來限制了試樣的尺寸,總的來說利用三維掃描設備模擬顆粒真實形狀的方法還需進一步發展,其中如何平衡模擬精度以及計算效率是其中的關鍵。

3 存在問題

盡管采用各類方法對不同顆粒形狀的巖土類散粒體已有一定的研究,但在考慮顆粒形狀進行離散元模擬及分析時,還存在以下問題。

(1)在工程適用性方面,目前缺乏統一的工程應用顆粒形狀分類標準,建議在充分認識顆粒形狀參數對土力學特性影響的基礎上構建出一套行之有效的顆粒形狀分類標準。在此基礎上進行大規模離散元數值分析,以期合理描述工程問題,采用離散元方法為工程問題提出合理的解決方案。

(2)進行顆粒材料真實形狀模擬時,成像設備的掃描精度要求較高,且利用數學函數模擬真實形狀對于使用者本身有門檻要求,不便于推廣應用。運用數字方程來模擬實際形狀的方法對內凹型的顆粒模擬效果不理想。

(3)參數標定一直是離散元模擬試驗的一大難點,無論是試錯法還是現場試驗法參數標定還是十分繁瑣。

(4)離散元法模擬真實形狀主要集中于塊石、礫石、卵石、土石混合體等材料,但真實形狀千萬種,單一種形狀就費時費力,實際工程中一個部位的模擬就需幾十甚至上百萬個顆粒單元,一次模擬就需數百個小時,若是提升計算效率模擬精度又無法有效保證,因此需要一種兼具效率與精度的廣泛適用性方法。

4 結論與展望

離散元法在巖土顆粒材料數值分析中已得到了一定的應用,尤其在機制分析、復雜應力路徑模擬和多尺度耦合研究方面。顆粒的形貌特征對土體力學性質產生影響已經得到了學者們的普遍認同,并取得一定的研究成果。

(1)離散元方法在巖土數值分析主要集中于細觀機制分析、復雜應力路徑模擬、多尺度耦合等方面,其中細觀機制分析為離散元巖土工程應用的熱點,但大部分還停留在定性分析階段,不少學者提出的定量分析方法具有一定的參考意義。

(2)模擬顆粒典型形狀由于模擬門檻較低,可便捷生成量化后的顆粒形狀有適用于機理研究。近年來模擬真實形狀已可應用小規模工程,但參數選取困難需、誤差不可避免等都是實際工程應用的難點。在提升計算效率方面可以從大數據神經網絡出發通過對算法的優化從根本上提升計算效率。

(3)散粒體的加卸載過程在巖土工程中無處不在,巖土顆粒材料在加卸載過程中的力學響應具有一定的特殊性,目前對于散粒體加卸載的研究主要集中于應力路徑,研究顆粒形狀對加卸載過程的影響仍然較少。在后續的研究中可以對這方面進行深入分析。