大規(guī)模顆粒系統(tǒng)的精確縮尺和粗?；x散元方法

趙婷婷， 馮云田

(1.太原理工大學(xué) 機械與運載工程學(xué)院，太原 030024；2.斯旺西大學(xué) 辛克維奇計算工程中心，英國斯旺西 SA1 8EN；3.寧波大學(xué) 沖擊與安全工程教育部重點實驗室，寧波 315211)

1 引 言

顆粒材料在自然界、工程應(yīng)用和日常生活中廣泛存在，由于其具有非連續(xù)、非均質(zhì)及各向異性的特性，使得以有限元方法為代表的傳統(tǒng)連續(xù)性數(shù)值計算方法無法準(zhǔn)確描述其力學(xué)行為。離散元方法[1]從20世紀(jì)70年代建立后不斷發(fā)展與完善，已成為探索顆粒材料物理力學(xué)性質(zhì)、解決不同領(lǐng)域工程問題的有效數(shù)值分析工具[2]。離散元方法的優(yōu)勢在于可以從微細觀尺度直接模擬顆粒之間的相互作用，進而反映顆粒系統(tǒng)的宏觀力學(xué)行為。在巖土工程領(lǐng)域，離散元方法可以描述典型巖土材料從微觀裂隙到宏觀破壞的全過程[3,4]；在工業(yè)工程領(lǐng)域，涉及到顆粒材料的儲存、混合、涂層以及運輸?shù)冗^程都可以用離散元方法來模擬[5,6]。離散元方法的優(yōu)勢也導(dǎo)致其在模擬工程尺度問題時會遇到計算資源不足的問題，當(dāng)采用真實顆粒粒徑和數(shù)量模擬實際問題時，現(xiàn)有的算法和計算機硬件水平難以有效支撐。真實系統(tǒng)的顆粒數(shù)量一般為萬億級別，現(xiàn)階段離散元模擬工作的顆粒數(shù)量通常為百萬至千萬的水平[7-9]，雖然目前已知單卡GPU已經(jīng)可以模擬1億規(guī)模的顆粒數(shù)量[10]，工程尺度應(yīng)用中面臨的超高計算量問題仍無法通過現(xiàn)有GPU技術(shù)有效解決。

相比于直接用大顆粒代替小顆粒的做法，近些年來已有很多研究者采用粗?；?Coarse-graining)理論來解決離散元方法在模擬工程尺度問題時計算量巨大的問題。粗粒化方法同樣使用較少數(shù)量的大顆粒代替系統(tǒng)中數(shù)量巨大的小顆粒，其出發(fā)點在于抓住大尺度的主要物理現(xiàn)象，對于小尺度的相對次要的物理現(xiàn)象可以平均化或者忽略。從問題的物理本質(zhì)來看，顆粒材料的宏觀力學(xué)行為主要由顆粒的集體行為決定，而不是單個顆粒獨自運動的軌跡，只要能保持顆粒材料的離散性質(zhì)，就可以反映其主要特性，這就為粗?；侄卧陔x散元方法中的應(yīng)用提供了基礎(chǔ)。為保證通過粗?；侄翁幚砗蟮念w粒集合能真實反映原始顆粒集合的物理特性，需要在兩系統(tǒng)間建立合理的等價關(guān)系。

國內(nèi)外研究者在粗?；碚撆c離散元方法結(jié)合方面開展了大量的研究工作，現(xiàn)有的粗?；椒〞ㄟ^縮放粒徑、調(diào)整參數(shù)及修改模型等手段使得放大后的顆粒集合體仍舊可以保持原始顆粒集合的性質(zhì)。文獻[11-13]提出了CG模型(Coarse Grain Model)用于流化床模擬，粗?；到y(tǒng)中的顆粒直徑是原始系統(tǒng)中顆粒直徑的h倍，粗?；到y(tǒng)中的一個顆粒代表原始系統(tǒng)中呈立方體規(guī)律排列的h3個顆粒。假設(shè)原始系統(tǒng)中呈立方體排布的顆粒集合中每一個顆粒的平動速度和轉(zhuǎn)動速度相同，且等于粗?；w粒的平動速度和轉(zhuǎn)動速度。CG模型還假設(shè)當(dāng)粗粒化顆粒發(fā)生二元碰撞時，原始集合體中的所有顆粒發(fā)生同步的二元碰撞，并將所有原始顆粒的接觸力進行疊加得到粗?；w粒的接觸力。在以上假設(shè)下，得到粗?；到y(tǒng)中接觸力的縮放系數(shù)為原始顆粒系統(tǒng)的h3倍。對于拖曳力和重力，同樣推導(dǎo)得到了h3的縮放系數(shù)。而對于范德華力，CG模型根據(jù)能量守恒的原則進行推導(dǎo)，得到縮放系數(shù)為h2。文獻[14-16]提出了CGSF(coarse-grained method for granular shear flow)用于模擬顆?；旌线^程中的剪切流動，同樣認(rèn)為粗?；w粒代表呈立方體規(guī)律排布的原始顆粒集合，該模型只針對顆粒剪切流動的情況，在推導(dǎo)粗粒化系統(tǒng)與原始系統(tǒng)之間四種能量守恒關(guān)系時(動能、彈性能和摩擦阻尼粘滯阻尼能量)，重點考慮顆粒的切向運動速度。在CGSF模型簡化顆粒排布幾何關(guān)系和運動形式的假設(shè)下，通過對滑動摩擦系數(shù)、線性剛度系數(shù)和恢復(fù)系數(shù)進行縮放來滿足能量守恒關(guān)系，得到的縮放規(guī)律較為復(fù)雜，特別是滑動摩擦系數(shù)的縮放關(guān)系還與粗粒化顆粒實時的角速度有關(guān)。CGSF模型對于滾筒中的顆?；旌线^程粗?；M的效果良好，但該模型缺乏廣泛的適用性。文獻[17,18]提出了SPA模型(Similar Particle Assembly)，將原始顆粒粒徑放大h倍得到粗粒化顆粒，該模型不再對粗?；w粒代表的原始顆粒排布做出假設(shè)，認(rèn)為粗?；到y(tǒng)的顆粒排布與原始系統(tǒng)相似。SPA模型對控制方程中各項的縮放規(guī)律缺乏嚴(yán)格的理論推導(dǎo)，通過假設(shè)粒徑對顆粒動力學(xué)行為具有決定性作用，直接將h3作為接觸力、液橋力和拖曳力的放大系數(shù)。此外，不同學(xué)者提出的粗?；Ｐ瓦€包括Imaginary Sphere模型[19]和Representative Particle模型[20]等。

可以看出，粗?；椒ㄔ陔x散元模擬中得到了越來越多的關(guān)注，但現(xiàn)有的粗?；Ｐ痛蠖嘀苯訌男枰M的問題出發(fā)，通過分析問題本身的特征提出一系列假設(shè)，進而得到粗粒化與原始系統(tǒng)的等價關(guān)系。通過這種方式得到的粗?；Ｐ停M管在特定應(yīng)用中取得了比較好的模擬效果，但很難推廣到其他問題中。而且由于假設(shè)的提出往往具有隨意性，使得無法對原始系統(tǒng)與粗粒化系統(tǒng)計算差距的產(chǎn)生原因以及規(guī)模進行進一步分析。文獻[21,22]從更一般的角度出發(fā)，提出了介于原始系統(tǒng)和粗粒化系統(tǒng)之間精確縮尺(Exact Scaling)系統(tǒng)，并且通過嚴(yán)格的理論推導(dǎo)，得到了在精確縮尺系統(tǒng)中顆粒集合各物理量應(yīng)該滿足的縮放關(guān)系。本文將在精確縮尺模型的基礎(chǔ)上，通過多尺度方法，建立粗?；到y(tǒng)和原始系統(tǒng)之間的縮放關(guān)系，得到離散元接觸模型中相關(guān)參數(shù)的縮放定律，并通過離散元算例進行驗證。

2 精確縮尺模型

原始系統(tǒng)、精確縮尺系統(tǒng)以及粗粒化系統(tǒng)之間的關(guān)系如圖1所示，為了便于說明，圖1中顆粒規(guī)則排列。原始系統(tǒng)和粗?；到y(tǒng)占據(jù)的幾何區(qū)域大小是相同的，粗?；到y(tǒng)中的顆粒直徑較原始系統(tǒng)中顆粒直徑放大一定的倍數(shù)；在精確縮尺系統(tǒng)中，顆粒直徑及幾何區(qū)域較原始系統(tǒng)同步放大相同的倍數(shù)，可以將粗?；到y(tǒng)看作精確縮尺系統(tǒng)的一個局部區(qū)域。需要說明的是，對于大尺度顆粒系統(tǒng)的模擬，精確縮尺方法和粗?；椒ú⒉皇莾煞N處于并列位置的方法，采用精確縮尺方法可以準(zhǔn)確推導(dǎo)出原始系統(tǒng)小粒徑顆粒集合體與放大后系統(tǒng)大粒徑顆粒集合體之間不同物理量的比例關(guān)系。但由于精確縮尺方法會將系統(tǒng)的計算區(qū)域同步放大，因而原始系統(tǒng)與精確縮尺系統(tǒng)的顆粒數(shù)量保持一致，從計算效率的角度來看，精確縮尺方法在大尺度顆粒系統(tǒng)的模擬方面不會帶來提高。但為了建立原始系統(tǒng)與粗?；到y(tǒng)之間的縮放關(guān)系，可以將精確縮尺系統(tǒng)作為橋梁，首先分析原始系統(tǒng)與精確縮尺系統(tǒng)之間各物理量需要滿足的比例關(guān)系。

圖1 原始系統(tǒng)、精確縮尺系統(tǒng)以及粗粒化系統(tǒng)

2.1 相似定律

由量綱分析可知，在只考慮物體的機械運動時，任意物理量q的量綱可由國際標(biāo)準(zhǔn)單位制下的基本變量組合長度[L]、質(zhì)量[M]和時間[T]推導(dǎo)得

[q]=LaLMaMTaT

(1)

其用向量形式的單位標(biāo)準(zhǔn)基本變量表示為

〈q〉=(aL,aM,aT)T

(2)

(3)

在由不同單位基本變量組合表示的系統(tǒng)中，物理量q′表示為

〈q′〉=R-1〈q〉

(4)

2.2 顆粒系統(tǒng)的縮放關(guān)系

(5)

在精確縮尺系統(tǒng)中，選取三個基本變量對應(yīng)的縮放系數(shù)，分別為長度h、時間h和 密度1，即基本單元轉(zhuǎn)換系數(shù)為Hb={h,h,1}，理論上三個基本變量的縮放系數(shù)可以任意選取，目前的取值組合可以為解釋原始系統(tǒng)與精確縮尺系統(tǒng)之間的等價關(guān)系帶來方便。由轉(zhuǎn)換矩陣的逆R-1及縮放系數(shù)Hb即可推得精確縮尺系統(tǒng)中任意物理量對應(yīng)的縮放系數(shù)

(6)

以顆粒系統(tǒng)中的物理量力F為例，在由標(biāo)準(zhǔn)變量組合表示的原始系統(tǒng)中，其量綱為

[F]=N=LMT-2,〈F〉=(1,1,-2)T

(7)

〈F′〉=R-1〈F〉=(4,-2,1)T

(8)

由式(6)可以得到在精確縮尺系統(tǒng)中F的縮放系數(shù)為

hF=Hb·^〈F′〉=h4·h-2·11=h2

(9)

根據(jù)以上推導(dǎo)，可以得到精確縮尺系統(tǒng)中各物理量的縮放系數(shù)，部分主要物理量列入表1，詳見文獻[24]。選擇當(dāng)前的基本單元組合以及對應(yīng)的縮放系數(shù)，可以使精確縮尺系統(tǒng)中的應(yīng)力、應(yīng)變、動能密度以及應(yīng)變能密度與原始系統(tǒng)相等，保證了精確縮尺系統(tǒng)與原始系統(tǒng)之間的等價關(guān)系。

表1 精確縮尺系統(tǒng)中部分物理量的縮放系數(shù)

在離散元計算中，精確縮尺系統(tǒng)與原始系統(tǒng)之間的等價關(guān)系可以通過兩種方式實現(xiàn)，(1) 文獻[24]詳細討論了對于離散元方法中不同種類的接觸模型，通過保證其在應(yīng)力應(yīng)變形式下的表達式與縮放系數(shù)無關(guān)對接觸參數(shù)進行處理；(2) 將離散元計算中涉及到的物理量完全按照量綱對應(yīng)的縮放系數(shù)進行放大或縮小。以上兩種方法在物理上是等價的，可以根據(jù)在已有離散元程序中實現(xiàn)的便捷程度進行自由選擇。

由于精確縮尺模型的時間變量較原始系統(tǒng)放大了h倍，對應(yīng)的時步Δt也會較原始系統(tǒng)放大h倍，在采用中心差分法進行時間積分時，兩系統(tǒng)需要的時步數(shù)是相同的，也就是說采用精確縮尺系統(tǒng)代替原始系統(tǒng)并不能從計算效率上帶來任何提高。提出精確縮尺系統(tǒng)的作用在于，可以從更一般的角度對粗?；到y(tǒng)中顆粒層面相關(guān)物理量的處理(包括接觸模型和顆粒相對速度等)給出可解釋的理論依據(jù)。

3 粗粒化模型

粗?；到y(tǒng)在放大顆粒粒徑的同時，保持系統(tǒng)總體區(qū)域與原系統(tǒng)一致，粗?；椒〞p少顆粒數(shù)量，粗粒化系統(tǒng)與原始系統(tǒng)之間不再具有幾何相似性，無法精確重現(xiàn)原始系統(tǒng)的物理性質(zhì)。但在精確縮尺系統(tǒng)中得到的相似定律可以應(yīng)用于粗?；到y(tǒng)，保證粗?；到y(tǒng)計算結(jié)果具有較高的精度。

取原始系統(tǒng)及粗?；到y(tǒng)中相同幾何區(qū)域的顆粒為研究對象，將圖2中黑色橢圓內(nèi)的顆粒集合看作代表性體積單元(RVE)。為保證粗?；到y(tǒng)的離散元計算結(jié)果能重現(xiàn)原始系統(tǒng)的物理力學(xué)性質(zhì)，兩系統(tǒng)對應(yīng)的RVE需要滿足幾何一致、質(zhì)量、動量以及能量的近似守恒等條件。

圖2 原始系統(tǒng)及粗?；到y(tǒng)中的代表性單元

(10)

在RVE邊界上的顆粒接觸數(shù)量滿足

(11)

(12)

動量守恒的條件要求滿足

(13)

能量守恒包括動能守恒、應(yīng)變能守恒以及能量耗散速率守恒三個方面，動能守恒可以由顆粒質(zhì)量以及速度之間的關(guān)系推導(dǎo)得到。

根據(jù)圖3所示顆粒系統(tǒng)中代表性單元的接觸受力情況，兩系統(tǒng)中平均柯西應(yīng)力的表達式[23]為

圖3 顆粒系統(tǒng)代表性單元受力分析

(14)

根據(jù)式(11)，可以得到兩系統(tǒng)中RVE單元邊界接觸力的縮放關(guān)系為

(15)

對于原始系統(tǒng)中的RVE單元，其總體平動控制方程為

(16)

粗粒化系統(tǒng)中RVE單元的平動控制方程為

(17)

已知在精確縮尺系統(tǒng)中力的縮放系數(shù)為h2(式(9))，與粗粒化系統(tǒng)細觀顆粒尺度力的縮放系數(shù)相同，進一步分析可以發(fā)現(xiàn)，精確縮尺模型中提出的對于不同種類離散元接觸模型的處理完全適用于粗?；到y(tǒng)的離散元計算，對于離散元計算中涉及到的無量綱系數(shù)(摩擦系數(shù)、泊松比和阻尼系數(shù)等)不需要做任何縮放。可以看出，原始系統(tǒng)和粗?；到y(tǒng)之間存在兩種尺度的縮放關(guān)系，即雙尺度粗?；?Two -scale coarse graining)，細觀顆粒層面相關(guān)物理量的縮放關(guān)系與精確縮尺模型得到的結(jié)果相同(如接觸力的縮放系數(shù)為h2)，宏觀顆粒集合層面相關(guān)物理量遵循另外一種縮放關(guān)系(如重力和拖曳力的縮放系數(shù)為h3)。

需要說明的是，以上原始系統(tǒng)與粗?；疪VE單元的控制方程只考慮了平動情況。對于轉(zhuǎn)動情況，由于系統(tǒng)總體轉(zhuǎn)動能無法簡單地直接將各顆粒單元的轉(zhuǎn)動能求和得到，所以轉(zhuǎn)動相關(guān)物理量的縮放規(guī)律更為復(fù)雜。根據(jù)目前的研究，需要根據(jù)不同算例中顆粒的實際運動情況具體分析轉(zhuǎn)動能的產(chǎn)生原因，進而得到對應(yīng)的轉(zhuǎn)動相關(guān)物理量的縮放系數(shù)。

對于時間變量，粗粒化系統(tǒng)和原始系統(tǒng)之間同樣存在兩個不同尺度的縮放關(guān)系，宏觀層面物理時間的縮放系數(shù)為1，細觀層面顆粒松弛時間的縮放系數(shù)為h，這就為離散元計算效率的提高帶來了極大的好處。當(dāng)顆粒粒徑放大h倍時，顆粒數(shù)量減少h3倍使得計算效率提高h3倍，計算時步減少h倍又可以使計算效率提高h倍，故粗粒化系統(tǒng)的計算時間是原始系統(tǒng)計算時間的1/h4。

4 算例分析

兩個算例證明，由精確縮尺模型推導(dǎo)得到的離散元接觸模型相關(guān)參數(shù)的縮放關(guān)系同樣適用于粗?；到y(tǒng)的離散元計算。算例的不同工況都采用的是粗?；椒?，即保持顆粒集合的宏觀幾何尺寸不變，只將顆粒粒徑進行放大。粗?；Ｐ屯ㄟ^保證任意兩顆粒單元之間接觸相似進而得到顆粒集合整體力學(xué)行為的相似，對于任意級配的顆粒集合都是適用的，為了便于比較不同顆粒粒徑對應(yīng)的計算結(jié)果，本文采用了單粒徑的顆粒集合形式。

4.1 筒倉側(cè)壁壓力計算

建立如圖4所示的筒倉模型，原始系統(tǒng)中的顆粒粒徑為5 mm，將顆粒粒徑放大2倍和3倍，分別采用線性模型和赫茲模型用于接觸力的計算。由文獻[22]可知，對于三維離散元計算，線性接觸模型是尺度相關(guān)模型，當(dāng)系統(tǒng)中顆粒粒徑放大h倍時，需要將剛度系數(shù)K同樣放大h倍用于計算；而赫茲模型是尺度無關(guān)模型[22]，用于不同尺度系統(tǒng)計算時，不需要改變接觸參數(shù)。計算表2所列的8種不同工況，觀察不同接觸系數(shù)對筒倉側(cè)壁壓力計算結(jié)果的影響。

圖4 筒倉模型及原始顆粒集合

表2 筒倉側(cè)壁壓力計算工況

線性接觸模型對應(yīng)的計算結(jié)果如圖5(a)所示，當(dāng)剛度系數(shù)隨著粒徑變化放大相同的倍數(shù)時，筒倉側(cè)壁壓力的計算結(jié)果與原始系統(tǒng)的計算結(jié)果差距較小，當(dāng)接觸剛度系數(shù)保持不變時，粗?；到y(tǒng)的計算結(jié)果無法反映原始系統(tǒng)的力學(xué)性質(zhì)，說明精確縮尺系統(tǒng)推導(dǎo)得到的線性接觸模型參數(shù)轉(zhuǎn)換關(guān)系適用于粗粒化系統(tǒng)的離散元計算。對于赫茲接觸模型，由于其具有尺度不變性，在不同尺度粗?；到y(tǒng)的計算中可以選取和原始系統(tǒng)相同的接觸模型參數(shù)，最終得到的計算結(jié)果誤差較小，如 圖5(b)所示。

圖5 不同計算工況下的筒倉側(cè)壁壓力

4.2 休止角計算

為了驗證粗?；Ｐ椭屑氂^層面的顆粒間接觸力應(yīng)該滿足h2的縮放關(guān)系，采用粗?；椒ㄓ嬎悴煌w粒尺度對應(yīng)的考慮粘聚力的顆粒材料休止角，粘聚力的計算公式[24]為

Fa=F0(1-δ/D0)

(18)

式中F0為最大引力，D0為引力范圍，δ為兩顆粒的重疊距離。

首先在圓筒內(nèi)生成顆粒試樣，將圓筒緩慢提升，顆粒在重力作用下自由滑落，形成穩(wěn)定結(jié)構(gòu)用以計算顆粒材料的休止角。將原始系統(tǒng)中的顆粒粒徑分別放大2倍和3倍，將粘聚力按照h，h2和h3的縮放系數(shù)進行放大，具體工況列入表3，將計算得到的最終休止角與原始顆粒系統(tǒng)的休止角進行對比，如圖6所示。圖6(a)為當(dāng)粒徑的縮放系數(shù)hd=2時，三種粘聚力縮放系數(shù)hf對應(yīng)的計算結(jié)果，可以看出，在粗?；到y(tǒng)的縮放比例不大時，粘聚力采用平方關(guān)系或者立方關(guān)系進行縮放，都能得到與原始系統(tǒng)相近的結(jié)果。當(dāng)粒徑進一步擴大，hd=3時的計算結(jié)果如圖6(b)所示，可以看出，當(dāng)粘聚力按照精確縮尺模型提出的相似定律擴大h2倍時，最終得到的休止角與原始系統(tǒng)差距最小。將不同工況對應(yīng)的計算時間同樣列入表3，可以看出，粗?；到y(tǒng)的計算時間近似等于原始系統(tǒng)計算時間的1/h4，證明了第3節(jié)對粗?；到y(tǒng)計算效率提高的結(jié)論。

表3 休止角計算工況

圖6 不同計算工況下的休止角

5 結(jié) 論

本文針對離散元模擬大規(guī)模顆粒系統(tǒng)時涉及到的縮放問題，建立了一種雙尺度縮放體系，并通過不同算例進行了驗證。采用量綱分析的方法，得到了顆粒系統(tǒng)各物理量在原始系統(tǒng)及精確縮尺系統(tǒng)之間的縮放關(guān)系，為離散元接觸模型中接觸參數(shù)的處理提供了理論依據(jù)。采用多尺度描述方法，建立了粗?；到y(tǒng)與原始系統(tǒng)的代表性單元(RVE)，根據(jù)不同系統(tǒng)RVE單元之間質(zhì)量守恒、動量守恒以及能量守恒關(guān)系，得到粗?；到y(tǒng)與原始系統(tǒng)之間宏觀和細觀兩種不同尺度的縮放關(guān)系。在細觀顆粒層面上，由精確縮尺模型得到的相似定律完全適用，即離散元接觸模型接觸參數(shù)應(yīng)該按照文獻[22]提出的方法進行處理，可以由筒倉側(cè)壁壓力及休止角的離散元計算結(jié)果進行驗證。粗?；碚摽梢詾楣こ坛叨却笠?guī)模離散元計算提供最有效的解決辦法，顆粒粒徑擴大h倍，離散元計算效率提高h4倍。但目前還沒有理論化的方法可以預(yù)測粗粒化系統(tǒng)與原始系統(tǒng)之間的計算誤差，需要在今后的工作中展開進一步的研究。此外，在今后的工作中，還需要對轉(zhuǎn)動相關(guān)物理量的縮放定律進行更深入的研究。