顆粒-顆粒接觸力的熱力學模型?

蔣亦民 劉佑

1)(中南大學物理與電子學院,長沙 410083)

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

1 引 言

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

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

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

2 顆粒作用力模型

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

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

圖1 兩個運動顆粒的坐標和力示意圖Fig.1.Coordinates and forces of two grains.

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

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

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

其中的kn,t,ηn,t,μ是材料常數.文獻中并沒有明確給出相應的彈簧勢能w彈的計算公式.這時如果期望分析計算力鏈的彈性勢能就沒有標準了.也許一個最好的做法是取

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

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

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

作為w彈的自變量(注意對本文考慮的二維運動,只討論矢量的兩個分量即可),即

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

這里的上標“熱”和“松”,分別表示圖2中左右兩種耗散機制的力貢獻.類似地可將彈性變形un,ut的運動方程寫成

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

相應地,有熱功率

其中的I是來自緩沖區w松的貢獻.最后,圖2中機械能的定義是兩個顆粒的動能與彈性能之和:

當沒有外力時,能量守恒要求機械能w機械、松動能w松和熱能之和保持常數,或它們的時間導數:

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

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

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

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

將取類似文獻[15,16]的彈性勢能模型

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

這里要求法向彈簧的長度變化總是小于零:un＜0,但剪切彈簧的ut可正可負.后面將看到un＜0可以由熱力學穩定來保障.

松動能的衰減可以寫成弛豫時間模型的形式

(9),(10)和(11),(12)式中的耗散項,可以用標準的Onsager非平衡熱力學來處理.仿照GSH的做法,有

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

彈簧變形量un,t的運動方程為:

松動能的運動方程為

另外熱功率為

將(28)—(32)與牛頓方程(2),(3)聯立求解,可計算兩個顆粒運動軌跡.

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

一個從GSH繼承來的重要性質,是非線性彈性勢能(19)有失穩現象,其對應的恰好是庫侖屈服[15].勢能(19)的力學穩定區域由不等式

給出,等號是穩定區域的邊界方程.另外在穩定區域里總是有un＜0,故(19)—(21)不會出現虛數(因為實際運動總是發生在穩定區域內).利用(20),(21)式可將(35)式寫成

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

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

3 碰撞恢復系數

作為一個簡單例子,取下面的參數模型

和

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

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

在分析長時間接觸的力鏈演化問題時,除顆粒初始位置和速度外,還需要un,t,w松的初始值.通常可以取初始法向彈簧變形un等于初始重疊量δ.另外任何接觸力模型(包括Hertz-Mindlin和Luding)都會出現Heaviside階梯函數Θ(x). 數值分析時可用近似的解析函數來代替,例如取

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

圖3 顆粒-顆粒對心碰撞過程的時間演化曲線 (a)位置坐標;(b)接觸力(c)彈簧勢能we,動能wk和機械能w;(d)松動能wrelaxation;是兩個顆粒的初始能量(動能)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.

本節計算的目的僅限于用一個簡單例子對熱力學建模方法做些具體解釋.實際顆粒的參數情況不一定是(37)—(41)式那樣簡單,參數值的標定也需要結合實驗仔細處理.另外出于方便解釋概念和方法,本文忽略了三維和顆粒轉動運動φ(t).但在分析實際問題時,它們是應該計入的.

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

4 討 論

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

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

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

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

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

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

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

圖5 變速箱的能量與力學關系(注意α可以是動態的,但m,k假設為常數)Fig.5.Energetics and mechanics of gearbox.Notably α may temporally change,but m,k are assumed constant.

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

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

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

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

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

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

5 結 論

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

感謝厚美瑛、孫其誠、程曉輝、劉曉星和張國華的討論.

[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不同的宇宙(長沙:湖南科學技術出版社)]

[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)[孫其誠,厚美瑛,金峰 2011顆粒物質物理與力學(北京:科學出版社)]

[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