程宏旸， Thomas Weinhart

(特文特大學(xué) 工學(xué)院，荷蘭恩斯赫德 7500AE)

1 引 言

顆粒材料如沙土、巖石、藥物粉末和農(nóng)業(yè)谷物等廣泛存在于自然界和工業(yè)應(yīng)用中。這些材料均由離散的可自由運(yùn)動(dòng)的顆粒集合而成。因此，其宏觀尺度力學(xué)性質(zhì)由顆粒間相互碰撞及微-細(xì)觀結(jié)構(gòu)演變決定[1,2]。

顆粒材料微觀、細(xì)觀結(jié)構(gòu)及其空間尺度隨外部荷載變化而演變，即使僅考慮固態(tài)顆粒介質(zhì)(如準(zhǔn)靜態(tài)下的剪脹[3]、塑性屈服[4]、組構(gòu)張量和各向異性[5,6])，其宏觀尺度上的建模也可能是相當(dāng)復(fù)雜的。盡管當(dāng)前宏觀模型可以很好地描述準(zhǔn)靜態(tài)下復(fù)雜顆粒行為，流態(tài)顆粒介質(zhì)的剪切增稠、稀化[7]及分層[8]等流變現(xiàn)象的宏觀模擬還是非常困難，目前主要依靠顆粒及細(xì)觀尺度模型來進(jìn)行研究。

微-宏觀耦合模擬與單純的微觀或宏觀模擬相比具有明顯優(yōu)勢。通過建立合理的耦合模型，可以保持宏觀方法較高的計(jì)算效率和微觀方法較高的計(jì)算精度。本文主要探討同步多尺度方法(concurrent multi-scale methods),即微觀與宏觀子模型在同一個(gè)計(jì)算域上求解并在時(shí)間上同步。這需要與基于微-宏觀關(guān)系在物質(zhì)點(diǎn)和表征單元之間傳遞的分層多尺度[9-11]耦合方法作出區(qū)分。同步多尺度方法主要分為兩類,一是表面耦合，可用于模擬顆粒與可變形單元之間的相互作用；二是體積耦合，可用于模擬顆粒材料連續(xù)-非連續(xù)和流態(tài)-固態(tài)轉(zhuǎn)換等復(fù)雜物理現(xiàn)象。

以有限元(FEM)和離散元(DEM)為例(圖1)，表面耦合通常依賴于DEM模型中可變形構(gòu)件外表面與顆粒材料相接處的界面ΓC=ΓF E∩ΓD E。該界面接收來自FEM網(wǎng)格的運(yùn)動(dòng)信息，并且將顆粒與該界面間的相互作用力映射到FEM模型的牽引邊界上。由于表面耦合通常通過狄拉克函數(shù)來實(shí)現(xiàn)[12]，施加到FEM域上的耦合力場很可能是非光滑的，導(dǎo)致難以處理的累積誤差。目前FEM-DEM界面耦合方法的應(yīng)用，如輪胎、土工膜、土工格柵與土顆粒之間的相互作用[12,19,20]，均采用基于狄拉克函數(shù)的微-宏觀映射。

體積耦合又稱為Arlequin多模型耦合框架，由Dhia等[13,14]提出，最初用于耦合具有不同網(wǎng)格分辨率的FEM模型，由Wellmann等[15]首次擴(kuò)展到FEM-DEM耦合方法上。體積耦合在FEM-DEM模擬上的應(yīng)用主要通過在體積重疊區(qū)域ΩC=ΩF E∩ΩD E(圖1)上施加動(dòng)態(tài)約束，使微觀和宏觀兩個(gè)子模型的解盡可能一致[15]。而如何從微觀顆粒運(yùn)動(dòng)信息中提取適應(yīng)于FEM求解空間的宏觀場則是需要解決的難點(diǎn)之一。Wellmann等[15]使用有限元形狀函數(shù)構(gòu)建微-宏觀傳遞法則，將FEM-DEM體積耦合方法應(yīng)用在單樁安裝問題上,即在樁體周圍變形較大區(qū)域使用DEM，遠(yuǎn)離樁土作用范圍使用FEM，DEM域至FEM域的過渡區(qū)間則采用體積耦合方法描述。

圖1 同步多尺度模擬示例(ΩD E為非連續(xù)體顆粒材料，ΩF E1和ΩF E2分別為用連續(xù)體方法描述的可變形結(jié)構(gòu)與連續(xù)體顆粒材料)

本文采用粗?；?CG)方法[16]，將非連續(xù)性的DEM顆粒數(shù)據(jù)首先映射到在連續(xù)體上定義的宏觀場上，然后將這一宏觀場耦合到FEM網(wǎng)格上。粗?；僮饕肓朔Q為粗?；瘜挾?CG width)的自定義參數(shù)，為均勻化后的宏觀場定義了一個(gè)可調(diào)整的空間尺度?；谟纱至；茖?dǎo)出的表面耦合和體積耦合項(xiàng)公式，將探討如何通過調(diào)整CG寬度減小數(shù)值誤差；以彈性塊沖擊顆粒床和離散-連續(xù)介質(zhì)間的波傳播兩個(gè)數(shù)值算例，展示粗?；瘜挾葘?duì)同步多尺度耦合系統(tǒng)的能量守恒特性的影響，并討論如何在保持相同能量耗散的條件下獲得更高的計(jì)算效率。

2 基于粗?；腇EM-DEM耦合方法

在同步多尺度耦合框架內(nèi)，計(jì)算域Ω劃分為多個(gè)子域；每個(gè)子域中的材料行為由不同子模型描述。根據(jù)特定的耦合方法，耦合區(qū)域的定義也有所不同，如ΩC=ΩF E∩ΩD E適用于體積耦合，而表面耦合區(qū)間則定義為ΓC=?ΩF E∩?ΩD E，如圖1所示。本節(jié)將介紹如何使用CG定義空間上連續(xù)的速度場和表面牽引力場，以及如何使用這些均勻化的宏觀場進(jìn)行表面和體積耦合。上標(biāo)DE和FE表示離散元和有限元子模型，下標(biāo){α,β,γ,…}和{i,j,k,…}則用來區(qū)分DEM和FEM模型的變量。

2.1 微-宏觀傳遞法則

(1)

(2)

(3)

2.1.1 接觸力與均勻化后的宏觀牽引力

(4)

(5)

(6)

上述推導(dǎo)證明文獻(xiàn)[12,17,18]廣泛使用的表面耦合力公式為粗?；砻骜詈系囊环N特殊情況。

2.1.2 顆粒運(yùn)動(dòng)與均勻化后的宏觀速度場

2.1.1節(jié)給出了使用CG來均勻化顆粒-結(jié)構(gòu)之間接觸力的一般公式。對(duì)于顆粒材料整體來說，更多情況下需要通過細(xì)觀結(jié)構(gòu)和顆粒運(yùn)動(dòng)得到相應(yīng)的宏觀場。本節(jié)介紹使用CG將顆粒尺度的速度或位移均勻化至宏觀FEM模型上。

考慮Np個(gè)顆粒，與體積耦合域內(nèi)的有限單元e在空間上重疊。由式(2)可以在耦合區(qū)間內(nèi)任意位置處獲得在空間上連續(xù)的速度場vC G為

(7)

式中Np為所有耦合顆粒的個(gè)數(shù)。注意，由于粗粒度函數(shù)的有效范圍c有限，只有在核函數(shù)影響范圍內(nèi)的顆粒才會(huì)考慮在均勻化計(jì)算中。

(8)

(9)

使用式(9)定義宏觀速度場需要求得每個(gè)形狀函數(shù)Ψi在任意耦合顆粒位置xα處的值[15]，而式(8)只需要在高斯積分點(diǎn)上操作，后者的代碼實(shí)現(xiàn)顯然更加容易。

2.1.3 基于粗粒化的微-宏觀傳遞的一般性

2.2 控制方程與耦合項(xiàng)

為了方便處理由顆粒材料碰撞、沖擊和剛體運(yùn)動(dòng)產(chǎn)生的大位移，本文采用完全拉格朗日公式建立有限元模型的控制方程，以廣義胡克定律描述連續(xù)體材料的應(yīng)力應(yīng)變關(guān)系，時(shí)間積分則采用隱式Newmark-beta法。目的是獲得更高的精度來測試FEM-DEM耦合算法[19]。FEM和DEM子模型的時(shí)間步長相同，均為粒間碰撞時(shí)間tc的1/20。

顆粒尺度的力學(xué)行為則通過牛頓第二定律以及顆粒間的接觸力模型求解。接觸力模型通常包括法向/切向彈簧-阻尼系統(tǒng)，以及控制顆粒間剪切強(qiáng)度的摩爾庫倫準(zhǔn)則。由于接觸力模型需要精確地計(jì)算顆粒間接觸面積或重疊深度，離散元模擬的時(shí)間積分常采用顯示Leap-frog方案。

式(5,6)給出了從DEM顆粒施加到FEM網(wǎng)格的節(jié)點(diǎn)耦合力公式。此外，F(xiàn)EM網(wǎng)格更新的位置與速度也需要通過內(nèi)插傳遞到其在DEM域的副本上，即為顆粒提供位移邊界的一組相互連接的三角形剛體上。這些三角形平面與顆粒之間的接觸則采用常規(guī)的DEM接觸算法來計(jì)算。

(10)

(11)

3 數(shù)值算例

基于粗粒化的多尺度耦合與傳統(tǒng)方法相比的優(yōu)勢在于其可自由定義的空間尺度，使耦合物理場變得更加均勻，微-宏觀傳遞具有非局部(nonlocal)特性?；诖至；腇EM-DEM表面和體積耦合方法,可以通過單球在懸臂梁上的滑動(dòng)和離散-連續(xù)體材料模型內(nèi)部的波傳播算例進(jìn)行驗(yàn)證。本節(jié)將通過彈性塊沖擊顆粒床和體積耦合深度不同的懸臂梁兩個(gè)算例，研究耦合系統(tǒng)因表面或體積耦合產(chǎn)生或耗散的能量隨時(shí)間的變化，并調(diào)查粗?；瘜挾葘?duì)耦合系統(tǒng)能量守恒的影響。兩個(gè)數(shù)值算例的材料參數(shù)列入表1；FEM和DEM子模型的時(shí)間步長均為顆粒碰撞時(shí)間tc的1/20。

表1 有限元與離散元模型的材料參數(shù)

3.1 彈性體沖擊顆粒床

首先，通過顆粒在懸臂梁上受重力作用的滑動(dòng)，驗(yàn)證基于粗?；谋砻骜詈戏椒āＨ鐖D2所示，顆粒受重力在懸臂梁上的運(yùn)動(dòng)軌跡與考慮在不同接觸點(diǎn)施加集中點(diǎn)力的懸臂梁彎曲(為方便獲得解析解，此處懸臂梁不受重力)解析解一致。與常規(guī)局部耦合方法相比，采用粗?；砻骜詈戏椒ǐ@得的顆粒運(yùn)動(dòng)軌跡明顯更光滑，由此可以推斷粗?；彩诡w粒-單元間的耦合力場更光滑。

圖2 顆粒受重力在彈性懸臂梁上的運(yùn)動(dòng)軌跡

考慮如圖3所示的數(shù)值算例，檢驗(yàn)可變形彈性塊在重力作用下沖擊顆粒材料時(shí)，是否因耦合作用產(chǎn)生多余的能量。彈性塊由8個(gè)8節(jié)點(diǎn)六面體單元構(gòu)成，在t=0時(shí)從顆粒上5d處受重力自由下落，在t=0.04 s時(shí)沖擊顆粒床。耦合模擬開始之前，所有顆粒已經(jīng)在重力作用下松弛至準(zhǔn)靜態(tài)，顆粒系統(tǒng)的動(dòng)能-彈性勢能比Ek/Ep低于0.01。顆粒材料的側(cè)面為周期性邊界，底部是半無限剛性墻，耦合模擬初始狀態(tài)如圖3所示。

圖3 FEM-DEM表面耦合模擬的初始狀態(tài)

首先，考慮是否使用粗?；M(jìn)行微-宏觀傳遞，即粗?；瘜挾葹?和10R的情況(c=10R)。FEM-DEM耦合模型于t=0.1 s時(shí)的模擬結(jié)果如圖4所示。如不采用粗粒化進(jìn)行微-宏觀傳遞(圖4(a))，即每個(gè)在耦合面上的接觸力按狄拉克函數(shù)處理時(shí)，彈性塊沖擊進(jìn)入顆粒床后不能達(dá)到力平衡狀態(tài)，顆粒與單元間的耦合力不足以抵抗彈性塊向下運(yùn)動(dòng)的趨勢。采用粗粒化進(jìn)行微-宏觀傳遞時(shí)，由于每一個(gè)顆粒-單元接觸力均在一定耦合半徑內(nèi)施加到了相鄰而非單個(gè)單元上，彈性體與顆粒床的相互作用在一定時(shí)間后已趨于穩(wěn)定。如圖4(b)所示，顆粒床和彈性塊速度遠(yuǎn)小于圖4(a)中不使用粗?；瘯r(shí)的速度。

圖4 FEM-DEM表面耦合模擬在t =0.1 s的狀態(tài)(|Ve|和|Vp|分別為8節(jié)點(diǎn)有限單元和離散元球形顆粒的速度絕對(duì)值)

3.1.1 耦合系統(tǒng)總能量

圖5(a)給出了FEM-DEM耦合系統(tǒng)z方向總動(dòng)量隨時(shí)間的變化。與基于CG的表面耦合方法相比，常規(guī)耦合方法不能有效降低彈性塊的動(dòng)量。尤其是在彈性塊沖擊顆粒床之后的回彈階段，正方向的動(dòng)量最大絕對(duì)值與負(fù)方向基本一致。采用CG表面耦合方法，系統(tǒng)的總動(dòng)量逐漸減小至零，彈性塊-顆粒耦合系統(tǒng)在t=0.1 s后快速進(jìn)入準(zhǔn)靜態(tài)。

圖5(b)展示了粗?；瘜挾确謩e為0和10R的耦合系統(tǒng)總能量隨時(shí)間的變化?？梢钥闯觯皇褂么至；僮鲿r(shí)，耦合系統(tǒng)總能量從彈性塊接觸顆粒床t=0.04 s時(shí)起不斷增加，在回彈結(jié)束后(t= 0.07s)略微減少，并在二次下落開始時(shí)(t= 0.08 s)再次增加。由于多余的能量多以應(yīng)變能形式進(jìn)入FEM子模型，故有限元計(jì)算在積累誤差到達(dá)一定程度后不再收斂。

圖5 采用常規(guī)耦合和CG耦合方法時(shí)FEM-DEM系統(tǒng)z方向總動(dòng)量和總能量隨時(shí)間的變化

使用粗?；椒詈蠒r(shí)，系統(tǒng)總能量在彈性塊接觸顆粒床之前和耦合系統(tǒng)達(dá)到力平衡狀態(tài)之后均不發(fā)生變化。由顆粒間摩擦產(chǎn)生的能量耗散僅發(fā)生在彈性體與顆粒材料間不斷調(diào)整微觀結(jié)構(gòu)至完成動(dòng)力松弛(t∈[0.04,0.14] s)這一階段。

3.1.2 粗?；瘜挾葘?duì)耦合系統(tǒng)能量的影響

為探究粗?；瘜挾葘?duì)耦合結(jié)果的影響，本文選取粗粒化寬度分別為5R，10R和20R。三種情況下，耦合模擬的可視化和彈性塊沖擊顆粒床后達(dá)到力平衡狀態(tài)的動(dòng)量和能量變化與圖4和圖5類似。限于篇幅，本文不展開討論。

粗?；瘜挾葘?duì)于耦合系統(tǒng)能量的影響主要體現(xiàn)在應(yīng)變能和動(dòng)能上。如圖6所示，耦合系統(tǒng)應(yīng)變能和動(dòng)能隨時(shí)間變化的每個(gè)峰值，隨粗粒化寬度增加而逐漸減小。由此可推斷，基于粗?；谋砻骜詈戏椒捎行У靥幚韽椥詨K沖擊顆粒床所產(chǎn)生的空間、時(shí)間尺度不均勻性，減小計(jì)算累積誤差，增加耦合計(jì)算穩(wěn)定性。

圖6 粗粒化寬度對(duì)耦合系統(tǒng)內(nèi)應(yīng)變能和動(dòng)能的影響

3.2 彈性波在離散-連續(xù)介質(zhì)中的傳播

采用FEM-DEM耦合模型對(duì)懸臂梁結(jié)構(gòu)的動(dòng)力特性進(jìn)行數(shù)值分析，如圖7所示。梁的右端FEM側(cè)節(jié)點(diǎn)固定，左端DEM側(cè)的顆粒無約束。FEM和DEM子模型均不受重力。顆粒之間拉伸與壓縮的彈簧系數(shù)相同；按表1選取的彈簧系數(shù)和幾何尺寸可以得到與FEM模型一致的材料參數(shù)，如彈性模量和泊松比等。

圖7 采用體積耦合方法模擬的懸臂梁(耦合深度DV S=10R)

在t=0時(shí)，最左側(cè)一層的顆粒受到初始速度v0=(0,-0.1,0) m/s的瞬時(shí)沖擊，隨后引起的彈性波將從顆粒系統(tǒng)傳遞到連續(xù)體系統(tǒng)，并在一定時(shí)間后反射回到顆粒系統(tǒng)。為了簡化分析，每一個(gè)體積耦合單元內(nèi)只有一個(gè)顆粒與之耦合。與文獻(xiàn)[18]不同，本文同時(shí)調(diào)整兩個(gè)耦合參數(shù)，即體積耦合深度(coupling depth)DV S=6R～14R和粗?；瘜挾?coarse-graining width)c=0～2.5ΔX。對(duì)于尺寸相同的DEM和FEM子模型來說，耦合計(jì)算域越大，耦合模型所能描述的整體計(jì)算域越小，計(jì)算效率也就越低。本節(jié)主要分析耦合系統(tǒng)的總能量隨耦合深度的變化，而粗?；瘜挾茸鳛榫鶆蚧奈ㄒ粎?shù)，其取值對(duì)降低耦合系統(tǒng)能量耗散也是非常重要的。

3.2.1 體積耦合系統(tǒng)內(nèi)彈性波傳播的時(shí)空演變

由于方向y的波傳播可以考慮為一維問題，本文取每一個(gè)x-z平面內(nèi)顆?；蚬?jié)點(diǎn)速度的平均值。以DV S=10R為例，圖8給出了各個(gè)x-z平面節(jié)點(diǎn)/顆粒平均速度在離散顆粒系統(tǒng)和連續(xù)體系統(tǒng)內(nèi)的時(shí)空演變。可以看出，彈性波由DEM子模型一端開始傳入，至體積耦合區(qū)間后，以相同的波速進(jìn)入FEM子模型。兩子模型內(nèi)波傳播的波速相同且耦合區(qū)域內(nèi)數(shù)值波動(dòng)較小。由此可驗(yàn)證粗?；w積耦合算法的準(zhǔn)確性。

圖8 顆粒系統(tǒng)和連續(xù)體系統(tǒng)內(nèi)彈性波時(shí)空演化

3.2.2 耦合參數(shù)對(duì)體積耦合系統(tǒng)能量耗散的影響

圖9給出了耦合系統(tǒng)(DV S=6R)總能量隨時(shí)間的變化。當(dāng)彈性波進(jìn)入耦合域ΩC,t=0.035 s時(shí)，耦合系統(tǒng)總能量開始耗散，并在波峰完全離開耦合域t=0.055 s后結(jié)束?？梢钥闯?，隨著粗?；瘜挾鹊脑黾樱詈舷到y(tǒng)的能量耗散與不使用粗?；闆r相比明顯減小。針對(duì)本節(jié)數(shù)值算例，取c=2.0ΔX或者c =2.5ΔX時(shí)，可使耦合系統(tǒng)總能量耗散降到最低。

圖9 不同F(xiàn)EM-DEM體積耦合系統(tǒng)(耦合深度DV S=6R，粗?；瘜挾萩 =0～2.5ΔX)的總能量隨時(shí)間的變化

為了更準(zhǔn)確地量化耦合系統(tǒng)的能量耗散，本文定義總能量耗散系數(shù)ΔEt/Et。從圖10可以看出，ΔEt/Et隨著耦合深度和粗?；瘜挾鹊脑黾又饾u減小。在粗?；瘜挾仍黾拥阶顑?yōu)值時(shí)，可以得到相對(duì)較小的能量減衰系數(shù)。耦合深度和粗粒化寬度均可以調(diào)整耦合域動(dòng)態(tài)約束的均勻性，提高罰函數(shù)的約束效應(yīng)。因此，增加耦合深度或者粗?；瘜挾染梢缘玫礁玫鸟詈闲Ч?。但增加耦合深度意味著在耦合范圍內(nèi)需要考慮更多的顆粒和單元，這會(huì)極大地降低計(jì)算效率。因此，更優(yōu)的做法是，在ΔEt/Et相同的前提下，選擇合理的粗粒化寬度，降低耦合顆粒和單元的數(shù)量。

圖10 FEM-DEM體積耦合系統(tǒng)的總能量耗散系數(shù)ΔEt/Et隨粗?；瘜挾群腕w積耦合深度的變化

4 結(jié) 論

本文介紹了如何使用粗粒化方法建立顆粒材料及顆粒-可變形構(gòu)件相互作用的多尺度模型，基于微-宏觀傳遞法則的推導(dǎo)，給出了基于粗?；腇EM-DEM表面耦合和體積耦合公式，解釋了為何常規(guī)局部耦合是粗?；詈系囊粋€(gè)極限情況，通過兩個(gè)數(shù)值算例展示了基于粗?；鸟詈戏椒ň哂懈玫姆€(wěn)定性和更高的效率。主要結(jié)論如下,基于粗?；腇EM-DEM耦合可以通過調(diào)整核函數(shù)影響范圍，使均勻化后的DEM宏觀場與FEM網(wǎng)格相適應(yīng)；非零的粗?；瘜挾仁咕鶆蚧哂蟹蔷植康奶匦裕磫卧c顆粒之間不需要直接接觸或重疊；較優(yōu)的粗?；瘜挾瓤梢越档蛿?shù)值誤差，使耦合模擬效率更高，尤其是在處理多個(gè)顆粒與有限單元耦合時(shí)可以避免使用較大的顆粒-單元界面彈簧，即很小的時(shí)間步長；體積耦合一般需要耦合域范圍足夠大以確保耦合區(qū)間較強(qiáng)的動(dòng)態(tài)約束。合理地增加粗?；瘜挾瓤梢栽诒３窒嗤偰芰亢纳⒌臈l件下，減小耦合域范圍，提高計(jì)算效率。

綜上所述，粗?；菍?shí)現(xiàn)顆粒材料表面與體積多尺度模擬非常方便的數(shù)學(xué)工具。該方法在體積耦合框架下可描述顆粒材料固態(tài)-流態(tài)之間的相互轉(zhuǎn)換或者材料發(fā)生的復(fù)雜物理變化，如3D打印等。在表面耦合框架下可解決諸多顆粒流與可變形結(jié)構(gòu)相互作用的工程問題，如單樁安裝、油井出砂、堵塞與巖石切割。