顆粒-顆粒接觸力的熱力學(xué)模型?

蔣亦民 劉佑

1)(中南大學(xué)物理與電子學(xué)院,長(zhǎng)沙 410083)

2)(蒂賓根大學(xué)理論物理研究所,德國(guó),蒂賓根 72076)

1 引 言

顆粒物質(zhì)(諸如沙堆等)與傳統(tǒng)晶體材料(例如半導(dǎo)體硅)的重要不同是多出一個(gè)介觀(或細(xì)觀)顆粒尺度.這使得從電子-原子核之間的二體庫侖靜電力出發(fā),計(jì)算其物理力學(xué)性質(zhì)的所謂“從頭算(ab initio)”思想,雖然原則概念上是存在的,但實(shí)際操作上沒有可能性.也就是說對(duì)硅晶體可以從頭計(jì)算,但對(duì)顆粒物質(zhì)最小尺度的分析計(jì)算只能從細(xì)觀顆粒出發(fā).另外為了計(jì)算數(shù)目龐大的顆粒系統(tǒng),還不能像赫茲那樣將顆粒處理成可變形的彈性固體[1],而是必須簡(jiǎn)化為只有六個(gè)運(yùn)動(dòng)學(xué)(kinematic)變量的、帶轉(zhuǎn)動(dòng)的“質(zhì)點(diǎn)粒子”.如果要求顆粒i的六個(gè)運(yùn)動(dòng)學(xué)量,即質(zhì)心位矢和角位矢{ri(t),φi(t)}遵守牛頓力和力矩方程的話,則需要建立{ri,φi}與其他顆粒j之間的作用力fij關(guān)系模型來閉合牛頓方程組,然后對(duì)其數(shù)值求解,可得到該質(zhì)點(diǎn)系統(tǒng)細(xì)觀結(jié)構(gòu)和力結(jié)構(gòu)的時(shí)間演化(我們將省略與本文內(nèi)容無關(guān)的邊界條件).如果進(jìn)一步選擇合適的物理體積元做這些細(xì)觀信息的統(tǒng)計(jì)平均(諸如粗?；骄?,還可得各種宏觀平均量及其漲落的時(shí)空變化以及各種關(guān)聯(lián)性質(zhì).這套稱作軟球離散元(DEM)模擬方法的一個(gè)優(yōu)勢(shì)是可以描述以長(zhǎng)程關(guān)聯(lián)為特征的彈性力鏈的呈現(xiàn)與消亡現(xiàn)象,又稱固液轉(zhuǎn)變或jamming轉(zhuǎn)變.也就是說軟球DEM適用于顆粒物質(zhì)的所有類固液氣態(tài)及其轉(zhuǎn)變[2].

軟球DEM面臨的顆粒-顆粒接觸作用力建模問題,在晶體物性的“從頭算”領(lǐng)域里是不存在的.“從頭算”認(rèn)為庫侖靜電力能夠解決所有凝聚態(tài)物理問題[3,4],但對(duì)顆粒物理而言接觸力建模才是它的一個(gè)基礎(chǔ)課題.由于這是個(gè)帶耗散的不可逆非保守力,原則上需要在熱力學(xué)概念基礎(chǔ)上建模,具體地講就是除了接觸力外,還應(yīng)該給出機(jī)械能功率與熱功率的比例.或者給出計(jì)算彈性勢(shì)能的表達(dá)式,因?yàn)榻佑|時(shí)的總能量變化功率可以將作用力乘以運(yùn)動(dòng)學(xué)速度得到:fij(vi?vj)(這里vi≡dri/dt),動(dòng)能的計(jì)算公式又是熟知的,故一旦明確了接觸過程彈性勢(shì)能的計(jì)算方式,問題也就徹底解決了.注意著名的赫茲法向接觸力fn～|ri?rj|3/2[1],是忽略耗散時(shí)的結(jié)果.有耗散時(shí),特別是有塑性耗散時(shí)(這在切向運(yùn)動(dòng)時(shí)必須考慮),接觸力是不能簡(jiǎn)單地表達(dá)為顆粒幾何運(yùn)動(dòng)學(xué)變量的代數(shù)函數(shù)型這樣的捆綁關(guān)系的.這個(gè)復(fù)雜性來自接觸力中有彈性貢獻(xiàn),而彈性是需要用額外的獨(dú)立變量來表征的概念[3,4]:即所謂的平移對(duì)稱破缺后呈現(xiàn)的Goldstone量(彈應(yīng)變或彈簧長(zhǎng)度變化).很多模型都注意到了有塑性時(shí)彈簧長(zhǎng)度變化與幾何運(yùn)動(dòng)學(xué)相對(duì)量{ri?rj,φi?φj}的關(guān)系需要用微積分來描述,特別是切向.例如Luding模型在處理切向力時(shí),就引入了一個(gè)切向彈簧長(zhǎng)度變量,它與運(yùn)動(dòng)學(xué)量的聯(lián)系是利用增量關(guān)系模型和一個(gè)歷史累積變量來定義[5].有些人不用彈簧長(zhǎng)度,直接建立切向力與運(yùn)動(dòng)學(xué)量的增量關(guān)系模型(參見文獻(xiàn)[2])(工程領(lǐng)域常常將某些微積分方程寫成增量關(guān)系形式,這個(gè)習(xí)慣很可能源自Truesdell[6]).長(zhǎng)期以來,接觸力建模研究一直處于這樣的唯象力學(xué)階段,并未去澄清熱力學(xué)量,即熱功率或彈性勢(shì)能的定義.值得提到的是,很多(例如Hertz-Mindlin)三維接觸力模型的可逆部分采用了虎克力乘以重疊量平方根的Boussinesq模型[7,8].這里雖然不涉及耗散,由于Boussinesq力不是保守的,我們?nèi)詿o法知道彈性勢(shì)能.

由于這些問題,雖然有少數(shù)DEM工作計(jì)算了彈性勢(shì)能與動(dòng)能比值[9,10],一個(gè)嚴(yán)格熱力學(xué)的接觸力模型仍值得建立.考慮到目前一個(gè)嚴(yán)格熱力學(xué)的顆粒物質(zhì)連續(xù)力學(xué)理論——顆粒固體流體動(dòng)力學(xué)(GSH)已經(jīng)存在[11?14],接觸力的熱力學(xué)模型完全可以用類似的方案得到(第二節(jié)).這時(shí)GSH中處理庫侖屈服、塑性、滯迴等復(fù)雜力學(xué)現(xiàn)象的成功方法,可以直接繼承到接觸力上.例如前面提到的Boussinesq力沒有勢(shì)能的問題,完全可以連同庫侖屈服一起,用一個(gè)類似GSH的勢(shì)能模型來解決[15].由于各種能量耗散已經(jīng)在細(xì)觀接觸的層面上明確定義,我們相信基于它的DEM計(jì)算能更好地描述長(zhǎng)時(shí)間接觸力鏈的熱損耗和彈性勢(shì)能及其演化(今后擬開展的工作).本文的第三節(jié)將針對(duì)二體力特有的碰撞恢復(fù)系數(shù)隨碰撞初速度變化問題,做些簡(jiǎn)單的數(shù)值計(jì)算,以表明熱力學(xué)方法能夠計(jì)入這個(gè)實(shí)驗(yàn)現(xiàn)象.目前常用的Hertz-Mindlin和Luding模型的恢復(fù)系數(shù)都是材料常數(shù),這個(gè)變化關(guān)系長(zhǎng)期以來一直未能被描述(參見文獻(xiàn)[2]).鑒于能量和耗散與接觸力的關(guān)系比較復(fù)雜,我們將盡可能地給出詳細(xì)解釋和討論,以方便讀者了解熱力學(xué)能量方法與以往的純力學(xué)唯象方法的區(qū)別(第四節(jié)).

2 顆粒作用力模型

如圖1,兩個(gè)質(zhì)量m1,m2,半徑a1,a2的球形顆粒,實(shí)驗(yàn)室坐標(biāo)下的質(zhì)心位矢是r1,r2,相互作用力f.記相對(duì)位矢r=r1?r2,(法向)重疊量δ=d?r,其中r=|r|,d=a1+a2.兩球位置的連線方向,即法向單位矢量n=r/r=(nx,ny).切向單位矢量t=(?ny,nx).記f的法向和切向投影為fn,t,有

為簡(jiǎn)單起見,本文假設(shè)圖1的所有矢量都位于二維平面,并且忽略轉(zhuǎn)動(dòng)(即局限于考慮沒有轉(zhuǎn)動(dòng)的二維運(yùn)動(dòng),這些簡(jiǎn)化不影響將討論的內(nèi)容和概念).這時(shí)兩球的運(yùn)動(dòng)將由牛頓方程

圖1 兩個(gè)運(yùn)動(dòng)顆粒的坐標(biāo)和力示意圖Fig.1.Coordinates and forces of two grains.

所支配.其中的v1,2=dr1,2/dt是速度,f1,2外力是兩球受到的外力(考慮為已知,例如重力或邊界墻力).另外記兩球的相對(duì)速度為v=v1?v2,它的法向投影vn=vn,切向投影vt=vt.兩球的切向相互位移是δt.用點(diǎn)表示時(shí)間導(dǎo)數(shù),例如˙δ=dδ/dt.由(2)和(3)式可得相對(duì)速度的運(yùn)動(dòng)方程:

其中m=m1m2/(m1+m2)是折合質(zhì)量.

顯然一旦有了力fn,t=fn,t(r1,r2,v1,v2)的模型公式,上面的牛頓方程組就閉合了(外力是已知的).所謂的軟球DEM模擬,就是在合適的邊界條件下,對(duì)這樣二體相互作用的大量顆粒系統(tǒng)進(jìn)行數(shù)值求解.注意fn,t的建模不是一件容易的事情,因?yàn)槌吮Ｊ氐膹椈尚妥饔昧ν?還需要考慮諸如庫侖屈服打滑(slip)或塑性(plasticity)、黏力(dashpot)、滯迴(hysteresis)等一系列的不可逆耗散現(xiàn)象.我們面臨的不是一個(gè)哈密頓保守力學(xué)問題,而是復(fù)雜的動(dòng)-靜摩擦力模型.另外根據(jù)赫茲的分析,法向力fn中的彈簧力部分肯定不是線性的.如引言里提到的,目前常用的DEM力模型的一個(gè)問題是能量的耗散情況不明確.例如含黏力Hertz-Mindlin模型:

其中的kn,t,ηn,t,μ是材料常數(shù).文獻(xiàn)中并沒有明確給出相應(yīng)的彈簧勢(shì)能w彈的計(jì)算公式.這時(shí)如果期望分析計(jì)算力鏈的彈性勢(shì)能就沒有標(biāo)準(zhǔn)了.也許一個(gè)最好的做法是取

沒有滑動(dòng)時(shí),它對(duì)δt的導(dǎo)數(shù)的確是ft公式(6)的彈性部分,但出現(xiàn)打滑就不好說了.另外它對(duì)δ的導(dǎo)數(shù)也不是fn公式(5)的彈性部分.總之對(duì)上述Hertz-Mindlin模型,我們覺得是不好定義彈簧勢(shì)能.實(shí)際上Hertz-Mindlin模型可解釋為用Boussinesq力描述法向和切向彈簧力,再加上黏力和庫侖屈服.由于Boussinesq力不是保守力,因而無從知道相應(yīng)的彈性勢(shì)能.為改善以往顆粒-顆粒作用力模型的這個(gè)能量缺陷,下面將模仿GSH的二級(jí)不可逆耗散圖像(見圖2),提出一套基于能量守恒的建模方法.

一般而言如果要求彈性能w彈只是幾何相對(duì)位置r(t)的代數(shù)函數(shù),兩球之間不會(huì)有滑動(dòng)或塑性.也就是說含塑性變形的運(yùn)動(dòng)學(xué)(kinematic)變量不能直接用來描述彈性,需要額外引入表征接觸處彈性變形或彈簧長(zhǎng)度變化的位移矢量:

圖2 顆粒-顆粒作用時(shí)的二級(jí)不可逆熱力學(xué)示意圖,圖中的三個(gè)能量之和是守恒的Fig.2.Two-stage-irreversibility of grain-grain interaction.Sum of the three energies is conserved.

作為w彈的自變量(注意對(duì)本文考慮的二維運(yùn)動(dòng),只討論矢量的兩個(gè)分量即可),即

彈 性 能w彈的 導(dǎo) 數(shù)πn=??w彈/?un和πt=??w彈/?ut將給出顆粒-顆粒相互作用力fn,t的彈性部分.可以不失一般地將它們寫成:

這里的上標(biāo)“熱”和“松”,分別表示圖2中左右兩種耗散機(jī)制的力貢獻(xiàn).類似地可將彈性變形un,ut的運(yùn)動(dòng)方程寫成

其中R松≥0是它的激發(fā)功率,I≥0是它向熱能弛豫的衰減功率.松動(dòng)能的激發(fā)功率是與πn,t的乘積,加上與vn,t的乘積:

相應(yīng)地,有熱功率

其中的I是來自緩沖區(qū)w松的貢獻(xiàn).最后,圖2中機(jī)械能的定義是兩個(gè)顆粒的動(dòng)能與彈性能之和:

當(dāng)沒有外力時(shí),能量守恒要求機(jī)械能w機(jī)械、松動(dòng)能w松和熱能之和保持常數(shù),或它們的時(shí)間導(dǎo)數(shù):

這可以直接計(jì)算證明.的確微分方程(16)有

代入(9)和(10)式,有

將(18)和(13),(15)式相加,即得(17)式.

與以往的模型對(duì)比,基于方程(7)—(17)理論構(gòu)架下的接觸相互作用除了給出彈性勢(shì)能和熱功率外,在處理彈性變形量與幾何運(yùn)動(dòng)學(xué)量的關(guān)系上,以及機(jī)械能與熱能的耗散關(guān)系上,引入了不那么僵硬的微分方程和能量緩沖區(qū)等措施.這將使得理論具有足夠的靈活性來描述前面提到的各種復(fù)雜的接觸力學(xué)現(xiàn)象(庫侖屈服、打滑、塑性、黏力、滯迴等).

將取類似文獻(xiàn)[15,16]的彈性勢(shì)能模型

其中k,ξ,c是材料常數(shù),c是描述濕黏力的常數(shù)[16].注意un和ut分別是法向和切向彈簧長(zhǎng)度變化,量綱是米.彈性能量w彈量綱焦耳.彈簧系數(shù)k量綱為這里考慮了非線性的赫茲3/2冪率彈簧力.勢(shì)能(19)給出的彈簧力是

這里要求法向彈簧的長(zhǎng)度變化總是小于零:un＜0,但剪切彈簧的ut可正可負(fù).后面將看到un＜0可以由熱力學(xué)穩(wěn)定來保障.

松動(dòng)能的衰減可以寫成弛豫時(shí)間模型的形式

(9),(10)和(11),(12)式中的耗散項(xiàng),可以用標(biāo)準(zhǔn)的Onsager非平衡熱力學(xué)來處理.仿照GSH的做法,有

是非對(duì)角(off-diagonal)遷移系數(shù).如果把上面的矩陣方程寫成非矩陣形式,有f松=?η松v,即類似地,有和和將它們代入前面的(7)—(17)式,得顆粒-顆粒相互作用力為:

彈簧變形量un,t的運(yùn)動(dòng)方程為:

松動(dòng)能的運(yùn)動(dòng)方程為

另外熱功率為

將(28)—(32)與牛頓方程(2),(3)聯(lián)立求解,可計(jì)算兩個(gè)顆粒運(yùn)動(dòng)軌跡.

以上接觸理論的材料參數(shù)是:這里除彈性勢(shì)能中的三個(gè)參數(shù)k,ξ,c必須是常數(shù)外,其他參數(shù)都是與耗散有關(guān)的遷移系數(shù),并且可以是松動(dòng)能w松、重疊量δ等變量的函數(shù).考慮到接觸力的復(fù)雜性,在這里保留了較多的遷移系數(shù),以便理論有足夠的靈活性來應(yīng)對(duì).

一個(gè)從GSH繼承來的重要性質(zhì),是非線性彈性勢(shì)能(19)有失穩(wěn)現(xiàn)象,其對(duì)應(yīng)的恰好是庫侖屈服[15].勢(shì)能(19)的力學(xué)穩(wěn)定區(qū)域由不等式

給出,等號(hào)是穩(wěn)定區(qū)域的邊界方程.另外在穩(wěn)定區(qū)域里總是有un＜0,故(19)—(21)不會(huì)出現(xiàn)虛數(shù)(因?yàn)閷?shí)際運(yùn)動(dòng)總是發(fā)生在穩(wěn)定區(qū)域內(nèi)).利用(20),(21)式可將(35)式寫成

失穩(wěn)意味著一旦觸及(35)或(36)式,接觸力的彈性部分將快速松動(dòng).這個(gè)現(xiàn)象可以方便地通過讓在屈服邊界急劇增大來描述.

除失穩(wěn)外,材料參數(shù)(34)式的進(jìn)一步建模還需要顧及重疊量δ的影響,因?yàn)榻佑|力、熱功率、松動(dòng)能等都是重疊時(shí)δ＞0才能被激發(fā)而出現(xiàn)的東西.一旦顆粒分離δ＜0,它們的激發(fā)源都會(huì)消失,同時(shí)也迅速衰減為零(假設(shè)沒有濕黏力c=0).下面試用一組簡(jiǎn)單的材料參數(shù)值,以討論碰撞恢復(fù)系數(shù)問題為例,對(duì)此做些具體說明.

3 碰撞恢復(fù)系數(shù)

作為一個(gè)簡(jiǎn)單例子,取下面的參數(shù)模型

和

其中Θ(x)是Heaviside階梯函數(shù),d=d1=d2和m=m1=m2分別是顆粒直徑和質(zhì)量(考慮兩個(gè)相同顆粒).(37)式中假設(shè)了沒有濕黏力,另外盡量略去了一些遷移系數(shù).模型(38)意味著彈簧變形的激發(fā)源,即(30)和(31)式右邊的第一項(xiàng)只在重疊時(shí)出現(xiàn),分離后它們是零.黏滯系數(shù)(39)式也是這樣.與GSH類似,(40)式假設(shè)切向彈簧變形的弛豫系數(shù)比例于(相當(dāng)于GSH里的顆粒溫度Tg).(41)式在滿足穩(wěn)定條件(35)時(shí)等于零,一旦不滿足穩(wěn)定條件,它將取適當(dāng)?shù)恼?起到盡快松弛切向彈簧變形、讓系統(tǒng)迅速返回穩(wěn)定區(qū)域的目的.注意上述模型的法向耗散與切向耗散有很大差異,后者比前者多考慮了(40)和(41)式兩個(gè)系數(shù).這顯然有其合理性,因?yàn)榍邢蚝纳⒌拇_比法向的復(fù)雜,需要更多的遷移系數(shù)來描寫.另外δ和un,t有長(zhǎng)度量綱,為了讓Heaviside函數(shù)的自變量是量綱1的,用顆粒粒徑d對(duì)它們做了約化.另外發(fā)生顆粒重疊和屈服的數(shù)學(xué)描寫總會(huì)涉及不連續(xù)性,需要用Heaviside函數(shù)來表達(dá).這里是將它們反映在遷移系數(shù)(38),(39)和(41)式中,并且后者與彈性勢(shì)能失穩(wěn)判據(jù)(35)式有關(guān).這個(gè)基于熱力學(xué)概念的描述方法,與直接在力模型中唯象描寫的傳統(tǒng)做法是有差異的.

圖3是上述參數(shù)值下數(shù)值計(jì)算的兩個(gè)等質(zhì)量和直徑的、沿著x方向相對(duì)運(yùn)動(dòng)的顆粒做對(duì)心碰撞的情況.碰撞前后的速度分別是±v0和±vf(都沿著x方向),并且都是分離的,fn,t,un,t,w松都是零.顆粒重疊作用時(shí),它們的激發(fā)、衰減、以及圖2中三個(gè)能量之間的轉(zhuǎn)換情況等,都能夠明確計(jì)算(見圖3).注意由于熱功率R≥0,圖3(c)的機(jī)械能(黑色曲線)總是隨時(shí)間減小或不變.如果忽略碰撞過程細(xì)節(jié),可以將這個(gè)接觸作用模型看作以e=|vf/v0|為恢復(fù)系數(shù)的非彈性碰撞.

在分析長(zhǎng)時(shí)間接觸的力鏈演化問題時(shí),除顆粒初始位置和速度外,還需要un,t,w松的初始值.通?？梢匀〕跏挤ㄏ驈椈勺冃蝩n等于初始重疊量δ.另外任何接觸力模型(包括Hertz-Mindlin和Luding)都會(huì)出現(xiàn)Heaviside階梯函數(shù)Θ(x). 數(shù)值分析時(shí)可用近似的解析函數(shù)來代替,例如取

本文的熱力學(xué)接觸力模型的恢復(fù)系數(shù)e=|vf/v0|不是材料常數(shù),而是初速度v0和(斜碰撞(oblique impact)時(shí)的)碰撞長(zhǎng)度的函數(shù).對(duì)心碰撞時(shí),本節(jié)考慮的模型的恢復(fù)系數(shù)隨初速度增加而減小(圖4i).這與實(shí)驗(yàn)定性符合[2],表明熱力學(xué)建模方法有描述復(fù)雜觀測(cè)細(xì)節(jié)的潛力.

圖3 顆粒-顆粒對(duì)心碰撞過程的時(shí)間演化曲線 (a)位置坐標(biāo);(b)接觸力(c)彈簧勢(shì)能we,動(dòng)能wk和機(jī)械能w;(d)松動(dòng)能wrelaxation;是兩個(gè)顆粒的初始能量(動(dòng)能)Fig.3.Temporal evolution of normal collision of two grains:(a)Relative position;(b)contact force f=(c)potential of springs we,kinematic en-ergy wk,mechanic energy w;(d)relaxation energywrelaxation.The notation is initial(kine-matic)energy of two grains.

本節(jié)計(jì)算的目的僅限于用一個(gè)簡(jiǎn)單例子對(duì)熱力學(xué)建模方法做些具體解釋.實(shí)際顆粒的參數(shù)情況不一定是(37)—(41)式那樣簡(jiǎn)單,參數(shù)值的標(biāo)定也需要結(jié)合實(shí)驗(yàn)仔細(xì)處理.另外出于方便解釋概念和方法,本文忽略了三維和顆粒轉(zhuǎn)動(dòng)運(yùn)動(dòng)φ(t).但在分析實(shí)際問題時(shí),它們是應(yīng)該計(jì)入的.

圖4 對(duì)心碰撞恢復(fù)系數(shù)隨初速度變化曲線Fig.4.Variation of normal restitution coefficient with collision velicity.

4 討 論

為方便討論本文熱力學(xué)模型與其他接觸力模型的概念差異,我們將彈性能方程(8)給出的力,和方程(28)—(31)寫成下面的矢量形式:

1)彈性力π總是要求為保守力.也就是說有明確定義的彈性勢(shì)能w彈(u).目前的DEM計(jì)算,雖然有些采用了保守的線性法向和線性切向彈簧,但有些是非保守的(例如Boussinesq力).也就是說對(duì)保守性的要求并不像熱力學(xué)那樣,是強(qiáng)制性的.非保守的彈簧模型會(huì)導(dǎo)致機(jī)械功率和熱功率之間的分配比例不能清楚定義的問題.

2)彈簧應(yīng)變u并不是純幾何的位矢r=∫vdt,而是由運(yùn)動(dòng)方程(43)描述的獨(dú)立變量.這在目前的DEM計(jì)算中都注意到了(特別是切向分量),但u的運(yùn)動(dòng)方程與(43)式有很大不同.

3)接觸力方程(44)中的黏力在當(dāng)前DEM計(jì)算中大都有類似的對(duì)應(yīng)(例如(5),(6)中的ηt,n項(xiàng)).但在對(duì)黏滯系數(shù)?η的處理方式上,沒有圖2所示的二級(jí)不可逆概念背景.

4)本文的接觸力和幾何運(yùn)動(dòng)學(xué)速度(f,v),與彈簧力和彈簧變形(π,u)之間的關(guān)系方程(43),(44)采用了典型的Onsager非平衡熱力學(xué)形式.但這里除了熟知的機(jī)械能和熱能外,還必須考慮一個(gè)中間緩沖能量w松去控制某些遷移過程和系數(shù).這個(gè)稱作“二級(jí)不可逆”的方案是顆粒物質(zhì)熱力學(xué)的基本特征.它不僅對(duì)宏觀連續(xù)力學(xué)層面的GSH有效[11?14],細(xì)觀層面的顆粒接觸熱力學(xué)也必須這樣處理.它能把庫侖屈服、彈性松弛、力與幾何變形之間的滯迴和棘輪(ratcheting)等復(fù)雜力學(xué)現(xiàn)象,直接自然地與遷移過程聯(lián)系起來,使得用熱力學(xué)語言來理解和分析它們成為可能.從物理角度看這顯然是合理的,因?yàn)轭w粒物質(zhì)的這些復(fù)雜力學(xué)現(xiàn)象都源自于耗散發(fā)熱.

5)方程(43),(44)右邊的第二項(xiàng)可分別理解為彈性弛豫和黏滯.右邊的第一項(xiàng)涉及到以往未曾遇到的非對(duì)角遷移系數(shù)?G或?α=?1??G.如果忽略其他遷移系數(shù),可將?α單獨(dú)看作一個(gè)具有改變兩端力和形變比值效果的“力學(xué)元件”.我們把它形象地稱作“變速箱(gearbox)”,其能量和力學(xué)關(guān)系情況如圖5所示.值得指出的是,如果變速比?G恒定(α=const.)并且u=(1?α)ε的話,可將圖中的彈簧勢(shì)能和力寫成幾何變形ε的函數(shù):w彈=(1?α)2kε2/2和f=?(1?α)2kε,它們滿足熱力學(xué)關(guān)系f=??w彈/?ε.這意味著如果不考慮α的動(dòng)態(tài)效應(yīng),變速箱是不需要的,因?yàn)樗淖饔弥皇钦{(diào)整與其聯(lián)接的彈簧的彈性系數(shù),完全可以將其歸入彈簧元件中.這也許是以往的各種力學(xué)元件都不曾提到變速箱的原因.但如果α是動(dòng)態(tài)的,在考慮具有熱力學(xué)背景的力學(xué)理論時(shí),變速箱與彈簧就必須區(qū)別開來了.變速箱這個(gè)非對(duì)角遷移系數(shù)是我們近年來為理解顆粒物質(zhì)宏觀彈塑轉(zhuǎn)變行為而提出的[17],它的統(tǒng)計(jì)物理和動(dòng)理學(xué)意義還有待深入研究.唯象地看它的力學(xué)效果是:“在不增加熱損失的前提下減弱彈性”．目前可以肯定的是,作為彈應(yīng)變u的驅(qū)動(dòng)源項(xiàng)的系數(shù)(見(43)式),它與彈性的出現(xiàn)與消亡密切相關(guān),而且是動(dòng)態(tài)的．本文為了簡(jiǎn)單,假設(shè)α隨著由重疊量控制的接觸彈性的出現(xiàn)與消亡,在0和1之間跳變,見模型(38)．當(dāng)然實(shí)際情況可以復(fù)雜許多,例如還與松動(dòng)能w松有關(guān),在0和1之間變化等．

6)GSH熱力學(xué)的一個(gè)重要特征是圖2所示的二級(jí)不可逆耗散現(xiàn)象:除了左邊的直接耗散為熱的過程外,機(jī)械能還可以通過右邊的緩沖能量w松進(jìn)行耗散,同時(shí)遷移系數(shù)也可以受其控制(即?G,?λ,?η可以是w松的函數(shù))．需要注意的是,庫侖屈服現(xiàn)象雖然源自于彈性勢(shì)能的靜力失穩(wěn)(見(35),(36)),我們還需要在某些遷移系數(shù)上采取相應(yīng)的措施來保障任何動(dòng)力學(xué)演化都不會(huì)進(jìn)入該非穩(wěn)區(qū)域,例如模型(41)．(因?yàn)槿魏斡^測(cè)到的實(shí)際系統(tǒng)都在熱力學(xué)穩(wěn)定區(qū)域內(nèi)運(yùn)動(dòng)著,不會(huì)進(jìn)入非穩(wěn)區(qū)域．這與廣為熟悉的范德瓦耳斯氣體理論類似,它的熱力學(xué)狀態(tài)空間里也有一個(gè)實(shí)際氣體不可能進(jìn)入的非穩(wěn)區(qū)域)．

圖5 變速箱的能量與力學(xué)關(guān)系(注意α可以是動(dòng)態(tài)的,但m,k假設(shè)為常數(shù))Fig.5.Energetics and mechanics of gearbox.Notably α may temporally change,but m,k are assumed constant.

由于兩個(gè)顆粒之間的很多相互作用行為是可以實(shí)驗(yàn)觀測(cè)的,任何接觸力模型原則上都應(yīng)該用這些細(xì)觀顆粒尺度的測(cè)量結(jié)果進(jìn)行驗(yàn)證和標(biāo)定．但以往的一些DEM工作發(fā)現(xiàn),這個(gè)細(xì)觀-宏觀聯(lián)系并不都是協(xié)調(diào)的,例如用線性彈簧計(jì)算的宏觀行為,卻比用更合理的赫茲彈簧要好(參見文獻(xiàn)[12]的第六章)．也許這些矛盾和困難與所采用的接觸力模型對(duì)熱力學(xué)概念的描述不夠全面有關(guān)．如果是這樣的話本文的熱力學(xué)建模方法將有助于問題的澄清．

目前顆粒物理研究方法可分為純幾何、DEM接觸力、和GSH熱力學(xué)三種．幾何法致力于僅從顆粒的運(yùn)動(dòng)學(xué)變量(r,φ)出發(fā),通過深入考察其幾何時(shí)空結(jié)構(gòu)與關(guān)聯(lián),以及幾何構(gòu)型熵(con fi guration entropy)或漲落耗散關(guān)系( fl uctuation-dissipation relation)等,來解釋顆粒物質(zhì)的各種物性(參見文獻(xiàn)[18—21])．也就是說期望構(gòu)建類似開普勒行星運(yùn)動(dòng)理論那樣的顆粒物理．傳統(tǒng)DEM方法的出發(fā)點(diǎn)是{r,φ,f},即比純幾何法增加了接觸力f,并且期望在牛頓方程的基礎(chǔ)上解釋顆粒物質(zhì)物性．由于接觸力f不是可以用哈密頓描寫的保守力,DEM一般都是直接數(shù)值求解大量顆粒的牛頓方程組．如果忽略f中的不可逆,僅保留保守的彈性力,也可以用傳統(tǒng)的平衡態(tài)統(tǒng)計(jì)物理來分析這個(gè)哈密頓顆粒系統(tǒng)[22]．另外對(duì)f中的不可逆部分,還可以嘗試用動(dòng)理學(xué)[23]或mode-coupling[24]等非平衡統(tǒng)計(jì)技術(shù)來分析．但隨著顆粒數(shù)密度增加,這些非平衡統(tǒng)計(jì)技術(shù)都面臨困難,目前還無法滿意地處理顆粒固體,特別是集體彈性行為的納入和描述．因此數(shù)值求解牛頓方程仍是當(dāng)前分析接觸作用系統(tǒng)的最好辦法(如引言里提到的,它能描述宏觀彈性的呈現(xiàn)于消亡)．注意幾何法和傳統(tǒng)DEM法都可以不具體涉及熱力學(xué)熵,可稱作“非熱(athermal)”類型理論．與之對(duì)比,熱力學(xué)方法則強(qiáng)制性地要求無論是連續(xù)力學(xué)的宏觀尺度,還是本文考慮的顆粒接觸力細(xì)觀尺度,熱力學(xué)能量守恒和熵增加原理都應(yīng)該得到明確的體現(xiàn)和遵守．當(dāng)然三種方法之間有共同相似或互補(bǔ)的內(nèi)容,但肯定還有不能互通的內(nèi)容,因?yàn)榕ｎD力不可能取代熱力學(xué),也不可能被純幾何學(xué)取代．考慮到耗散發(fā)熱是顆粒物質(zhì)的基本性質(zhì),物理熱力學(xué)應(yīng)該是最為合理和可靠的方法．

值得提到的是,固體顆粒間的法向接觸力一般沒有塑性打滑等現(xiàn)象,模型(37)中取很可能是普遍合理的．但在Luding模型里法向滯迴是他的特色[5]．為方便今后與其對(duì)比,我們?cè)贠nsager矩陣(23)—(26)里保留了這幾個(gè)遷移系數(shù)．遷移系數(shù)那些是無價(jià)值(irrelevant)可以去掉的,那些是有價(jià)值(relevant)需要保留的,及其相關(guān)的動(dòng)力學(xué)意義等,目前還不能確定．這些問題需要結(jié)合DEM計(jì)算和實(shí)驗(yàn)情況來逐步澄清．模型(37)—(41)的參數(shù)數(shù)值僅作為參考,沒有普適性．它們都是材料參數(shù),應(yīng)該依據(jù)實(shí)驗(yàn)來標(biāo)定．另外將本文理論向三維推廣時(shí),只需要對(duì)文中的力、速度、位移等矢量補(bǔ)充一個(gè)切向分量即可．對(duì)球形顆粒,由于對(duì)稱性兩個(gè)切向的遷移系數(shù)是一樣的,理論參數(shù)仍由(34)式列出．如果考慮(37)—(41)式的模型,這個(gè)三維理論給出的二維運(yùn)動(dòng)將退化為本文的公式．這時(shí)它給出的對(duì)心碰撞也就是本文的圖3和圖4．顆粒轉(zhuǎn)動(dòng)在一些接觸力模型有所考慮,例如Luding模型[5]．但也有不少DEM工作省略了這個(gè)復(fù)雜性,例如文獻(xiàn)[9]．本文出于簡(jiǎn)單和方便討論熱力學(xué)概念的原因也省略了顆粒轉(zhuǎn)動(dòng)．它的計(jì)入要復(fù)雜許多,建議參照Luding的做法．

常用接觸力模型與本文建議的熱力學(xué)方法有很多對(duì)應(yīng)的和相似的地方,例如彈簧和黏滯阻尼．算例模型(37)—(41)式中的彈簧系數(shù)k(對(duì)應(yīng)于剛度系數(shù)),也是常用模型中經(jīng)常出現(xiàn)的常數(shù)．但熱力學(xué)肯定還有自己的不同內(nèi)容和改善的地方,例如將屈服和打滑放到了遷移系數(shù)中、給出了機(jī)械能和熱功率的具體表達(dá)式、要求伴隨任何遷移耗散過程的熱功率都有“二次正定多項(xiàng)式”的形式(見(32),(33)式)等．注意“二次正定多項(xiàng)式”的要求來自熱力學(xué)穩(wěn)定性,它也會(huì)改善模型的數(shù)值穩(wěn)定性．常用模型一般只有黏性的耗散熱是平方正定的形式,其他復(fù)雜的不可逆現(xiàn)象往往沒有這個(gè)性質(zhì),例如屈服．顯然熱力學(xué)內(nèi)容的最大優(yōu)勢(shì)將體現(xiàn)在研究能量和耗散問題方面,這在最近已經(jīng)開始受到重視(例如Yu小組的工作[25,26])．

本文建模方法與常用接觸力模型的仔細(xì)對(duì)比,特別是具體改善和優(yōu)越的地方,將在今后的定量對(duì)比工作中具體體現(xiàn)出來．這里可分為物理和工程兩個(gè)方面來開展．物理方面,僅關(guān)心“二體力”這個(gè)基礎(chǔ)層面即可,也就是說只是對(duì)比不同方法和模型的“兩個(gè)顆粒動(dòng)力學(xué)”情況即可,這樣便于討論和澄清與基本概念相關(guān)的優(yōu)劣對(duì)比,同時(shí)避開大規(guī)模數(shù)值編程的麻煩．工程方面,主要考慮大量顆粒的DEM計(jì)算情況的對(duì)比,涉及大規(guī)模數(shù)值編程或商業(yè)軟件．另外“二體力”和“DEM”兩邊都應(yīng)該盡可能多地考慮與實(shí)驗(yàn)的對(duì)比．

在考慮計(jì)入顆粒轉(zhuǎn)動(dòng)因素時(shí),首先是增加描述顆粒轉(zhuǎn)動(dòng)角的運(yùn)動(dòng)方程(即牛頓力矩方程),在機(jī)械能公式(16)中增加顆粒轉(zhuǎn)動(dòng)動(dòng)能,以及可能的扭力勢(shì)能等．然后是滾動(dòng)力(rolling)和扭轉(zhuǎn)力(torsion)的建模問題．由于熱力學(xué)的影響主要體現(xiàn)在不可逆耗散部分,建模時(shí)涉及的純幾何運(yùn)動(dòng)學(xué)(kinematic)部分,以及滾動(dòng)力和扭轉(zhuǎn)力的可逆部分,基本可直接挪用Luding的相應(yīng)公式[5]．但滾動(dòng)力和扭轉(zhuǎn)力的不可逆部分,需要用與本文類似的能量守恒和Onsager遷移系數(shù)等方法來處理．

5 結(jié) 論

雖然耗散發(fā)熱是顆粒間接觸運(yùn)動(dòng)時(shí)的重要現(xiàn)象,當(dāng)前DEM采用的接觸力模型仍屬于唯象力學(xué)類型,對(duì)熱力學(xué)要求的機(jī)械功率和熱功率等的定義并不十分清楚．本文的主要結(jié)論是:基于熱力學(xué)的接觸力建模方法可以參照已有的GSH概念來實(shí)現(xiàn)．我們相信從唯象力學(xué)發(fā)展到熱力學(xué)模型,不僅僅從物理概念上更加合理,對(duì)諸如恢復(fù)系數(shù)行為等的一些實(shí)驗(yàn)現(xiàn)象的描述,也會(huì)有所改善．當(dāng)然在嘗試將本文內(nèi)容用于DEM數(shù)值計(jì)算之前,還需要先將其推廣到三維和計(jì)入顆粒轉(zhuǎn)動(dòng)運(yùn)動(dòng)．這些將是今后需要陸續(xù)開展的工作．

感謝厚美瑛、孫其誠(chéng)、程曉輝、劉曉星和張國(guó)華的討論.

[1]Hertz H 1881J.Reine Angew.Math.92 156

[2]Thornton C 2015Granular Dynamics,Contact Mechanics and Particle System Simulations—A DEM Study(Particle Technology Series,Volume 24)(eBook DOI 10.1007/978-3-319-18711-2)(Switzerland:Springer International Publishing AG)

[3]Laughlin R B,Pines D 2000Proc.Natl.Acad.Sci.USA97 28

[4]Laughlin R B 2004A Different Universe(Changsha:Hunan Science and Technology Press)(in Chinese)[王文浩譯2008不同的宇宙(長(zhǎng)沙:湖南科學(xué)技術(shù)出版社)]

[5]Luding S 2008Granular Matter10 235

[6]Truesdell C 1972Rational Thermodynamics(Berlin:Springer-Verlag)

[7]Boussinesq J 1873C.R.Hebd.Seances Acad.Sci.77 1521

[8]Bonneau L,Andreotti B,Clément E 2007Phys.Rev.E75 016602

[9]Sun Q C,Jin F,Zhou G D 2013Granular Matter15 119

[10]Kumar N 2014Ph.D.Dissertation(Enschede:University of Twente.The Netherlands ISBN:978-90-365-3634-9,DOI:10.3990/1.9789036536349)

[11]Jiang Y M,Liu M 2009Granular Matter11 139

[12]Sun Q C,Hou M Y,Jin F 2011Physics and Mechanics of Granular Matter(Beijing:Science Press)(in Chinese)[孫其誠(chéng),厚美瑛,金峰 2011顆粒物質(zhì)物理與力學(xué)(北京:科學(xué)出版社)]

[13]Jiang Y M,Liu M 2014Acta Mech.225 2363

[14]Jiang Y M,Liu M 2015Eur.Phys.J.E38 15

[15]Jiang Y M,Liu M 2003Phys.Rev.Lett.91 144301

[16]Jiang Y M,Liu M 2004Phys.Rev.Lett.93 148001

[17]Jiang Y M,Liu M 2007Phys.Rev.Lett.99 105501

[18]Torquato S,Stillinger F H 2010Rev.Mod.Phys.82 2633

[19]Edwards S F,Moun field C C 1996Physica A226 1

[20]Edwards S F,Moun field C C 1996Physica A226 12

[21]Edwards S F,Moun field C C 1996Physica A226 25

[22]Parisi G,Zamponi F 2010Rev.Mod.Phys.82 789

[23]Luding S 2009Nonlinearity22 R101

[24]Hayakawa H,Otsuki M 2008Prog.Theor.Phys.119 381

[25]Kuang S B,Zou R P,Pan R H,Yu A B 2012Ind.Eng.Chem.Res.51 14289

[26]Hou Q F,Kuang S B,Yu A B 2017Chem.Engineer.Sci.161 67