關(guān)于采用粗?；岣哳w粒材料多尺度模擬守恒特性的研究

程宏旸， Thomas Weinhart

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

1 引 言

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

顆粒材料微觀、細觀結(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)象的宏觀模擬還是非常困難，目前主要依靠顆粒及細觀尺度模型來進行研究。

微-宏觀耦合模擬與單純的微觀或宏觀模擬相比具有明顯優(yōu)勢。通過建立合理的耦合模型，可以保持宏觀方法較高的計算效率和微觀方法較高的計算精度。本文主要探討同步多尺度方法(concurrent multi-scale methods),即微觀與宏觀子模型在同一個計算域上求解并在時間上同步。這需要與基于微-宏觀關(guān)系在物質(zhì)點和表征單元之間傳遞的分層多尺度[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)格的運動信息，并且將顆粒與該界面間的相互作用力映射到FEM模型的牽引邊界上。由于表面耦合通常通過狄拉克函數(shù)來實現(xiàn)[12]，施加到FEM域上的耦合力場很可能是非光滑的，導(dǎo)致難以處理的累積誤差。目前FEM-DEM界面耦合方法的應(yīng)用，如輪胎、土工膜、土工格柵與土顆粒之間的相互作用[12,19,20]，均采用基于狄拉克函數(shù)的微-宏觀映射。

體積耦合又稱為Arlequin多模型耦合框架，由Dhia等[13,14]提出，最初用于耦合具有不同網(wǎng)格分辨率的FEM模型，由Wellmann等[15]首次擴展到FEM-DEM耦合方法上。體積耦合在FEM-DEM模擬上的應(yīng)用主要通過在體積重疊區(qū)域ΩC=ΩF E∩ΩD E(圖1)上施加動態(tài)約束，使微觀和宏觀兩個子模型的解盡可能一致[15]。而如何從微觀顆粒運動信息中提取適應(yīng)于FEM求解空間的宏觀場則是需要解決的難點之一。Wellmann等[15]使用有限元形狀函數(shù)構(gòu)建微-宏觀傳遞法則，將FEM-DEM體積耦合方法應(yīng)用在單樁安裝問題上,即在樁體周圍變形較大區(qū)域使用DEM，遠離樁土作用范圍使用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ù)，為均勻化后的宏觀場定義了一個可調(diào)整的空間尺度?；谟纱至；茖?dǎo)出的表面耦合和體積耦合項公式，將探討如何通過調(diào)整CG寬度減小數(shù)值誤差；以彈性塊沖擊顆粒床和離散-連續(xù)介質(zhì)間的波傳播兩個數(shù)值算例，展示粗粒化寬度對同步多尺度耦合系統(tǒng)的能量守恒特性的影響，并討論如何在保持相同能量耗散的條件下獲得更高的計算效率。

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

在同步多尺度耦合框架內(nèi)，計算域Ω劃分為多個子域；每個子域中的材料行為由不同子模型描述。根據(jù)特定的耦合方法，耦合區(qū)域的定義也有所不同，如ΩC=ΩF E∩ΩD E適用于體積耦合，而表面耦合區(qū)間則定義為ΓC=?ΩF E∩?ΩD E，如圖1所示。本節(jié)將介紹如何使用CG定義空間上連續(xù)的速度場和表面牽引力場，以及如何使用這些均勻化的宏觀場進行表面和體積耦合。上標(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)證明文獻[12,17,18]廣泛使用的表面耦合力公式為粗?；砻骜詈系囊环N特殊情況。

2.1.2 顆粒運動與均勻化后的宏觀速度場

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

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

(7)

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

(8)

(9)

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

2.1.3 基于粗?；奈?宏觀傳遞的一般性

2.2 控制方程與耦合項

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

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

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

(10)

(11)

3 數(shù)值算例

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

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

3.1 彈性體沖擊顆粒床

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

圖2 顆粒受重力在彈性懸臂梁上的運動軌跡

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

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

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

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

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

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

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

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

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

3.1.2 粗?；瘜挾葘︸詈舷到y(tǒng)能量的影響

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

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

圖6 粗?；瘜挾葘︸詈舷到y(tǒng)內(nèi)應(yīng)變能和動能的影響

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

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

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

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

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

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

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

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

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

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

為了更準(zhǔn)確地量化耦合系統(tǒng)的能量耗散，本文定義總能量耗散系數(shù)ΔEt/Et。從圖10可以看出，ΔEt/Et隨著耦合深度和粗?；瘜挾鹊脑黾又饾u減小。在粗粒化寬度增加到最優(yōu)值時，可以得到相對較小的能量減衰系數(shù)。耦合深度和粗?；瘜挾染梢哉{(diào)整耦合域動態(tài)約束的均勻性，提高罰函數(shù)的約束效應(yīng)。因此，增加耦合深度或者粗?；瘜挾染梢缘玫礁玫鸟詈闲Ч５黾玉詈仙疃纫馕吨隈詈戏秶鷥?nèi)需要考慮更多的顆粒和單元，這會極大地降低計算效率。因此，更優(yōu)的做法是，在ΔEt/Et相同的前提下，選擇合理的粗?；瘜挾龋档婉詈项w粒和單元的數(shù)量。

圖10 FEM-DEM體積耦合系統(tǒng)的總能量耗散系數(shù)ΔEt/Et隨粗粒化寬度和體積耦合深度的變化

4 結(jié) 論

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

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