張興剛，曹 磊

(1.貴州大學(xué) 物理學(xué)院,貴州 貴陽 550025;2.中國電建集團貴陽勘測設(shè)計研究院,貴州 貴陽 550081)

顆粒介質(zhì)的靜力學(xué)問題是顆粒物質(zhì)研究中基礎(chǔ)且重要的方面，它與實際應(yīng)用密切相關(guān)(例如，工程中地基的承載能力、砂石骨料的受力結(jié)構(gòu)，材料科學(xué)中顆粒材料的力學(xué)性能，等等)，對于其它復(fù)雜體系(例如，懸浮液、泡沫、膠體等)的理解也有重要的意義[1-2]。人們在很早的時候就發(fā)現(xiàn)，筒倉中的糧食與水在靜力學(xué)方面有很不相同的表現(xiàn)，糧倉底部所受壓力會隨糧食堆積高度的增加而趨于飽和，這就是顆粒物質(zhì)靜力學(xué)中一個著名的效應(yīng)——糧倉效應(yīng)[1，3]。早在1895年，德國工程師Janssen就提出了基于一定假設(shè)的連續(xù)介質(zhì)模型，對糧倉效應(yīng)進行解釋。然而，顆粒物質(zhì)是一種典型的復(fù)雜體系，目前還沒有關(guān)于顆粒介質(zhì)靜力學(xué)的系統(tǒng)、成熟的理論，許多基本問題仍需深入探究。因此，自上世紀(jì)九十年代的顆粒物質(zhì)研究熱潮以來，仍有許多研究者對糧倉效應(yīng)及相關(guān)的問題進行研究[4-13]。例如，最近有關(guān)反糧倉效應(yīng)[11]和磁性顆粒體系糧倉效應(yīng)[12]的工作。已有的研究工作以實驗和計算機模擬為主，理論方面也主要是對Janssen模型、OSL模型等宏觀唯象模型[4]的討論，對于糧倉效應(yīng)產(chǎn)生的微觀機理以及唯象參數(shù)與微觀結(jié)構(gòu)的聯(lián)系還不是很清楚。為了深入認識顆粒介質(zhì)所服從的力學(xué)規(guī)律，很有必要加強顆粒體系的微觀結(jié)構(gòu)、力網(wǎng)、力傳遞等問題的研究，基于微觀的幾何和力學(xué)分析建立體系的微觀特征與宏觀性質(zhì)之間的聯(lián)系[13-15]。

對于圓盤的晶格堆積，文獻[6]中通過結(jié)構(gòu)和力學(xué)平衡的分析得到了相鄰兩層顆粒間力傳遞的規(guī)則，并且利用力傳遞規(guī)則得到了隨著堆積層數(shù)增加容器底部總壓力會趨向飽和的結(jié)果。對于非晶格的堆積結(jié)構(gòu)，不借助計算機模擬很難從理論上直接進行結(jié)構(gòu)和力學(xué)的分析。標(biāo)量q模型是研究顆粒堆積中力的幾率分布的一個著名的模型[15]。為了便于從理論上研究筒倉底部壓力分布與頂部壓力分布的關(guān)系，可以采用q模型的基本思想，將筒倉中隨機堆積的顆粒體系抽象簡化為格點系統(tǒng)。本文以二維筒倉中隨機堆積的顆粒體系為物理背景，建立了一個由吸收系數(shù)p和側(cè)向傳遞系數(shù)q決定的格點系統(tǒng)模型，該模型的解析求解涉及到三對角形式的傳遞系數(shù)矩陣A(p,q)的相似對角化。文中通過相似變換將A(p,q)對稱化，并且提出基于二階差分方程的方法對矩陣的特征值和特征向量進行分析和計算；最后，將相關(guān)的數(shù)學(xué)結(jié)果用于討論糧倉效應(yīng)的有關(guān)問題。

1 力傳遞的數(shù)學(xué)模型

為討論問題方便，這里主要研究二維的情形，將糧倉效應(yīng)中的物理體系抽象為圖1所示的帶有邊界的格點系統(tǒng)。在格點系統(tǒng)中，一個格點等效于一小團顆粒，第m行n列的格點用符號P(m,n)表示，格點的總行數(shù)M≥2，總列數(shù)N≥5，并且每個格點都受到重力m0g的作用。類似于標(biāo)量q模型，這里主要關(guān)注每個作用力在豎直方向上的分量，因此在下面的討論中，所說的作用力一般是指力在豎直方向上分力的大小。頂部各個格點所受的壓力可用行向量f0=(f01，f02,…,f0N)描述，底部各個格點對容器底的壓力用fM=(fM1，fM2,…,fMN)描述。為了給出底部壓力分布fM與頂部壓力分布f0的關(guān)系，需要結(jié)合實際物理體系的特征定義恰當(dāng)?shù)南噜弮蓪痈顸c間力傳遞的規(guī)則。這里以隨機堆積為物理背景，定義了圖1所示的一種典型的力傳遞規(guī)則。規(guī)定同一行的格點間沒有豎直方向的分力，不相鄰的兩行格點間也沒有相互作用。對于任意格點P(m,n)，稱其所受的重力、頂部對其產(chǎn)生的壓力、以及第m-1行格點對其產(chǎn)生的作用力為輸入力；將P(m,n)對第m+1行格點、及其對容器的作用力稱為輸出力；對于每一個格點，只有上述五種豎直方向分力的存在，為了保證格點在豎直方向的力平衡，總輸入力必須等于總輸出力。圖1中用實線箭頭比較完整地畫出了第一行格點的各個輸入力及輸出力的情況(重力沒有畫)，圖中的虛線箭頭上也標(biāo)出了每種輸出力對應(yīng)的系數(shù)。記P(m,n)所受的總輸入力為wmn，總輸入力向量wm=(wm1，wm2,…,wmn)，重力向量b=(m0g，m0g,…,m0g)，有，

圖1 格點系統(tǒng)模型中力的傳遞Fig.1 Rules of force propagation in the lattice system

w1=f0+b=(f01+m0g，f02+m0g,…,

f0N+m0g)

(1)

圖1中定義的力傳遞規(guī)則，用數(shù)學(xué)語言表達就是當(dāng)1≤m≤M-1時，有力傳遞方程

wm+1=wmA+b

(2)

其中，N階矩陣

A(p,q)=

(3)

稱其為傳遞系數(shù)矩陣，它是一個三對角矩陣[16]。傳遞系數(shù)矩陣只由p和q這兩個參量確定，因此為了表述方便，稱這個模型為pq模型。模型中的p稱為吸收系數(shù)，反映了由于容器側(cè)壁與顆粒的摩擦所導(dǎo)致的格點總輸入力被器壁吸收的程度，它滿足0≤p<1；若容器側(cè)壁與顆粒無摩擦，則p=0，有摩擦?xí)r0

由圖1中定義的力傳遞規(guī)則，有fM=wM；因此由(1)、(2)兩式可以導(dǎo)出fM與f0的關(guān)系。假設(shè)存在可逆矩陣Q使A=QΛQ-1，并且假設(shè)A-E可逆(E是N階單位矩陣)，不難導(dǎo)出

fM=f0QΛM-1Q-1+bQ(Λ-E)-1(ΛM-E)Q-1

(4)

于是從理論上計算fM的關(guān)鍵在于考慮傳遞系數(shù)矩陣A的相似對角化以及矩陣A-E的可逆性。

2 矩陣A(p,q)的相似對角化

D-1AD=B(p,q)=

(5)

(6)

為了給出矩陣P與Λ，需要求解B的特征值與特征向量。設(shè)xT=(x1,x2,…,xN)T是B的特征向量，λ是相應(yīng)的特征值，那么應(yīng)該有,

(B-λE)xT=0

(7)

傳統(tǒng)的方法是先由B的特征方程計算特征值，再將特征值代入方程組(7)求相應(yīng)的特征值。與傳統(tǒng)的計算方法不同，這里采用直接分析齊次線性方程組(7)的方法，通過同時進行特征值與特征向量的計算從而簡化理論分析過程。將(5)代入(7)可得帶有邊值條件的差分方程，

(8)

這個差分方程的邊值條件比較復(fù)雜。引入新的變量yT=(y1,y2,…,yN)T，定義可逆變換yT=DxT；由(8)可得到關(guān)于yT的差分方程，

(9)

由于B是對稱矩陣，其特征值都是實數(shù)，因此只需在λ∈的范圍內(nèi)討論方程(9)。可以根據(jù)二階常系數(shù)差分方程的理論求(9)中差分方程的通解。該差分方程對應(yīng)的特征方程為

(10)

(1)第一種情況λ′2=1：方程(10)有二重解r1=r2=λ′，此時(9)中差分方程的通解為

yn=c1λ′n+c2nλ′n,1≤n≤N

(11)

以y1和yN為邊值，那么(11)中的待定系數(shù)c1與c2應(yīng)該滿足，

(12)

這個方程組中系數(shù)矩陣行列式不為0，容易得到向量(c1,c2)T的唯一解，將其解代入(11)式有，

(13)

利用(13)得到y(tǒng)2、yN-1關(guān)于y1、yN的表達式，考慮到λ′2=1并且將它們代入(9)中的邊值條件有，

(14)

記這個關(guān)于向量(y1,yN)T的方程組的系數(shù)矩陣行列式為det(p,q,N,λ′)，顯然λ′2=1所對應(yīng)λ要成為特征值需要滿足方程，

det(p,q,N,λ′)=0

(15)

為方便起見，分兩種情況進行討論：

①當(dāng)λ′=1時，方程(15)可寫為

p(N-1)[2q(1-p)+p(N-1)]=0

(16)

②當(dāng)λ′=-1時，方程(15)可寫為

p(1-2q)(N-1)[p(N-1)-2q(1+p(N-2))]=0

(17)

利用λ=q(λ′-1)+1，(16)、(17)式，以及上述討論過程不難得出下面的結(jié)論。

可將λ′=-1代入(14)式，給出(y1,yN)T的通解；然后利用(13)式得出向量yT，再利用xT=D-1yT并且單位化后就可得到結(jié)論2中相應(yīng)于λ=1-2q的單位特征向量。

(2)第二種情況λ′2≠1：方程(10)有兩個不同的根r1、r2，此時(9)中差分方程的通解為

可以用復(fù)函數(shù)簡化問題的討論，設(shè)自變量z定義在復(fù)平面上，復(fù)余弦與復(fù)正弦函數(shù)分別定義為

(19)

記集合S1={z|lmz=0},S2=Pz|Rez=kπ,k∈},S3=S1∩S2，顯然S1就是復(fù)平面上的實軸，那么有，

定理1當(dāng)且僅當(dāng)z∈S1∪S2時，函數(shù)cosz取實數(shù)值；z∈S3時，cos2z=1；z∈S1-S3時，cos2z<1；z∈S2-S3時，cos2z>1。

定理2函數(shù)cosz與sinz的零點都在實軸上。

運用復(fù)數(shù)的知識不難證明這兩個定理，這里不再贅述。利用根r1、r2的形式并且結(jié)合定理1，可以設(shè)λ′=cosz，不難得出r1=eiz、r2=e-iz；需要注意的是由于λ′∈、λ′2≠1，因此只需在z∈(S1∪S2)-S3的范圍內(nèi)討論問題。于是(9)中差分方程的通解變?yōu)?/p>

yn=c1einz+c2e-inz，1≤n≤N

(20)

以y1和y2為邊值，那么(20)中的待定系數(shù)c1與c2應(yīng)該滿足，

(21)

利用定理2易知該方程組中系數(shù)矩陣行列式在整個z∈(S1∪S2)-S3的范圍內(nèi)都不為零，容易得到向量(c1,c2)T的唯一解，將其解代入(20)式有，

(22)

利用(22)得到y(tǒng)N-1、yN關(guān)于y1、y2的表達式，將它們代入(9)中的邊值條件有，

(23)

這個關(guān)于(y1,y2)T的方程組的系數(shù)矩陣行列式為

(24)

顯然當(dāng)λ′2≠1時矩陣B的特征值就由方程

det(z,p,q,N)=0

(25)

中z的解的情況決定。這是一個復(fù)雜的三角方程，可以通過令z′=ez從而將其轉(zhuǎn)化為關(guān)于z′的高次代數(shù)方程；求出該高次代數(shù)方程的所有解便可得到B(p,q)的所有特征值，進一步可利用(22)式求出所有相應(yīng)的特征向量，從而解決A(p,q)的相似對角化問題。但是，一般情況下只能求該方程的數(shù)值解，因此為了得出解析的結(jié)果只能對一些特殊的情況進行討論。明顯p=1/2、q=1是最簡單、在物理上也屬于比較典型的情形，因為它相當(dāng)于筒倉系統(tǒng)中顆粒與器壁的摩擦不是很小也不是很大，并且顆粒團由于成拱效應(yīng)明顯而幾乎不向其正下方傳遞力的情況。此時方程(25)為

sin(N+1)z=0

(26)

(27)

(28)

3 應(yīng)用與討論

由結(jié)論1，當(dāng)容器側(cè)壁與顆粒有摩擦(即p≠0)時，1不是B(p,q)的特征值，易知這時矩陣A(p,q)-E是可逆的。因此p≠0時可以利用(4)式給出容器底部壓力分布fM與頂部壓力分布f0的關(guān)系。利用上面介紹的方法可以求出B(p,q)的所有特征值和特征向量，于是由(4)、(5)、(6)式可得，

(29)

可以將這個公式用于討論糧倉效應(yīng)等力學(xué)問題。在糧倉效應(yīng)的有關(guān)實驗中，通常關(guān)注容器底部測得的有效質(zhì)量Me(即底部的總壓力與重力加速度g的比值)與容器中顆粒總質(zhì)量Mt的關(guān)系。對于p=1/2、q=1這種典型的情況，如果容器頂部無壓力，那么可由(27)、(28)、(29)式得到，

(30)

其中，fMn是向量fM的第n個分量。在實際問題中更關(guān)心M和N遠大于1時Me隨Mt變化的大致情況。明顯(30)式的求和部分只有當(dāng)k取奇數(shù)時才有貢獻，而且隨著k的增大求和部分會迅速地減小并趨于0；于是對于近似地討論問題可以取一級近似，得到，

(31)

設(shè)M=h(N+1)2，其中h可間接地反映堆積的高度。當(dāng)N>>1時(31)式簡化為

(32)

(33)

利用Janssen模型[1,4]可以得到Me=Ms(1-e-Mt/Ms)，其中Ms叫作飽和質(zhì)量；顯然(33)式的結(jié)果與Janssen模型得到的結(jié)果很接近，這說明本文中的pq模型也可以解釋糧倉效應(yīng)，它在一定程度上正確反映了筒倉系統(tǒng)中力傳遞的特征。利用(29)式，還可以比較方便地討論容器底部的應(yīng)力分布、點載荷的響應(yīng)等方面的力學(xué)問題。當(dāng)然，如果要深入了解p和q對這些力學(xué)問題的影響，就需要求解方程(25)；根據(jù)前面的數(shù)學(xué)討論，在一般情況下應(yīng)該結(jié)合數(shù)值計算才能解決這個問題。

4 結(jié)論

為了從微觀上理解糧倉效應(yīng)，我們針對二維隨機堆積建立了一個格點系統(tǒng)模型——pq模型。該模型的解析求解涉及到三對角形式的傳遞系數(shù)矩陣A(p,q)的相似對角化，這也是本文的重點。我們通過相似變換將A(p,q)對稱化，并且提出基于二階差分方程的方法比較有效地解決了該對稱矩陣的特征值和特征向量的計算。對于p=1/2、q=1這種典型的情況，我們解析地導(dǎo)出了糧倉底部的有效質(zhì)量隨顆?？傎|(zhì)量的變化關(guān)系；其它p，q取值的情況也可轉(zhuǎn)化為數(shù)值求解高次代數(shù)方程(25)的問題。本文討論的pq模型只是格點系統(tǒng)模型中一種可能的力傳遞形式，對于不同的堆積結(jié)構(gòu)和環(huán)境條件，可能具有不同的力傳遞規(guī)則，在數(shù)學(xué)上一般表現(xiàn)為不同的傳遞系數(shù)矩陣。本文提出的理論方法對于其它傳遞系數(shù)矩陣的情形有借鑒的價值，研究的結(jié)果也能加深對微觀力傳遞規(guī)則與宏觀壓力分布之間關(guān)系的認識。