王嗣強(qiáng) 季順迎

(大連理工大學(xué),工業(yè)裝備結(jié)構(gòu)分析國(guó)家重點(diǎn)實(shí)驗(yàn)室,大連 116023)

1 引 言

顆粒物質(zhì)是一種復(fù)雜的能量耗散體系,其摩擦和黏滯作用可使材料所受沖擊載荷的能量有效衰減[1?4].在外載荷的作用下,顆粒間的重新排列調(diào)整接觸力的傳輸方向和接觸時(shí)間,將局部的沖擊載荷在空間擴(kuò)展和時(shí)間延長(zhǎng),進(jìn)而降低沖擊強(qiáng)度[5,6].同時(shí),顆粒間的非彈性碰撞引起的飛濺現(xiàn)象,可將系統(tǒng)吸收的沖擊能量轉(zhuǎn)化為顆粒自身的機(jī)械能,從而達(dá)到耗能緩沖的效果,并已得到相應(yīng)實(shí)驗(yàn)的驗(yàn)證[7,8].此外,顆粒物質(zhì)受沖擊后通常會(huì)產(chǎn)生塑性變形,局部化的應(yīng)變場(chǎng)會(huì)導(dǎo)致固體狀態(tài)和流體狀態(tài)同時(shí)存在并相互轉(zhuǎn)化,是顆粒介質(zhì)類固-液相變的具體體現(xiàn)[9].因此,開(kāi)展非球形顆粒物質(zhì)緩沖耗能的研究有助于揭示顆粒材料的內(nèi)在作用機(jī)理及其主要影響因素,并為緩沖減振技術(shù)的提升和工程應(yīng)用提供很好的借鑒.

顆粒材料的緩沖性能與其幾何形態(tài)有密切關(guān)系,而在自然條件或工業(yè)生產(chǎn)中均是由非球形顆粒構(gòu)成的顆粒系統(tǒng),其在單元排列、動(dòng)力過(guò)程和運(yùn)動(dòng)形態(tài)等方面與球形顆粒有很大的差異[10,11].大量的實(shí)驗(yàn)和數(shù)值模擬主要集中在球形顆粒上,而對(duì)非球形顆粒系統(tǒng)的宏細(xì)觀特性研究則相對(duì)較少[12,13].近年來(lái),顆粒形狀對(duì)系統(tǒng)宏觀力學(xué)性能和動(dòng)態(tài)響應(yīng)過(guò)程的影響引起廣泛關(guān)注[14,15].為合理地描述具有復(fù)雜幾何形態(tài)的顆粒材料,黏接顆粒單元[16?18]、擴(kuò)展多面體單元[19?21]、超二次曲面[22?24]等不同的非球形單元及其構(gòu)造方法不斷發(fā)展和完善起來(lái).其中,超二次曲面單元可構(gòu)造凸面體和凹面體的幾何形狀,并容易地改變函數(shù)參數(shù)進(jìn)而控制顆粒的長(zhǎng)寬比和表面尖銳度.然而,目前超二次曲面單元的接觸力都采用線性近似計(jì)算模型,接觸剛度也通常采用經(jīng)驗(yàn)值,在很大程度上影響計(jì)算精度.基于球體的赫茲接觸模型[25]和基于橢球體表面曲率的接觸模型[26]都為超二次曲面單元非線性接觸力的計(jì)算提供了很好的研究思路.本文在此基礎(chǔ)之上考慮超二次曲面單元在不同接觸模式下接觸點(diǎn)處的曲率半徑,進(jìn)而建立合理的非線性接觸模型.

目前,沖擊載荷作用下顆粒材料的緩沖耗能研究主要通過(guò)對(duì)沖擊坑的深度、沖擊持續(xù)時(shí)間和沖擊物及底板的受力等討論不同沖擊物的速度、角度、質(zhì)量和粒徑以及顆粒床的填充厚度、密集度和顆粒參數(shù)如密度、粒徑、彈性模量、恢復(fù)系數(shù)、摩擦系數(shù)等因素的影響[27?30].此外,實(shí)驗(yàn)研究[31]和離散元計(jì)算[32]表明:顆粒厚度是影響球形顆粒材料緩沖性能的重要因素,并存在一個(gè)臨界厚度;不同表面摩擦系數(shù)及顆粒形態(tài)間緩沖性能存在明顯差異.采用離散單元法(DEM)有助于對(duì)沖擊過(guò)程中顆粒間的接觸作用[33]、沖擊力的傳遞途徑及能量耗散[34]進(jìn)行細(xì)觀分析,從而揭示非球形顆粒材料緩沖耗能的內(nèi)在機(jī)理.

本文基于超二次曲面方程構(gòu)造非球形顆粒形態(tài),采用離散單元法(DEM)建立不同顆粒形狀填充的顆粒床在球形沖擊物作用下的數(shù)值模型,并通過(guò)圓柱體沖擊的理論結(jié)果和球體沖擊的實(shí)驗(yàn)結(jié)果對(duì)計(jì)算模型進(jìn)行驗(yàn)證.在不同顆粒層厚度、顆粒長(zhǎng)寬比及表面尖銳度下,通過(guò)對(duì)底板上沖擊力的數(shù)值計(jì)算,探討顆粒形狀對(duì)顆粒材料緩沖特性的影響規(guī)律,并通過(guò)顆粒材料初始排列的密集度揭示非球形顆粒緩沖性能的內(nèi)在機(jī)理.

2 非球形顆粒的離散元模型

針對(duì)非球形顆粒單元的幾何特點(diǎn),這里采用基于連續(xù)函數(shù)包絡(luò)的超二次曲面方程描述顆粒的復(fù)雜形態(tài),并發(fā)展相應(yīng)的非線性接觸模型計(jì)算單元間的接觸作用,對(duì)沖擊載荷下沖擊物對(duì)底板的作用力進(jìn)行數(shù)值模擬.

2.1 基于連續(xù)函數(shù)包絡(luò)的超二次曲面單元

基于二次曲面方程擴(kuò)展得到的超二次曲面模型是描述非球形顆粒的一種普遍方法,其函數(shù)形式可定義為[35]

式中,a,b和c分別表示顆粒沿主軸方向的半軸長(zhǎng);n1和n2表示顆粒表面的尖銳度參數(shù).對(duì)于n1=n2=2得到球體或橢球體顆粒,n1?2且n2=2得到柱體顆粒,n12且n22得到塊體顆粒.對(duì)于n1,n2→∝時(shí),從理論上可以描述具有尖角的顆粒形態(tài),但這種模型往往受搜索算法的限制而難以實(shí)現(xiàn).目前,通過(guò)超二次曲面方程構(gòu)造非球形顆粒的普遍方法是給定參數(shù)n1和n2的合理取值范圍從而討論顆粒形狀(表面尖銳度)對(duì)宏細(xì)觀特性的影響[15,36,37].因此,本文中n1,n2∈[2,8],且尚需對(duì)這種近似模型與具有尖角的真實(shí)模型進(jìn)一步對(duì)比和分析.圖1顯示在超二次曲面方程中改變不同的半軸長(zhǎng)和表面尖銳度參數(shù)得到不同的單元模型.

圖1 三維超二次曲面單元模型Fig.1.Super-quadric elements in three dimensions.

2.2 單元間的接觸判斷方法

與球形顆粒簡(jiǎn)單且高效的接觸判斷相比,非規(guī)則顆粒間的搜索時(shí)間占整個(gè)計(jì)算時(shí)間的比率會(huì)顯著增加,其計(jì)算效率往往受顆粒形狀、邊界條件、搜索算法等因素的影響而具有顯著差異.本文針對(duì)非球形顆粒間接觸判斷的復(fù)雜性,考慮顆粒間不同接觸模式下的顆粒取向,采用相應(yīng)的球形包圍盒和方向包圍盒(oriented bounding box,OBB),如圖2所示.這兩種方法可有效減少潛在的接觸對(duì),提高程序的計(jì)算效率.球形包圍盒作為單元間第一次接觸的粗判斷,其接觸理論基于球體與球體間的接觸計(jì)算,包圍球體半徑可表示為如果兩球體間球心的距離小于其半徑和,則接觸,否則不接觸.OBB作為單元間第二次接觸的粗判斷,其接觸理論基于分離軸算法并考慮顆粒取向和長(zhǎng)寬比[38].如果兩個(gè)幾何體在空間向量上的投影發(fā)生重疊,則接觸,否則不接觸.

圖2 球形包圍盒和OBBFig.2.Bounding spheres and oriented bounding boxes.

非線性牛頓迭代方法作為單元間第三次接觸判斷,該方法將求解單元間最短距離的優(yōu)化問(wèn)題轉(zhuǎn)化為非線性方程組進(jìn)行求解,如圖3所示.考慮非球形顆粒接觸時(shí)表面法向平行且反向,同時(shí)基于中間點(diǎn)到兩個(gè)顆粒表面的幾何勢(shì)能相等建立非線性方程組[39,40]:

式中,X=(x,y,z)T,Fi(X)和Fj(X)分別表示全局坐標(biāo)下顆粒i和j的函數(shù)方程,?F(X)表示顆粒表面沿x,y,z方向的外法向,λ2表示兩個(gè)顆粒間的外法向平行且反向.

圖3 超二次曲面單元間的接觸檢測(cè)Fig.3.Contact detection between super-quadric particles.

由此,牛頓迭代公式可寫作:

式中,X=X+dX,λ=λ+dλ.如果計(jì)算結(jié)果X0滿足Fi(X0)<0且Fj(X0)<0,則表明兩個(gè)單元間發(fā)生接觸,接觸法向定義為n=?Fi(X0)/∥?Fi(X0)∥.進(jìn)一步考慮顆粒表面點(diǎn)Xi和Xj滿足:Xi=X0+αn和Xj=X0+βn,且未知參數(shù)α和β可通過(guò)一元非線性牛頓迭代得到:和βk+1=因此,法向重疊量可表述為δn=Xi?Xj.

與單元間的接觸搜索類似,單元與固體邊界的接觸判斷通常將邊界表示為函數(shù)形式,通過(guò)聯(lián)立超二次曲面方程進(jìn)而建立非線性方程組進(jìn)行求解.本文中主要固體邊界是剛性圓筒,其函數(shù)形式可表示為:式中,D0為圓筒直徑,H0為圓筒高度,如圖4所示.考慮超二次曲面單元i與圓筒接觸時(shí)表面法向平行且同向,同時(shí)基于中間點(diǎn)到顆粒表面和圓筒表面的幾何勢(shì)能相等建立非線性方程組:

圖4 超二次曲面單元與圓筒邊界的接觸檢測(cè)Fig.4.Contact detection between a super-quadric particle and a cylinder.

由此,相應(yīng)的牛頓迭代公式可表示為

式中,X=X+dX,λ=λ+dλ.如果計(jì)算結(jié)果滿足則表明超二次曲面單元與圓筒邊界發(fā)生接觸,接觸法向定義為這里,圓筒表面點(diǎn)可表示為:考慮顆粒表面點(diǎn)且未知量α通過(guò)一元非線性牛頓迭代公式得:因此,單元與邊界的法向重疊量可表述為

一般而言,求解四元非線性方程組的牛頓迭代算法是影響單元搜索效率及整體計(jì)算時(shí)間的主要因素,且與顆粒表面尖銳度和長(zhǎng)寬比密切相關(guān).隨著尖銳度參數(shù)n1和n2及長(zhǎng)寬比偏離1的指數(shù)增加,單元間的搜索效率顯著降低;當(dāng)顆粒接近于球形時(shí),其計(jì)算效率最高[23].

2.3 超二次曲面單元間的非線性接觸模型

針對(duì)超二次曲面單元間不同的接觸模式,在傳統(tǒng)球形非線性接觸模型[25]的基礎(chǔ)上,考慮單元間局部接觸點(diǎn)處的平均曲率半徑[41],進(jìn)而建立相對(duì)合理的非線性接觸模型.因此,單元間的接觸力主要是由彈性力和黏滯力組成,同時(shí)考慮法向力、切向力和滾動(dòng)摩擦引起的附加力矩.其中,法向作用力Fn可寫作:

式中,Kn為兩個(gè)接觸單元的法向有效剛度;δn為法向重疊量;Cn為法向阻尼系數(shù);A為顆粒間的黏滯參數(shù);vn,ij為法向相對(duì)速度.這里法向有效剛度Kn、黏滯參數(shù)A分別為

式中,Kmean為接觸點(diǎn)處的局部平均曲率,其可寫作[41]:

其中,?F,?2F分別表示超二次曲面方程的一階和二階導(dǎo)函數(shù).

單元間的切向接觸力Fs主要由彈性力和黏滯力組成,同時(shí)考慮Mohr-Coulomb摩擦定律.則切向彈性力為

式中,μs為滑動(dòng)摩擦系數(shù),δt為切向重疊量,δt,max=μs(2?υ)/2(1?υ)·δn.

這里,切向黏滯力為

式中,Ct為切向黏滯系數(shù),vt,ij為單元間切向相對(duì)速度.

在單元發(fā)生相對(duì)轉(zhuǎn)動(dòng)時(shí),由滾動(dòng)摩擦引起的力矩Mr可表述為[42]

式中,μr為滾動(dòng)摩擦系數(shù),即本文中為法向的單位相對(duì)轉(zhuǎn)速,可表示為

目前,關(guān)于非球形顆粒時(shí)間步長(zhǎng)的理論研究相對(duì)較少,通常是在球形顆粒理論計(jì)算的基礎(chǔ)上給出合理的時(shí)間步長(zhǎng)[15,20,26].因此,本文中時(shí)間步長(zhǎng)的選擇是考慮非線性接觸模型的球形顆粒單元[43],即

式中,Rmin為顆粒的最小半徑,ρ為顆粒密度,G為剪切模量,υ為泊松比.為了保證計(jì)算精度和運(yùn)行效率,Rmin為超二次曲面單元三個(gè)方向半軸長(zhǎng)的最小值,即Rmin=min(a,b,c),且選擇dt6 0.1tmax作為計(jì)算步長(zhǎng)的依據(jù).

圖5 考慮等效曲率的非線性接觸模型Fig.5.Non-linear force model with considering the equivalent radius of the particle shape curvature.

3 顆粒材料沖擊過(guò)程的理論與實(shí)驗(yàn)驗(yàn)證

為驗(yàn)證超二次曲面單元間算法和接觸模型的有效性,分別對(duì)單個(gè)圓柱體沖擊的理論結(jié)果和球體沖擊的實(shí)驗(yàn)結(jié)果與離散元計(jì)算值進(jìn)行對(duì)比驗(yàn)證,從而更加合理地模擬沖擊力的變化規(guī)律.

3.1 圓柱體沖擊過(guò)程的理論驗(yàn)證

通過(guò)超二次曲面模型構(gòu)造一個(gè)理想圓柱體顆粒,給定初始沖擊速度下考慮回彈速度和回彈角速度隨沖擊角度的變化關(guān)系,同時(shí)假定整個(gè)過(guò)程無(wú)摩擦和重力作用,如圖6所示.沖擊后的角速度和速度可表述為[44]

式中,m為圓柱體的質(zhì)量;為沖擊前圓柱體形心的速度;ε為接觸點(diǎn)處的回彈系數(shù);φ為接觸點(diǎn)與顆粒形心的連線與圓柱體表面的夾角;r為接觸點(diǎn)到顆粒形心的距離;θ為顆粒表面與邊界的夾角,即沖擊角度;Iyy為顆粒繞y軸的轉(zhuǎn)動(dòng)慣量.在不同角度的沖擊過(guò)程中,初始沖擊速度恒定為1 m/s,回彈系數(shù)ε恒定為0.85.當(dāng)顆粒與邊界接觸碰撞后,立刻移除邊界從而避免在沖擊角度過(guò)小或過(guò)大時(shí)產(chǎn)生的二次碰撞.離散元模擬的主要參數(shù)列于表1中.

圖6 圓柱體沖擊的示意圖,x-z平面的投影Fig.6.Scheme of cylinder-wall impact,x-z projection.

圖7(a)和圖7(b)顯示了0?—90?范圍內(nèi)不同沖擊角度下無(wú)量綱的回彈速度和角速度與理論結(jié)果[44]的對(duì)比.值得注意的是,當(dāng)沖擊角度分別為0?和90?時(shí),圓柱體與平面的接觸不再為點(diǎn)接觸.因此,回彈速度的理論解僅與回彈系數(shù)相關(guān),而回彈角速度的理論解為零.從圖7可以看出,無(wú)量綱的回彈速度約為0.85,無(wú)量綱的回彈角速度約為0.離散元的計(jì)算結(jié)果與理論結(jié)果[44]基本符合,表明超二次曲面模型和算法適用于模擬非球形顆粒的沖擊過(guò)程.

表1 圓柱體沖擊離散元模擬的主要計(jì)算參數(shù)Table 1.DEM simulation parameters of cylinder impact.

圖7 圓柱體沖擊的解析結(jié)果[44]與離散元計(jì)算值對(duì)比(a)無(wú)量綱的回彈速度(b)無(wú)量綱的回彈角速度Fig.7.Comparison of cylinder-wall impact between analytical results[44]and DEM simulation results:(a)The dimensionless post-impact translational velocity(b)the dimensionless post-impact angular velocity

3.2 顆粒材料沖擊實(shí)驗(yàn)的對(duì)比驗(yàn)證

為驗(yàn)證超二次曲面離散元計(jì)算模型的可靠性,分別對(duì)顆粒層厚度為0.5和2.0 cm的球體沖擊過(guò)程進(jìn)行數(shù)值模擬,并與文獻(xiàn)[31]的實(shí)驗(yàn)結(jié)果進(jìn)行對(duì)比.離散元模擬中球形顆粒單元的粒徑按均勻概率密度函數(shù)在[4.0 mm,5.0 mm]內(nèi)隨機(jī)分布,且均值為4.5 mm,主要的離散元計(jì)算參數(shù)列于表2中.

當(dāng)顆粒層厚度分別為H=0.5和2.0 cm時(shí),圓筒底部受沖擊力的時(shí)程如圖8(a)和圖8(b)所示.當(dāng)顆粒層厚度H=0.5 cm時(shí),沖擊力峰值與沖擊持續(xù)時(shí)間與實(shí)驗(yàn)結(jié)果基本符合,但當(dāng)顆粒層厚度H=2.0 cm時(shí),沖擊力峰值低于實(shí)驗(yàn)結(jié)果,沖擊持續(xù)時(shí)間則大于實(shí)驗(yàn)值.同時(shí),實(shí)驗(yàn)中沖擊力呈現(xiàn)多次波動(dòng)現(xiàn)象,而離散元結(jié)果較為穩(wěn)定.這主要是由于實(shí)驗(yàn)測(cè)量的沖擊過(guò)程中,沖擊物引起圓筒底板產(chǎn)生一定的振動(dòng),進(jìn)而出現(xiàn)明顯的波動(dòng)現(xiàn)象.但在離散元模擬中底板通常設(shè)成剛性板,因此沖擊力的振蕩幅度不明顯.盡管離散元計(jì)算的數(shù)值結(jié)果與實(shí)驗(yàn)結(jié)果存在一定的偏差,但沖擊力的變化規(guī)律是一致的,進(jìn)一步表明基于超二次曲面模型的離散元方法可以合理地模擬顆粒系統(tǒng)的沖擊緩沖過(guò)程.

表2 球體沖擊離散元模擬的主要計(jì)算參數(shù)Table 2.DEM simulation parameters of sphere impact.

圖8 不同顆粒厚度下底板受力變化過(guò)程的實(shí)驗(yàn)值[31]和離散元計(jì)算值對(duì)比 (a)顆粒厚度為0.5 cm;(b)顆粒厚度為2.0 cmFig.8.Comparison of impact loads between experiment results[31]and DEM simulation results under different granular thicknesses:(a)H=0.5 cm;(b)H=2.0 cm.

4 顆粒層厚度與形狀對(duì)緩沖性能的影響

顆粒層厚度和顆粒形狀是影響緩沖特性的重要因素.對(duì)于圓筒內(nèi)球形顆粒材料的沖擊實(shí)驗(yàn)[31]和離散元結(jié)果[32]表明,底板所受的沖擊力隨球形顆粒層厚度的增加而減小,當(dāng)厚度大于臨界厚度時(shí),沖擊力峰值變化趨于穩(wěn)定.同時(shí),非球形顆粒呈現(xiàn)出低流動(dòng)性、高互鎖性等不同于球形顆粒的特殊性質(zhì)[11].為獲得不同顆粒形狀間較好的沖擊過(guò)程,這里令沖擊物的初始速度V0=5 m/s,在顆粒層表面上的初始高度H2=30 cm,摩擦系數(shù)μs=0.4,法向阻尼系數(shù)Cn=0.1,顆粒密度ρ=2500 kg/m3,顆粒半徑r0=5 mm,其余計(jì)算參數(shù)見(jiàn)表2.

4.1 顆粒層厚度的影響

通過(guò)超二次曲面方程構(gòu)造球體、圓柱體、正方體三種不同的顆粒形狀,在剛性圓筒中隨機(jī)產(chǎn)生,通過(guò)落雨法生成不同厚度的顆粒床,并在重力作用下顆粒系統(tǒng)保持穩(wěn)定平衡狀態(tài),如圖9所示.

當(dāng)無(wú)顆粒填充時(shí),圓筒底部所受的沖擊力峰值P0=4.67 kN,當(dāng)顆粒層厚度H=0.8,1.5,2.5,4.5,9.5 cm時(shí),計(jì)算得到的沖擊力時(shí)程如圖10(a)—(c)所示.從圖10(a)—(c)可以發(fā)現(xiàn),球體、圓柱體、正方體都呈現(xiàn)類似的沖擊特性,即隨著顆粒層厚度的增加,沖擊力峰值P逐漸減小且沖擊持續(xù)時(shí)間逐漸延長(zhǎng).將三種形狀在不同顆粒層厚度下的沖擊力峰值進(jìn)行統(tǒng)計(jì),如圖10(d)所示.可以發(fā)現(xiàn),當(dāng)顆粒層厚度H<7.0 cm時(shí),相比圓柱形和正方形顆粒,球形顆粒的緩沖性能最好,而圓柱形顆粒的緩沖效果優(yōu)于正方形顆粒.當(dāng)顆粒層厚度H>7.0 cm時(shí),顆粒層厚度和顆粒形狀對(duì)緩沖特性的影響不再顯著,且趨于穩(wěn)定.這里稱該厚度為臨界顆粒層厚度Hc,即Hc=7.0 cm.由于不同顆粒形態(tài)的緩沖特性受材料參數(shù)、尺度效應(yīng)、沖擊能量等不同因素的影響,因此在不同條件下臨界顆粒層厚度通常具有一定差異.

圖9 不同形狀顆粒材料緩沖特性的離散元模型 (a)球體;(b)圓柱體;(c)正方體Fig.9.Discrete element models of granular materials with different shapes for bu ff er properties study:(a)Sphere;(b)cylinder;(c)cube.

圖10 不同顆粒層厚度下底板力的沖擊力時(shí)程曲線 (a)球體;(b)圓柱體;(c)正方體;(d)三種形狀的沖擊力峰值變化Fig.10.Impact loads on bottom plate versus time under various thicknesses of granular with different shapes:(a)Sphere;(b)cylinder;(c)cube;(d)comparison of impact load peaks among three particle shapes.

4.2 不同顆粒尖銳度下密集度對(duì)沖擊力的影響

為研究顆粒表面尖銳度對(duì)顆粒材料緩沖性能的影響,這里保證不同形態(tài)的單個(gè)顆粒質(zhì)量相等,從而獲得不同顆粒形狀影響下的緩沖性能.由于超二次曲面方程產(chǎn)生的顆粒形狀具有幾何對(duì)稱性,這里給出1/4的顆粒表面輪廓隨函數(shù)尖銳度參數(shù)從2.0到8.0的變化過(guò)程,如圖11(a)所示.通過(guò)表面尖銳度參數(shù)的連續(xù)變化實(shí)現(xiàn)顆粒形狀從球體到正方體和圓柱體到正方體的兩種形態(tài)轉(zhuǎn)變,如圖11(b)所示.可以發(fā)現(xiàn),從球體到正方體變化的過(guò)程中保持參數(shù)n1=n2,而圓柱體到正方體變化的過(guò)程中保持參數(shù)n1=8.通過(guò)改變參數(shù)n2,從而得到在參數(shù)n2從2到8變化過(guò)程中13種具有不同表面尖銳度的顆粒形態(tài).

圖11 不同顆粒表面尖銳度參數(shù)對(duì)顆粒形狀的影響Fig.11.In fluences of blockiness parameters on the particle shapes.

圖12為從球體和圓柱體到正方體的兩種形狀變化過(guò)程中表面尖銳度參數(shù)分別為2.0,3.0,5.0,8.0時(shí)底板所受的沖擊力時(shí)程.在兩種形狀變化過(guò)程中,對(duì)不同尖銳度下沖擊力峰值p和顆粒系統(tǒng)初始密集度C0進(jìn)行統(tǒng)計(jì),如圖13所示.可以發(fā)現(xiàn),沖擊力峰值和初始密集度隨表面尖銳度的增加而顯著增加.在顆粒隨機(jī)堆積和沖擊過(guò)程中,表面尖銳度的主要作用是將顆粒間的點(diǎn)接觸逐步轉(zhuǎn)變?yōu)槊娼佑|,增加顆粒間的接觸面積,同時(shí)產(chǎn)生更加穩(wěn)定且密實(shí)的顆粒系統(tǒng).此外,增加顆粒表面尖銳度能阻止顆粒間的相對(duì)滑動(dòng)和滾動(dòng),使沖擊物在顆粒材料中的運(yùn)行距離相對(duì)減小,從而產(chǎn)生較高的底部沖擊力.以上結(jié)果表明:相比具有尖銳頂點(diǎn)和平面的顆粒形態(tài),光滑表面的球形顆粒具有更好的緩沖效果.同時(shí),顆粒間的面接觸會(huì)提高顆粒系統(tǒng)的密實(shí)度,從而降低系統(tǒng)的緩沖性能.

4.3 不同顆粒長(zhǎng)寬比下密集度對(duì)沖擊力的影響

為進(jìn)一步研究顆粒長(zhǎng)寬比對(duì)緩沖特性的影響,通過(guò)改變超二次曲面方程得到不同的長(zhǎng)寬比σ=c/a(=b),分別取σ=0.4,0.6,0.8,1.0,1.5,2.0,2.5和3.0得到8種不同長(zhǎng)寬比的圓柱形和正方形顆粒形狀,如圖14所示.

圖12 不同形狀變化過(guò)程中表面尖銳度對(duì)沖擊力時(shí)程的影響 (a)球體-正方體;(b)圓柱體-正方體Fig.12.In fluences of blockiness parameters on the impact load versus time for the different changing processes:(a)Sphere-cube;(b)cylinder-cube.

圖13 不同表面尖銳度下沖擊力峰值和初始密集度的變化 (a)沖擊力峰值;(b)初始密集度Fig.13.Impact load peaks and initial granular concentration under various blockiness parameters:(a)Impact load peaks;(b)initial concentration of granular material.

圖14 不同長(zhǎng)寬比的圓柱形和長(zhǎng)方形顆粒 (a)圓柱形顆粒;(b)長(zhǎng)方形顆粒Fig.14.Cylinder-like particles and cube-like particles with various aspect ratios:(a)Cylinder-like particles;(b)cubelike particles.

圖15為圓柱形和長(zhǎng)方形顆粒的長(zhǎng)寬比σ分別為0.4,1.0,1.5和2.5時(shí)底板的沖擊力時(shí)程曲線.同時(shí),分別對(duì)兩種形狀不同長(zhǎng)寬比的沖擊力峰值和初始密集度進(jìn)行統(tǒng)計(jì),如圖16所示.可以發(fā)現(xiàn),增加或減小長(zhǎng)寬比會(huì)使底部沖擊力峰值和顆粒系統(tǒng)的初始密集度降低,從而提高顆粒材料的緩沖性能.同時(shí),對(duì)于相同長(zhǎng)寬比的顆粒形態(tài),長(zhǎng)方形顆粒的底部沖擊力峰值都高于圓柱形顆粒.這主要是由于長(zhǎng)方形顆粒間的主要接觸模式為面接觸,從而產(chǎn)生更加密實(shí)且相對(duì)穩(wěn)定的顆粒系統(tǒng),縮短了沖擊的持續(xù)時(shí)間,使得長(zhǎng)方形顆粒具有較弱的緩沖性能.

在沖擊過(guò)程中,顆粒長(zhǎng)寬比的主要作用是調(diào)整緊密排列的顆粒系統(tǒng),使顆粒間存在更多孔隙,顆粒在沖擊載荷作用下有很大的自由移動(dòng)空間.同時(shí),顆粒長(zhǎng)寬比降低了系統(tǒng)的穩(wěn)定性,增加顆粒間的相對(duì)滑動(dòng)和轉(zhuǎn)動(dòng),延長(zhǎng)沖擊時(shí)間.此外,對(duì)于較高或較低長(zhǎng)寬比的顆粒系統(tǒng),顆粒間的自鎖可以調(diào)整沖擊力的傳輸方向并實(shí)現(xiàn)分散響應(yīng),從而將局部沖擊力在空間進(jìn)行擴(kuò)展,具有較好的緩沖效果.

圖15 不同圓柱形顆粒和長(zhǎng)方形顆粒長(zhǎng)寬比對(duì)沖擊力時(shí)程的影響 (a)圓柱形顆粒;(b)長(zhǎng)方形顆粒Fig.15.In fluences of aspect ratios on the impact load versus time for the different shapes:(a)Cylinder-like particles;(b)cube-like particles.

圖16 不同長(zhǎng)寬比下沖擊力峰值和初始密集度的變化 (a)沖擊力峰值;(b)初始密集度Fig.16.Impact load peaks and initial granular concentration under various aspect ratios:(a)Impact load peaks;(b)initial concentration of granular material.

5 結(jié) 語(yǔ)

基于超二次曲面方程構(gòu)造非球形顆粒形態(tài),在傳統(tǒng)非線性接觸模型的基礎(chǔ)上考慮局部接觸點(diǎn)處的平均曲率,采用離散元方法對(duì)顆粒材料在球形沖擊物作用下的緩沖特性進(jìn)行了數(shù)值分析.通過(guò)與圓柱體沖擊的理論結(jié)果和球體沖擊的實(shí)驗(yàn)結(jié)果進(jìn)行對(duì)比驗(yàn)證,表明基于超二次曲面模型的離散元方法適用于模擬非球形顆粒的沖擊過(guò)程.在此基礎(chǔ)之上,進(jìn)一步研究了顆粒層厚度、顆粒表面尖銳度和長(zhǎng)寬比對(duì)底部沖擊力的影響.研究結(jié)果表明:顆粒層厚度是影響不同形態(tài)顆粒材料緩沖特性的重要因素.當(dāng)顆粒層厚度HHc時(shí),緩沖性能對(duì)顆粒形狀和顆粒層厚度不再敏感.同時(shí),相比具有尖銳頂點(diǎn)和平面的顆粒形態(tài),表面光滑的球形顆粒具有更好的緩沖效果.顆粒間的面接觸會(huì)提高顆粒系統(tǒng)的初始密集度并增強(qiáng)系統(tǒng)的穩(wěn)定性,從而降低顆粒材料的緩沖性能.此外,增加或減小圓柱形和長(zhǎng)方形顆粒的長(zhǎng)寬比,有利于顆粒間的相對(duì)滑動(dòng)或轉(zhuǎn)動(dòng),同時(shí)對(duì)于較高或較低長(zhǎng)寬比的顆粒形態(tài),顆粒間的自鎖都會(huì)提高系統(tǒng)的緩沖性能.

參考文獻(xiàn)

[1]Katsuragi H,Durian D J 2007Nat.Phys.3 420

[2]Kondic L,Fang X,Losert W,O’Hern C S,Behringer R P 2012Phys.Rev.E85 011305

[3]Nordstrom K,Lim E,Harrington M,Losert W 2014Phys.Rev.Lett.112 228002

[4]Bester C S,Behringer R P 2017Phys.Rev.E95 032906

[5]Seguin A,Bertho Y,Gondret P,Crassous J 2009Europhys.Lett.88 44002

[6]Clark A H,Petersen A J,Kondic L,Behringer R P 2015Phys.Rev.Lett.114 144502

[7]Deboeuf S,Gondret P,Rabaud M 2008Environ.Sci.Technol.42 8459

[8]Deboeuf S,Gondret P,Rabaud M 2009Phys.Rev.E79 041306

[9]Ye X Y,Wang D M,Zheng X J 2012Phys.Rev.E86 061304

[10]Lu G,Third J R,Müller C R 2015Chem.Eng.Sci.127 425

[11]Zhong W,Yu A,Liu X,Tong Z,Zhang H 2016Powder Technol.302 108

[12]Zhu H P,Zhou Z Y,Yang R Y,Yu A B 2007Chem.Eng.Sci.62 3378

[13]Zhu H P,Zhou Z Y,Yang R Y,Yu A B 2008Chem.Eng.Sci.63 5728

[14]Elskamp F,Kruggel-Emden H,Hennig M,Teipel U 2017Granular Matter19 46

[15]Zhao S,Zhang N,Zhou X,Zhang L 2017Powder Technol.310 175

[16]Kruggel-Emden H,Rickelt S,Wirtz S,Scherer V 2008Powder Technol.188 153

[17]Li C Q,Xu W J,Meng Q S 2015Powder Technol.286 478

[18]Zeng Y,Jia F,Zhang Y,Meng X,Han Y,Wang H 2017Powder Technol.313 112

[19]Galindo-Torres S A,Pedroso D M 2010Phys.Rev.E81 061303

[20]Toson P,Khinast J G 2017Powder Technol.313 353

[21]Govender N,Wilke D N,Pizette P,Abriak N E 2018Appl.Math.Comput.319 318

[22]Lu G,Third J R,Müller C R 2012Chem.Eng.Sci.78 226

[23]Cui Z Q,Chen Y C,Zhao Y Z,Hua Z L,Liu X,Zhou C L 2013Chin.J.Computat.Mech.30 854(in Chinese)[崔澤群,陳友川,趙永志,花爭(zhēng)立,劉驍,周池樓 2013計(jì)算力學(xué)學(xué)報(bào)30 854]

[24]Cleary P W,Sinnott M D,Morrison R D,Cummins S,Delaney G W 2017Miner.Eng.100 49

[25]Di Renzo A,Di Maio F P 2004Chem.Eng.Sci.59 525

[26]Liu S D,Zhou Z Y,Zou R P,Pinson D,Yu A B 2014Powder Technol.253 70

[27]Goldman D I,Umbanhowar P 2008Phys.Rev.E77 021308

[28]Vet S J D,Bruyn J R D 2012Granular Matter14 661

[29]Clark A H,Petersen A J,Behringer R P 2014Phys.Rev.E89 012201

[30]Clark A H,Kondic L,Behringer R P 2016Phys.Rev.E93 050901

[31]Ji S Y,Li P F,Chen X D 2012Acta Phys.Sin.61 184703(in Chinese)[季順迎,李鵬飛,陳曉東2012物理學(xué)報(bào)61 184703]

[32]Ji S Y,Fan L F,Liang S M 2016Acta Phys.Sin.65 104501(in Chinese)[季順迎,樊利芳,梁紹敏2016物理學(xué)報(bào)65 104501]

[33]Sun Q C,Wang G Q 2008Acta Phys.Sin.57 4667(in Chinese)[孫其誠(chéng),王光謙 2008物理學(xué)報(bào) 57 4667]

[34]Peng Z,Jiang Y M,Liu R,Hou M Y 2013Acta Phys.Sin.62 024502(in Chinese)[彭政,蔣亦民,劉銳,厚美瑛2013物理學(xué)報(bào)62 024502]

[35]Barr A H 1981IEEE Comput.Graph.Appl.1 11

[36]Stenzel O,Salzer M,Schmidt V,Cleary P W,Delaney G W 2014Granular Matter16 457

[37]Delaney G W,Cleary P W 2010Europhys.Lett.89 34002

[38]Portal R,Dias J,de Sousa L 2010Arch.Mech.Eng.57 165

[39]Wellmann C,Lillie C,Wriggers P 2008Eng.Computat.25 432

[40]Podlozhnyuk A,Pirker S,Kloss C 2017Comp.Part.Mech.4 101

[41]Goldman R 2005Comput.Aided Geomet.Desig.22 632

[42]Gan J Q,Zhou Z Y,Yu A B 2017Powder Technol.320 610

[43]Kremmer M,Favier J F 2001Int.J.Numer.Meth.Eng.51 1407

[44]Kodam M,Bharadwaj R,Curtis J,Hancock B,Wassgren C 2010Chem.Eng.Sci.65 5863